Deployment optimization for a long-distance wireless backhaul network in industrial cyber physical systems

Abstract

Industrial wireless networks are an important component of industrial cyber physical systems, and their transmission performance directly determines the quality of the entire system. During deployment, the nodes of an industrial wireless network can be deployed in only some specific regions due to physical environment restrictions in the factory; thus, occlusions are not always effectively circumvented and network performance is reduced. Therefore, this article focuses on the layout problem of the industrial backhaul network: a WiFi long-distance, multi-hop network. The optimization objectives were network throughput and construction cost, and the network delay was used as a constraint. For small networks, we propose a hierarchical traversal method to obtain the optimal solution, whereas for a large network, we used a hierarchical heuristic method to obtain an approximate solution, and for extremely large networks, we used a parallel interactive local search algorithm based on dynamic programming. Then, if the original network layout cannot meet the transmission demands due to traffic bursts, we propose a network bandwidth recovery method based on the Steiner tree to recover the network’s performance. Finally, the results of a simulation showed that the algorithms proposed in this article obtain an effective solution and that the heuristic algorithm requires less computing time.

Keywords

Cyber physical systems industrial wireless network network optimization node deployment heuristic algorithm

Introduction

As wireless communication technology has developed, wireless technologies based on protocols such as WirelessHART,¹ WIA-PA,² and ISA100.11a³ have been widely applied in industrial cyber physical systems (CPS).⁴ However, because of restrictions to capital, technology, and equipment, among others, there may be a situation in which multiple network technologies coexist in a given period, thus forming a heterogeneous industrial wireless control network.

IEEE 802.11–based long-distance WiFi networks (WiFi-based long-distance (WiLD) multi-hop networks)⁵ can achieve long-distance signal transmission because the nodes incorporate high-power IEEE802.11 wireless cards and high-gain directive antennas. In addition, a single-hop link can be tens or hundreds of kilometers in length and feature high bandwidth, low cost, and wide coverage; therefore, it can provide backhaul access for the wireless network to achieve long-distance transmission of data.

This article considers the heterogeneous industrial monitoring network used in industrial CPS, as shown in Figure 1. The wireless network (such as WirelessHART, WIA-PA, or ZigBee) applied in each factory is considered a subnet of a heterogeneous network, and each subnet is connected via an industrial backhaul network^6,7 through IEEE 802.11 long-distance WiFi to constitute the wide-area heterogeneous industrial monitoring network. A control command given by the general control center is transmitted to the control plant in the subnet through the backhaul network and the subnet gateway, in that order. Meanwhile, the backhaul network provides an interconnection between the industrial control network and the Internet.

Figure 1.

A heterogeneous industrial monitoring network architecture in industrial CPS.

However, during the practical deployment of a network, WiLD nodes can be deployed in only specific regions due to physical environment restrictions in the factory. Meanwhile, a high tower cannot be erected in some factories, and an antenna can be erected only at the top of a building or existing lamp pole. As a result, the WiLD network’s antenna fails to effectively elude some occlusions, thus hindering the network’s transmission performance. Therefore, the layout of the network should be optimized. Figure 2 shows a partial map of Fushun Petrochemical Factory, China National Petroleum Corporation, in which the circled positions are metal tank areas or pipeline areas whose occlusion cannot be directly circumvented. In this case, it is necessary to consider the proper locations and number of nodes, the angle of the antenna, and other factors to optimize the network bandwidth and guarantee network performance.

Figure 2.

Example of occlusion in the factory.

This problem can be included in the scope of facility location. It is most similar to the P-median problem in facility location but with the following differences:

The traditional P-median problem is used to determine the minimum sum of the weights from demand points to the facility location, whereas in the network model of this article, consideration is given not only to the weights from demand points (gateway node) to the facility (exchange board) location but also to the weights between facilities through which the traffic flow passes.

Only a single target, namely, the weighted distance, is considered in the traditional P-median model, whereas the model in this article targets the optimization of the overall bandwidth and system overhead to plan the deployment of network nodes.

Thus, in this article, we used a network architecture based on WiLD multi-hop networks to optimize the node layout of the industrial CPS. The optimization objectives were network throughput and construction cost, and the network delay was used as a constraint. The main contributions of the article are as follows:

We propose a hierarchical P-median optimization model based on the traditional P-median model. There are two methods in this model. First, for a small network, a hierarchical traversal method based on the hierarchical P-median optimization model was adopted to obtain the optimal solution. Second, for a large network, a hierarchical heuristic method based on the hierarchical P-median optimization model was used to obtain an approximate solution.

To cope with the situation when the original network layout cannot meet the transmission demands due to traffic bursts, we propose a network bandwidth recovery method based on the Steiner tree (BRM-ST) to recover the network’s performance.

A test was carried out via simulation, and a comparison was made between the result obtained using the algorithm proposed in this article and that obtained using the traditional genetic algorithm (GA)–based heuristic method. The result showed that the proposed algorithm obtains an effective solution and that the heuristic algorithm requires less computing time.

The remainder of this article is organized as follows. We summarize the related works in section “Related work.” We formulate the problem and provide the system model in section “System model and formulation of the problem.” Following the model description, two methods based on hierarchical P-median optimization model and a parallel interactive local search algorithm based on dynamic programming (PILS-DP) are proposed in section “Solutions.” In section “Network node deployment for traffic bursts,” we propose network node deployment methods for traffic bursts. The experimental results are given in section “Experiment and analysis.” Finally, section “Conclusion” concludes the article.

Related work

Numerous studies of facility location models have been conducted by domestic and foreign scholars. Maximum Covering Location⁸ is often adopted when the decision-maker fails to satisfy all demands using the available resources (namely, fails to fully satisfy the demands under the given construction input and distance). Set Covering Location⁹ is used to obtain the lowest possible facility construction cost while satisfying all demands under a given emergency deadline. The P-median model¹⁰ determines the P facility locations and minimizes the weighted distance from the demand point to the facility (distance can also be expressed as traffic/transport time).

The traditional facility location problems include the p-center and p-coverage problems, which are NP-complete problems.¹¹ The existing research is more focused on the p-center problem. And some approximation algorithms such as branch-and-bound,¹² branch-and-price,¹³ Lagrangian relaxation,¹⁴ integral programming,¹⁵ linear programming relaxation,¹⁶ and randmized rounding¹⁷ are widely used to solve the p-center problem and no capacity limits facility location problems. Besides, the GA,¹⁸ neural network,¹⁹ and simulated annealing²⁰ are used to solve the large-scale p-center problem.

The current wireless network connection recovery methods are roughly divided into two approaches: deploying new relay nodes to the failure nodes or failure regions to recover the connection and moving some of the original nodes to recover the connection.

The method of moving some original nodes to recover the connection does not require the deployment of new nodes, but it does require the nodes to be mobile.

For recovery through deploying new nodes, some heuristic methods have been proposed. A previous study²¹ proposed a method in which the failure node’s neighboring positions are replaced to recover the network connectivity. However, this method cannot handle cases in which multiple nodes fail at the same time, and some nodes cannot move in some partitions.

In Lloyd and Xue,²² the minimum spanning tree algorithm based on a single-tiered relay node placement (MST-1TRNP) uses the Kruskal or the Prim algorithm to find the minimum spanning tree. Then relay nodes are deployed along the edge of the tree based on the length of the relay nodes’ communication radius. However, this method has higher complexity, and the average node degree is low after recovery. To improve the average node connectivity, the simply Connected Spiderweb (1C-SpiderWeb) algorithm is proposed.²³ 1C-SpiderWeb completes the relay node deployment by forming a network structure similar to a spider’s web. Although the average node connectivity increases, this approach requires a much higher number of relay nodes. A cell-based optimized relay node placement (CORP)²⁴ algorithm divides the network into units of equal size. Then, based on the grid positions in the partition, it deploys relay nodes to recover the connection. The network topology balances the network data transmission capacity and reduces the data transmission delay between partitions. However, the complexity of this approach is high, and it is not suitable for large-scale networks.

System model and formulation of the problem

The WiLD network is expressed as $G = (V, E)$ , where V and E refer to the sets of nodes and edges in the network, respectively, and $e_{ij} \in E$ denotes the link between node $v_{i} \in V$ and $v_{j} \in V$ . Suppose there is a finite number (N) of point locations at which to place WiLD nodes in the factory. The point locations capable of housing a WiLD node are denoted as $p_{i} (i = 1, 2, \dots, N)$ , and $p_{i} = {\begin{matrix} 1, place node \\ 0, others \end{matrix}$ . Of these locations, M must be assigned nodes that execute the functions of gateway nodes; specifically, $p_{i} = 1, i \in [1, M], M \leq N$ for these M nodes. In this article, assume that the traffic flow $f_{k} \in F$ of the perceptual and audio/video information collected by each sensor enters the WiLD network via a gateway node and then reaches the control center node $σ (σ \in V)$ after being transmitted by a WiLD relay node. Feedback information is sent to the relevant gateway node by the control center node $σ$ after data analysis and processing and then transmitted to the lower actuator node to execute the specific operation.

Each node in the network is equipped with one or several directive antennas. Here, we assume that the antennas are of the same type and have an identical beam width $θ$ in this article. To avoid interference between adjacent antennas, a buffer area with an included angle $Δ_{θ}$ is set between adjacent antennas of the same node, and the upper bound on the number of directive antennas that can be configured for each node is $k = ⌊ 360 \circ / (θ + Δ_{θ}) ⌋$ , as shown in Figure 3. Based on the conclusion in Zenghua et al.,²⁵ the included angle of antennas for the same node should not be less than $30 \circ$ . Thus, the process of determining the number and angles of directive antennas of a certain node is the process of determining the connectivity between this node and the surrounding nodes. We define $a_{ij} = {\begin{matrix} 1, \exists e_{ij} \in E \\ 0, others \end{matrix}$ , and the number of antennas on each node is $v_{i}$ , $a_{i} = \sum_{j \neq i \in [1, N]} a_{ij}$ , where $0 \leq a_{i} \leq ((2 π) / (θ + (π / 6)))$ . Suppose the basic construction cost of each node is the cost of the node without an antenna, $α$ , plus the cost added by each antenna, $β$ . Then, the construction cost of each node can be expressed as $c_{i} = α + β \cdot a_{i}$ . Therefore, the total cost of node construction in the network is $\sum_{i = 1}^{N} p_{i} c_{i}$ .

Figure 3.

Antenna configuration diagram.

Occlusion is measured by the bandwidth of the received signal. Suppose $B_{t}$ is the signal transmission bandwidth of the transmitting terminal and $B_{r}$ is the actual receiving bandwidth of the receiving terminal. Under the same transmission power, the larger the value of $B_{r}$ , the smaller the occlusion and the better the channel quality. Conversely, for a smaller value of $B_{r}$ , the occlusion will be more serious, and the channel quality will be poorer. When complete occlusion causes connection failure between nodes, $B_{r} = 0$ . Suppose that all links are symmetric links. The transmission bandwidth of signal in link $e_{ij} \in E$ is $B_{ij} : = B_{r}^{ij} = B_{r}^{ji}$ . $w^{f_{k}}$ is the end-to-end signal transmission bandwidth of the path that traffic $f_{k}$ goes through from source node $v_{i} (i \in [1, M])$ to destination node $σ$ . The traffic transmission bandwidth on the entire path is the minimum transmission bandwidth of a link on this path, namely, $w^{f_{k}} = min_{\forall e_{ij}^{f_{k}}, i, j \in [1, N]} w_{ij}^{f_{k}} \cdot x_{ij}$ , where $ϕ_{k} \leq w_{ij}^{f_{k}} \leq B_{ij}$ is the transmission rate of traffic $f_{k}$ on link $e_{ij} \in E$ and $ϕ_{k}$ is the minimum transmission bandwidth demand of traffic $f_{k}$ . $x_{ij} = 1$ indicates that node $v_{i}$ can be connected to node $v_{j}$ ; otherwise, $x_{ij} = 0$ . When $x_{ij} = 1$ , nodes $v_{i}$ and $v_{j}$ can be connected by link $e_{ij} \in E$ directly, or indirectly connected via other nodes.

This article aims to determine the placement locations of WiLD nodes and the number and directions of their antennas in the network to reduce the influence of occlusion on the network and maximize the network throughput while satisfying the traffic delay condition and bandwidth demand and minimize the overhead.

Optimization objective

The optimization objective is to obtain the maximum and minimum values, namely, $u_{1} = max \sum_{f_{k} \in F} w^{f_{k}}$ and $u_{2} = min \sum_{i = 1}^{N} p_{i} c_{i}$ , respectively. Therefore, we suppose $u (p, a) = (\sum_{f_{k} \in F} w^{f_{k}}) / (\sum_{i = 1}^{N} p_{i} c_{i})$ . Then, the optimization objective is $max u (p, a)$ .

Constraint

Real-time performance constraint

According to the analysis of Jintao et al.,²⁶ the delay of each hop is fixed when using the 2P MAC protocol, and we denote it as $τ$ . Then, the end-to-end delay of traffic $f_{k}$ can be expressed as $τ \cdot K_{f_{k}}$ , where $K_{f_{k}} = \sum_{i \neq j, i, j \in [1, N]} e_{ij}^{f_{k}}$ is the hop count of the path of $f_{k}$ . $e_{ij}^{f_{k}}$ is the link that traffic $f_{k}$ passes through. Suppose the end-to-end delay demand of traffic $f_{k}$ is $μ_{k}$ . Then $τ \cdot K_{f_{k}} \leq μ_{k}$ .

Bandwidth demand constraint

The bandwidth provided by the transmission path of the traffic must meet the bandwidth demand of the traffic, namely, $min_{\forall e_{ij}^{f_{k}}, i, j \in [1, N]} w_{ij}^{f_{k}} \cdot x_{ij} \geq d_{i}$ .

Link capacity constraint

The total transmission bandwidth occupied by each traffic on each link should not exceed the bandwidth provided by this link. Thus, $\forall i, j \in [1, N], j \neq i, \sum_{\forall f_{k} \in F} w_{ij}^{f_{k}} \leq B_{ij}$ .

Gateway node setting

For each gateway node M in the network, $p_{i} = 1, i \in [1, M]$ .

Node setting constraint

If node $v_{j}, j \in [1, N]$ can be connected with node $v_{i}, i \in [1, N]$ , then node $v_{i}$ must be the node that has been set. That is, $x_{ij} \leq p_{i}, x_{ij}, p_{ij} \in {0, 1}$ .

In conclusion, the optimization problem in this article can be considered a special type of the universal capacity limitation facility location (universal facility location) problem. The problem here differs from the universal facility location in that, in this article, traffic overhead is not only relevant to the node from the client side (gateway) to the location under selection, the overhead during the relay between exchange boards and from an exchange board to a data center node should also be considered. It has been proven that the universal facility location problem is a NP-hard problem; consequently, the problem in this article is also an NP-hard problem.

Solutions

Each directional antenna is replaced by an omni-directional antenna with the same transmission range. Point-to-point communication is adopted between nodes; therefore, the interference between nodes can be ignored. For each node, all nodes that can be connected to it when communicating at maximum power are considered to be nodes connected to this node. Thus, a topological graph of network connections can be constructed. The maximum bandwidth of an information transmission between two nodes connected with an edge is taken as the weight of that edge. According to the analysis of the network operation mode in the above section, the traffic flow bandwidth demand from the gateway node to the data center node is larger than the traffic flow bandwidth demand in the reverse direction. Therefore, if the bandwidth demand from a gateway node to the data center node can be fulfilled in the solution, the total bandwidth demand of the whole network is also satisfied.

The goal of the P-median model is to separately find proper locations for P facilities and designate for each demand point the given facility under the demand set of the given quantity and location. In addition, the P-median model identifies a set of candidate facility locations, thus obtaining the lowest transport cost between the facility and demand point. The traditional P-median model considers only the single-hop cost between each demand point and facility and is a single-hop single-objective model, whereas in the WiLD network, the end-to-end cost between two demand points should be considered, and it is therefore a multi-hop model. Both bandwidth demand and construction cost objectives are considered in this article; therefore, the proposed model is a multi-objective model. Consequently, to obtain a solution for the WiLD network, a hierarchical solution method based on the P-median model is presented in this article.

Hierarchical traversal solving method based on the P-median model

The total number of nodes in the network is N. Of these N nodes, M gateway nodes and a data center node must be deployed. Therefore, the total number of nodes in need of placement is $N - M - 1$ at most. The defined variable is $P \in [1, N - M - 1]$ , and the P-median model is established. If the total number of relay nodes cannot be determined, the improved P-median model can be applied to obtain the solution. The hierarchical traversal solving method based on the P-median model (P-HTM) is applied starting from $P = N - M - 1$ . The deployment location of each node is traversed under different P values to determine the maximum network flow under this deployment state until the maximum flow cannot meet the bandwidth demand; $u (p, a)$ is obtained at this time. Finally, values of $u (p, a)$ are compared under different P values to obtain $max u (p, a)$ of the network.

This algorithm is implemented by adding the restriction condition $p = \sum_{i = 1}^{N - M - 1} p_{i}$ on the uncertain node setting number and adopting a traversal method. The traversal method can obtain the exact solution, but its scouting speed is relatively slow. The P-median model-based relay-node location problem is a discrete problem; therefore, when the number of facility selection points is not large, the computing time of the traversal algorithm is limited.

The calculation of $u (p, a)$ under each P value in this algorithm is performed as shown in Algorithm 1. The basic steps of the algorithm are as follows:

Algorithm 1 Available bandwidth calculation
Require: $G = (V, E)$ , N, M
Ensure: $u (p, a)$
1: for ∀ gateway $v_{i} \in V$ do
2: connect $v_{s}$ to $v_{i}$ with links $e_{si}$ ;
3: $B_{si} \leftarrow + \infty$ ;
4: end for
5: for $e_{ij} \in E$ do
6: $w_{ij} \leftarrow 0$ ;
7: end for
8: while ∃ augmenting path $r_{i}$ from $v_{s}$ to $σ$ do
9: $r_{s σ} \leftarrow$ the path starting from $v_{s}$ to $σ$ and with minimum hop count;
10: $BwA (r_{s σ}, B_{ij})$ ;
11: for node from $σ$ to sdo
12: if ∃ bandwidth surplus point $v_{j}$ then
13: $r_{j σ} \leftarrow$ the path starting from $v_{j}$ to $σ$ and with minimum hop count;
14: $BwA (r_{j σ}, B_{ij})$ ;
15: end if
16: end for
17: calculate $\sum_{f_{k} \in F} w^{f_{k}}$ and $\sum_{i = 1}^{N} p_{i} c_{i}$ ;
18: end while

Initialization

The traffic bandwidth passing through each link is $0$ at the beginning. Set a virtual node $v_{s}$ in the network and add a link with bandwidth capacity $+ \infty$ between $v_{s}$ and each gateway node.

Construction of augmenting path

Judge whether there is an augmenting path between $v_{s}$ and the control center node $σ$ starting from $v_{s}$ . If there is such an augmenting path, find a path with the lowest number of intermediate nodes (the shortest augmenting path) from $v_{s}$ to $σ$ starting from $v_{s}$ along the link with the most residual bandwidth. When multiple paths have the same intermediate node number, select the augmenting path that can generate the largest augmenting capacity.

Link bandwidth occupation adjustment

Set $δ$ as the minimum value of residual bandwidth of each link on the selected augmenting path and increase $δ$ for bandwidth occupation of all links in the augmenting path.

Conduct reverse query for each node in each path along the selected augmenting path, in order, with $σ$ as the starting node. Suppose that node $v_{i}$ on this augmenting path satisfies the following conditions:

The input link of $v_{i}$ has residual bandwidth, but the output link has no residual bandwidth.

There is no such point with residual bandwidth at the input link but without residual bandwidth at the output link before $v_{i}$ .

There is residual bandwidth at the first section of the link on the augmenting path.

There is a path from $v_{i}$ to $σ$ .

Then, $v_{i}$ is defined as the surplus point. Keep the path from $v_{s}$ to $v_{i}$ in the original augmenting path unchanged, select the shortest path from $v_{i}$ to $σ$ along $v_{i}$ , and conduct link bandwidth occupation adjustment for path $v_{s} \to v_{i} \to σ$ again. The calculation function for link bandwidth is shown in Algorithm 2.

Algorithm 2 Function $BwA (r_{j σ}, c_{ij})$
Require: $G = (V, E)$ , $r_{s σ}$ , $B_{ij}$
Ensure: $w_{ij}$
1: for $e_{ij} \in E$ in path $r_{s σ}$ do
2: $δ \leftarrow min {B_{ij} - w_{ij}}$ ;
3: $w_{ij} \leftarrow w_{ij} + δ$ ;
4: end for

Regarding the calculation of $max u (p, a)$ , calculate the corresponding $u (p, a)$ under different P values and determine the maximum value $max u (p, a)$ .

Algorithm complexity analysis

During the execution process for the algorithm, each P should include searching the augmenting path and adjusting the link bandwidth occupation. When searching the augmenting path via depth-first search, if the last node in the current path can be accessed with link bandwidth extension, the longest length of the augmenting path is N; that is, node $σ$ is reached after advancing N steps to obtain the augmenting path. During reverse query, the number of reverse backoff steps is N at most. Therefore, the number of queries of the augmenting path is $| E |$ at most, and at most N nodes are queried each time. Thus, the total number of query operations is $N + | E | \cdot N$ at most and the complexity is $o (| E | \cdot N)$ . The complexity of conducting link bandwidth occupation adjustment in each iteration is $o (N)$ . Thereby, when the initial value of P is $N - M - 1$ , the overall complexity of the algorithm is $o (| E | \cdot N^{2} \cdot (N - M - 1))$ .

P-median model-based hierarchical heuristic solving method

To resolve the location problem of nodes when the network size is relatively large, a P-median model-based hierarchical heuristic solving method (P-HHM) is presented. First, let $P = P_{max} = N - M - 1$ and take the bandwidth that can be provided by the link as the weight to obtain the spanning tree that maximizes the total bandwidth, as shown in Algorithm 3. Then, subtract $(p_{max} - p)$ nodes in initial spanning tree based on different P values. Connect the residual nodes that are disconnected from the spanning tree to their adjacent nodes in the spanning tree. If there are several such adjacent nodes, select the nodes with the largest connection link bandwidth for connection. Calculate the maximum transmission bandwidth that the network can achieve under this P value until the obtained maximum transmission bandwidth fails to meet the bandwidth demand of traffic. Then, calculate $max u (p, a)$ under different P values to obtain the final result.

Algorithm 3 Spanning tree function $ST (G, W)$
Require: $G = (V, E)$ , $g_{i} \in V$
Ensure: minimum spanning tree edge set LE
1: ${a_{ij}^{k}} \leftarrow sort (E), v_{i}, v_{j} \in V, k \in [1, \| E \|]$ ;
2: $LE \leftarrow Φ$ ;
3: $LV \leftarrow Φ$ ;
4: while $V \neq Φ$ do
5: if $v_{i} \notin LV \cup v_{j} \notin LV$ then
6: $LE \leftarrow LE + max {a_{ij}^{k}}$ ;
7: $V \leftarrow V - {v_{i}} - {v_{j}}$ ;
8: $LV \leftarrow LV + {v_{i}} + {v_{j}}$ ;
9: end if
10: end while

The calculation process of the P-median-based heuristic solving method is as shown in Algorithm 4, and its main steps are as follows.

Algorithm 4 P-HHM
Require: $G = (V, E)$ , $g_{i} \in V$
Ensure: $max u (p, a)$
1: $LV \leftarrow ST (G, W)$ ;
2: ${ℜ_{i}} \leftarrow route (LV)$ ;
3: calculate $u (p, a)$ ;
4: while $P \geq 1$ do
5: if $w \geq \sum_{k} ϕ_{k}$ then
6: for $\forall v_{j} \in ℜ_{i}$ do
7: ${ℜ_{i}} \leftarrow {ℜ_{i}} - v_{j}$ ;
8: ${ℜ_{i}} \leftarrow route ({ℜ_{i}})$ ;
9: calculate $u (p, a)$ ;
10: end for
11: $P \leftarrow P - 1$ ;
12: end if
13: end while

Construction of spinning tree

The bandwidth that can be provided by the link is taken as the link weight based on the traditional “kruskal” algorithm to obtain the largest bandwidth spanning tree with data center node $σ$ as the root, to obtain the available path $ℜ_{i}$ from each gateway node $v_{i}, i \in [1, M]$ in G to the data center node $σ$ and to obtain $u (p, a)$ .

Deletion of nodes

Set $P = P - 1$ , subtract a node from the spanning tree of the first step, and then determine whether this should connect the branch node that was connected with the removed node to the node that has the largest bandwidth among all adjacent nodes ( $N - 1$ iterations are required for this process). Determine $u (p, a)$ in $N - 1$ conditions.

Bandwidth and overhead calculation

As the value of P decreases, $max \sum_{f_{k} \in F} w^{f_{k}}$ and $max \sum_{i = 1}^{N} p_{i} c_{i}$ will continue to decrease; when $max \sum_{f_{k} \in F} w^{f_{k}}$ reaches a certain threshold value, the algorithm ends to obtain $max u (p, a)$ .

Complexity analysis of the algorithm

In the construction process of the spanning tree based on the traditional “kruskal” algorithm, traversing is required for $| E |$ edges. The “While” loop of Algorithm 4 should be executed $N - M - 1$ time(s) at most. In each “While” loop, the nested “for” loop should be executed P time(s), and the maximum value of P is $N - M - 1$ . Therefore, the time complexity of the algorithm should not exceed $o ((N - M - 1)^{2})$ .

Parallel interactive local search algorithm based on dynamic programming

To resolve the location problem of nodes when the network size is relatively large, we propose a PILS-DP. The main ideas of this algorithm are as follows: the layout locations in the network are segmented according to geographic location; local searching is performed in each segment to obtain the optimal solution with boundary nodes under different conditions; and the globally optimal solution for the entire network is obtained by a dynamic programming algorithm based on the locally optimal solutions in each segment.

Network segmentation

Definition: cut set

Split the $N - M - 1$ nodes of the node set V in the network $G = (V, E)$ into m non-empty sets $V_{1}, V_{2}, \dots, V_{m}$ such that $V_{1} + V_{2} + \dots + V_{m} = V$ and $V_{i} \cap V_{j} = \emptyset, \forall i, j \in {1, 2, \dots, m}$ . For $\forall V_{i} \subset V$ , if $\exists e_{ab} \in E$ and $v_{a} \in V_{i}$ such that $v_{b} \in V - V_{i}$ , then we call the network node set $V_{i}$ a cut set of V.

Take the minimum hop count $h_{min} = min {h_{i, σ}}, i \in [1, M]$ from each gateway node $v_{i}$ to the data center node $σ$ . Divide all nodes except the data center node $σ$ evenly into $h_{min}$ cut sets and place all gateway nodes in the same cut set.

Local search in each cut set

For each cut-set $V_{i}$ , we use the depth-first search method to find any feasible solution that can meet the constraints first and obtain the node settings in this situation. For the nodes in each cut set, we define three operations: add, remove, and exchange, and the computation rules are shown below. If there is a feasible solution ${\tilde{V}}_{i} \in Π ({\hat{V}}_{i})$ that can satisfy $u_{{\tilde{V}}_{i}} > u_{{\hat{V}}_{i}}$ , then update ${\hat{V}}_{i} : = {\tilde{V}}_{i}$ , and

\begin{matrix} Π ({\hat{V}}_{i}) : = {\hat{V} + {v_{a}}, v_{a} \in V - \hat{V}} \\ \cup {\hat{V} - {v_{b}}, v_{b} \in \hat{V}} \\ \cup {\hat{V} + {v_{a}} - {v_{b}}, v_{a} \in V - \hat{V}, v_{b} \in \hat{V}} \end{matrix}

Definition: the boundary nodes

In the cut set $V_{i}$ of the network, if $\exists e_{ab} \in E$ such that $v_{a} \in V_{i}$ and $v_{b} \in V - V_{i}$ , we call $v_{a}$ the boundary node of the cut set $V_{i}$ .

In the process of local search, record the optimal solution in each situation that contains boundary nodes as the state variable in the feasible solution space during the period of inter-cut-set search.

Algorithm 5 PILS-DP
Require: $G = (V, E)$ , $g_{i} \in V$ , N, M
Ensure: state variable of each $V_{i}$
1: $h_{min} \leftarrow min {h_{i, σ}}, i \in [1, M]$ ;
2: divide N into $h_{min}$ cut set evenly and all gateways are in one cut set;
3: for ∀ cut set $V_{i}$ do
4: get a feasible deployment solution $\hat{V}$ using depth-first search;
5: $Π ({\hat{V}}_{i}) \leftarrow {\hat{V} + {v_{a}}, v_{a} \in V - \hat{V}} \cup {\hat{V} - {v_{b}}, v_{b} \in \hat{V}} \cup {\hat{V} + {v_{a}} - {v_{b}}, v_{a} \in V - \hat{V}, v_{b} \in \hat{V}}$ ;
6: if $\exists {\tilde{V}}_{i} \in Π ({\hat{V}}_{i})$ , $u_{{\tilde{V}}_{i}} > u_{{\hat{V}}_{i}}$ then
7: ${\hat{V}}_{i} \leftarrow {\tilde{V}}_{i}$ ;
8: end if
9: record the optimal solution in each situation with edge node as the state variable in each set cut;
10: end for

Dynamic programming inter-cut-set search

The basic idea to solve the problem with dynamic programming is as follows.

The $h_{min}$ cut sets obtained from the network segment make up the subproblem of $h_{min}$ phases. Each deployment decision of the phase is the subproblem solution in this phase.

In each phase, the feasible solutions are obtained by local search in the cut set for the different states (set to k, where $k = 1, 2, 3, \dots$ ). The state variable is denoted by $s_{k}$ , and the value of is the k possible alternative feasible solution.

In each state $s_{k}$ , the decision-making variable is $ξ_{k} (s_{k})$ . The values of $ξ_{k} (s_{k})$ are the elements in the state variable set.

The index function of the cost measured in the $h_{i} (i = 1, 2, \dots, h_{min})$ phase is $v_{k, n}$ . $v_{k, n} (s_{k}, x_{k})$ is the accumulation function from state $s_{k}$ to the decision-making variable $x_{k} (s_{k})$ and all subsequent stages.

The optimal value of the index function is the optimal index function, denoted by $f_{k} (s_{k})$ , and $f_{k} (s_{k}) = OPT {v_{k} (s_{k}, ξ_{k}) + f_{k + 1} (s_{k + 1})}$ . This is the basis for decision-making under the state $s_{k}$ .

According to the overall optimal solution and the reverse solving ideas of dynamic programming, we get the decision $ξ_{k} (s_{k})$ of each phase, and we can then obtain the ordered overall optimal decision $p_{1, h_{min}}$ , that is, $p_{1, h_{min}} = {ξ_{1} (s_{1}), ξ_{2} (s_{2}), \dots, ξ_{n} (s_{n})}$ .

The complexity of the algorithm

A depth-first search is carried out in each cut set. In the worst case, it is not until the last node that the feasible solution is obtained when using the depth-first search method to search for the feasible solution satisfying the constraints. Thus, the total computational complexity of the search is $o ((N - M - 1) / (h_{min})^{2})$ .

The time complexity of the dynamic programming algorithm is $o (\sum_{i \in [1, h_{min}]} k_{i} h_{i})$ .²⁷ Because $\sum_{i \in [1, h_{min}]} k_{i} h_{i} \leq h_{min} \sum_{i \in [1, h_{min}]} k_{i} \leq h_{min} (N - M - 1)$ , and in the worst cases, $h_{min} = 1$ , the overall complexity of PILS-DP is $o ((N - M - 1)^{2})$ .

Network node deployment for traffic bursts

When the original network layout cannot meet the transmission demands due to traffic bursts, we can recover the network’s performance by increasing the number of nodes or the number of antennas in some nodes to increase the network bandwidth and meet the traffic transmission demands. Thus, we propose a network BRM-ST. The basic idea of this method is that when the bandwidth is insufficient, we use the residual bandwidth that another path can provide to transmit; if we cannot find another path with sufficient residual bandwidth to meet the demands, we consider the link that cannot meet the bandwidth demands to be disconnected, and the network is partitioned. Then, we structure the network partition connectivity based on the Steiner tree²⁸ to meet the bandwidth transmission demands.

One possible approach is to select the node with the minimal cost path to the Steiner tree and add the node to the tree through the minimal cost path. However, the time complexity of iterating to find the possible paths is too high. To reduce the time complexity, we use the basic idea of the progressive increase in the Dijkstra algorithm to define the node search order and skip the feasible paths with higher cost to achieve a fast search. The algorithm mainly includes the following three stages: preparation, Steiner tree construction, and fast searching.

Preparation

The locations where nodes can be placed are assumed to have nodes equipped with omni-directional antennas and maximum power. When the transmission bandwidth of a link is less than the demand of the traffic burst, the link is seen as disconnected, and the network is partitioned. All partitions segi form the set D, ordered from large to small order according to distance ( $Dist (segi, σ)$ ) from the data center node $σ$ . All the locations that can accept a node between each partition form the set Q, and the number of such locations is n.

Steiner tree construction

Define the path weight of the path $P_{i}$ as $g_{i} = λ / w_{i}$ , where $w_{i}$ is the path bandwidth, whose value is the minimum of all the link bandwidths in this path; that is, $w_{i} = min_{\forall e_{ij} \in P_{i} 0} {w_{ij}}$ , and $λ$ is a constant.

Using the Dijkstra shortest path tree (Dijkstra-SPT) algorithm,¹⁰ we can get the minimum path weighted Steiner tree T rooted in $σ$ for the current network connection.

Fast searching

In the fast-searching stage, we introduce two variables $δ$ and $π$ for each node in the network. $δ (a)$ records the minimum-cost path from node $v_{a}$ to T, and $π (a, δ)$ records the parent node of $v_{a}$ in the minimum-cost path tree. The search starts from the data center node $δ$ and proceeds to gateway nodes that have traffic burst demands.

The minimum cost from node $v_{a}$ to the Steiner tree is defined as follows: for a given Steiner tree T rooted in node $δ$ , the cost of node $v_{a}$ reaching T through the path $P (a, T)$ is the cost without the shared link cost, that is, $cost (a, P) = \sum_{e \in P, e \notin T} cost (e)$ . Its value represents the added value to T when connecting to T through the path P.

The path weight set from $v_{a}$ to T for weights less than or equal to $λ / d_{a}$ is denoted as $S_{P} (a, T) = {P | g (P) \leq (λ / d_{a})}$ . In addition, the path cost from $v_{a}$ to T less than or equal to $λ / d_{a}$ is $δ (a) = min {cost (a, P) | P \in S_{P} (a, T)}$ . Thus, we can get the minimum feasible cost from $v_{a}$ to T, that is $δ (a, T) = min {δ (a)}$ , and the corresponding path is the minimal cost path from $v_{a}$ to T.

The recursion relation of $δ (a)$ is

δ (v_{a}) = {\begin{matrix} 0, v_{a} = σ, g_{i} = 0 \\ min {\begin{matrix} δ (w, g_{i} - g (w, v_{a})) \\ + cost (w, v_{a}) \\ | w \in adjacency (v_{a}) \end{matrix}}, v_{a} \neq s, g_{i} > 0 \\ \infty, others \end{matrix}

In the fast-searching stage, we search the node with minimum $δ$ (which has the minimum-cost path to the current tree) in turn. If the node is the gateway node, the node joins the path tree through the minimum-cost path. The main idea of the fast-searching stage is as follows.

Initialization. The current time delay constraint, Steiner tree $T = \emptyset$ , $δ (σ, 0) = 0$ , $δ (a, g_{i}) = \infty$ , $v_{a} \neq σ$ , $π (a, σ) = - 1$ (line 1 in Algorithm 6).

Algorithm 6 Fast searching
Require:Q, ${g_{i} \geq 0}$
Ensure: T
1: $T \leftarrow \emptyset$ , $δ (v_{a}) \leftarrow \infty$ , $π (v_{a}, g_{i}) \leftarrow - 1$ ;
2: choose $v_{a} \in Q$ with the $min cost (v_{a}, T)$ ;
3: if $v_{a}$ is the gateway node with growing bandwidth demand then
4: $v_{b} \leftarrow v_{a}$ , $k \leftarrow g_{i}$ ;
5: while $v_{b} \notin T$ do
6: $T \leftarrow T + {v_{b}}$ ;
7: $δ (v_{b}, k) \leftarrow 0$ ;
8: if $v_{b}$ has been searched then
9: $Q \leftarrow Q + {v_{b}}$ ;
10: $v_{c} \leftarrow π (v_{b}, k), k \leftarrow k - g (v_{c}, v_{b}), v_{b} \leftarrow v_{c}$ ;
11: end if
12: end while
13: end if
14: if $v_{a} = σ$ then
15: $T \leftarrow T + {v_{a}}$ ;
16: end if
17: for $\forall v_{d} \in adjacent (v_{a})$ do
18: $l \leftarrow g_{i} + g (v_{a}, v_{d})$ ;
19: if $v_{d} \notin T & δ (v_{d}, l) > δ (v_{a}, g_{i}) + cost (v_{a}, v_{d})$ then
20: $δ (v_{d}, l) \leftarrow δ (v_{a}, g_{i}) + cost (v_{a}, v_{d})$ ;
21: $π (v_{d}, l) = v_{a}$ ;
22: end if
23: if $δ (v_{d}, l)$ is the minimum deployment cost from $v_{d}$ to T & l is the minimum feasible path weight then
24: $Q \leftarrow Q + {v_{d}}$ ;
25: end if
26: end for
27: $Q \leftarrow Q - {v_{a}}$ ;

Choose the node with the minimum cost to T from Q. If multiple nodes are optimal, choose the node with maximum path bandwidth; denote it as $v_{a}$ . The minimum feasible path weight is $g_{i}$ (line 2 in Algorithm 6).

If $v_{a}$ is the gateway node, let $v_{b} = v_{a}$ , $k = g_{i}$ . If $v_{b} \notin T$ , and the corresponding edge $e (b, π (b, k))$ will join T, and let $δ (v_{b}, k) = 0$ . If $v_{b}$ has been searched, then $v_{b}$ will rejoin Q, along with $v_{c} = π (b, k), k = k - g (c, b), v_{b} = v_{c}$ and repeat until $v_{b} \in T$ (lines 3–16 in Algorithm 6).

For each neighbor node $v_{d}$ of $v_{a}$ (denoted as $v_{d} \in adjacent (v_{a})$ ), let $l = g_{i} + g (a, d)$ . If $v_{d} \notin T, δ (d, l) > δ (a, g_{i}) + cost (a, d)$ , update the cost information from $v_{d}$ to T, and let $δ (v_{d}, l) = δ (v_{a}, g_{i}) + cost (v_{a}, v_{d}), π (d, l) = v_{a}$ . If $v_{d}$ has been searched, $δ (d, l)$ is the minimum deployment cost and l is the minimum feasible path weight, then $v_{d}$ will rejoin Q (lines 17–26 in Algorithm 6).

Then, delete $v_{a}$ from the set Q (line 27 in Algorithm 6). If multiple nodes have traffic burst demands, repeat the above steps until T covers all these nodes and the algorithm ends.

Finally, check the results and delete the duplicate nodes and antennas that can communicate in the feasible bandwidth to further reduce the number of node deployments.

The complexity of the algorithm

In the Steiner tree construction stage, we use the Dijkstra-SPT algorithm to establish the spanning tree. Thus, the complexity is $o (n^{2})$ . In the fast-searching stage, we choose the most appropriate nodes from the node set for the next extension, which is to find the node with the minimum distance from the vector. The complexity depends on the data structure of the selected algorithms. In this article, we use an array to implement it; therefore, the time complexity is $o (n)$ . In the process of adding the node through the minimal cost path to Steiner tree T, the time complexity is $o (l_{a})$ , where $l_{a}$ is the number of nodes in the path. In the process of the search loops, the total number of loops is the total number of nodes to search. Although there are some nodes that are searched again, previous experiments showed that the average search frequency of a node is less than $2$ , so the time complexity of the loops is $o (n)$ . Overall, the time complexity of BRM-ST is $o (n^{2})$ , where n is the number of remaining locations where the nodes can be deployed.

Experiment and analysis

This section describes a simulation experiment conducted to verify the performance of the algorithm proposed in this article. We used MATLAB to simulate the performance of our algorithms. A network is randomly generated in a 2 km × 2 km area, points in the network are randomly selected for the data center node and the gateway node, and traffic is sent to the data center node via the gateway node. The maximum transmission bandwidth and minimum network construction overhead of the network under different conditions are obtained by changing the maximum number of selectable point locations, the transmission power of the nodes, and the traffic bandwidth demand in the network. The maximum transmission bandwidth of traffic between nodes is set to 30 M. The bandwidth between nodes decreases linearly as the distance between the nodes increases. Because there is no similar algorithm that considers the same two objectives (network throughput and construction cost), we compare our algorithms’ results with the existing GA. We simulate the computation times, the maximum bandwidth provided under different numbers of selectable point locations, and the number of nodes required under different traffic bandwidth transmission demands and the maximum transmission bandwidth with these algorithms. In addition, a comparison is made between the P-HTM and P-HHM algorithms presented in this article and GA. Because as we know that this is the first article that considers two objects of nodes layout with IEEE 802.11 long-distance 2P MAC, there are no techniques available in the literature, and we just compare our methods with the traditional GA. The test results are shown in Figures 4 and 5.

Figure 4.

Computing time of each algorithm under different optional locations.

Figure 5.

Calculation results of P-HTM, P-HHM, PILS-DP, and GAs. (a) The maximum bandwidth provided under different optional locations. (b) The number of nodes required under different traffic bandwidth demands. (c) The number of nodes required under the different maximum bandwidth of nodes (traffic bandwidth demand is 120).

Figure 4 shows the computing time of each algorithm for solving the problem based on different numbers of selectable point locations. It can be observed from the figure that the computing times of the P-HHM and PILS-DP algorithms are lower than that of GA and that of the P-HTM is higher than GA. In addition, as the number of selectable point locations increases, the computing times of the P-HTM algorithm and GA algorithm increase rapidly, but that of the P-HHM and PILS-DP algorithms increase slowly. Figure 5(a) shows the maximum bandwidth that the network can provide under different numbers of selectable point locations. Figure 5(b) shows the number of nodes required for network transmission under different traffic bandwidth transmission demands. Figure 5(c) shows the total number of nodes required by network with different maximum transmission bandwidths that can be provided by each node when the transmission bandwidth demand of network traffic is 120. It can be seen from Figure 5 that the results from the P-HTM, P-HHM, PILS-DP, and GAs are basically identical. Therefore, the simulation result shows that an exact solution can be obtained via P-HHM and PILS-DP algorithms, which also take less computing time. While with the parallel computing feature, PILS-DP can further reduce the computing time.

To further verify the performance of the algorithm, a network is set up in Fushun Olefin Factory, Sinopec Group for a test, as shown in Figure 6. Plants and equipment in the factory are distributed within the scope of 1.5 km × 1.5 km, and occlusions include buildings, metal tanks, and pipelines in the factory. There are 14 points in the factory that can be selected to set up nodes, and two of the nodes should be connected to gateway nodes.

Figure 6.

Actual node deployment in the factory.

A test is conducted under different bandwidth demands of two gateway nodes, and the required number of nodes is obtained by applying different algorithms, as shown in Figure 7. Figure 7 shows that the solution of the network model can be obtained via both the P-HTM and P-HHM algorithms.

Figure 7.

Node number required under different traffic bandwidth demands in the factory (unit: M).

In transmission performance recovery experiments, we compare our proposed algorithm based on a Steiner tree with the traditional minimum spanning tree–based algorithm in terms of the number of relay nodes deployed and the average node degree when meeting the transmission demands, as shown in Figures 8 –10.

Figure 8.

The impact of the partition number on the number of relay nodes deployed.

Figure 9.

The impact of the communication radius on the number of relay nodes required.

Figure 10.

The impact of the size of the partition on the number of relay nodes required.

Figure 8 describes the impact on the number of relay nodes’ deployment of the partition number. We can see from the figure that with the increase in the number of partitions, the required number of relay nodes of the two algorithms increases. This trend corresponds to the actual situation. An increase in the number of partitions will lead to the use of more relay nodes when connectivity is recovered. However, we can also see clearly from the figure that the proposed BRM-ST algorithm requires fewer relay nodes than the traditional minimum spanning tree–based algorithm, and with an increase in the number of partitions, this advantage will be more obvious. This phenomenon occurs because with an increase in the number of partitions, the Steiner tree can connect a much higher number of partitions, which will reduce the number of relay nodes required.

Figure 9 describes the impact of the communication radius on the number of relay nodes required when the partition number is 9. We can see from the figure that as the communication radius increases, the number of relay nodes required is reduced. This trend also corresponds with the actual situation because with the increase in the communication radius, the number of relay nodes required will decrease when the path length is fixed. At the same time, we can see from the figure that the number of relay nodes required by the BRM-ST algorithm proposed in this article is less than that required by the algorithm based on the minimum spanning tree. Of course, this advantage is not obvious because the increasing communication radius will lead to a decrease in the path length.

Figure 10 describes the impact of the size of the partition on the number of relay nodes required when the partition number is 9 and the communication radius is 20 km. We can see from the figure that with an increase in partition size, the number of relay nodes required increases. This result occurs because the increase in partition size will lead to an increase in the distance between partitions, thus the number of relay nodes required will increase. We can also see clearly from the figure that the number of relay nodes required by the proposed algorithm BRM-ST is less than that by the algorithm based on the minimum spanning tree, and with the increase in the partition size, these advantages are more obvious. This result is because the sum of the distances required by BRM-ST is smaller than that by the minimum spanning tree, and with the increase in the partition size, this advantage will be more obvious.

Conclusion

The goal of this article is to solve the layout problem of the P-median model for a heterogeneous industrial monitoring network based on a WiLD multi-hop network architecture. The traditional P-median model cannot be applied to optimize a multi-hop network or a multi-objective network; thus, a hierarchical P-median optimization method based on the traditional P-median model was proposed in this article. The network node layout was optimized by taking the network throughput and construction cost as optimization objectives and using network delay as a constraint. To resolve the node location problem when the network size is relatively large, a PILS-DP is proposed. Then, when the original layout of the network cannot meet the transmission demands due to traffic bursts, we recover the network’s performance using a network BRM-ST. A simulation experiment was performed using a network constructed in a practical environment. Furthermore, we compared the solution obtained by the algorithms proposed in this article with that obtained by the traditional GA-based heuristic method, thus verifying the performance of the proposed algorithms.

Footnotes

Acknowledgements

The authors thank the 12th International Conference on Wireless Algorithms, Systems, and Applications (WASA 2017) for providing an opportunity to preliminarily present the deployment problem and its preliminary solution. The authors also thank the anonymous reviewers for their constructive advice.

Handling Editor: Suparna De

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (grant nos 61502474, 61501447, and 61233007) and Youth Innovation Promotion Association of the Chinese Academy of Sciences.

References

Zhang

Feng

Weight-based link scheduling for convergecast in WirelessHART network. Int J Distrib Sens N 2016; 12(7): 1–6.

Wang

Zhang

Wang

et al . A TR069 WAN management protocol for WIA-PA wireless sensor networks. In: 25th Wireless and Optical Communication Conference (WOCC), Chengdu, China, 21–23 May 2016, pp.1–4. New York: IEEE.

Serizawa

Fujiwara

Yano

et al . Reliable wireless communication technology of adaptive channel diversity (ACD) method based on ISA100.11a Standard. IEEE T Ind Electron 2017; 64(1): 624–632.

Hussain

Ahmed

Sarma

et al . An efficient TDMA MAC protocol for multi-hop WiFi-based long distance networks. Wireless Pers Commun 2016; 86(4): 1971–1994.

Gai

Qiu

Zhao

et al . Dynamic energy-aware cloudlet-based mobile cloud computing model for green computing. J Netw Comput Appl 2016; 59: 46–54.

Kong

Yuan

Xiao

et al . Stochastic optimization model and solution algorithm for blending procedure of process industrial. Inform Control 2016; 45(1): 40–44.

Gai

Qiu

Ming

et al . Spoofing-jamming attack strategy using optimal power distributions in wireless smart grid networks. IEEE T Smart Grid 2017; 8: 2431–2439.

Yinhong

. The algorithms for the bi-level programming location model based on the demand assigning. In: 10th International conference on fuzzy systems and knowledge discovery (FSKD), Shenyang, China, 23–25 July 2013, pp.999–1004. New York: IEEE.

Diaby

Karwan

MH.

Advances in combinatorial optimization: linear programming formulations of the traveling salesman and other hard combinatorial optimization problems. Singapore: World Scientific, 2016.

10.

Colmenar

Greistorfer

Martł

et al . Advanced greedy randomized adaptive search procedure for the obnoxious p-median problem. Eur J Oper Res 2016; 252(2): 432–442.

11.

Megiddo

Supowit

KJ.

On the complexity of some common geometric location problems. SIAM J Comput 1984; 13(1): 182–196.

12.

Daskin

Maass

KL.

The p-median problem. In: Daskin

Maass

(eds). Location science. New York: Springer, 2015, pp.21–45.

13.

Ceselli

Righini

A branch-and-price algorithm for the capacitated p-median problem. Networks 2005; 45(3): 125–142.

14.

Kochetov

Large neighborhood local search for the p-median problem. Yugosl J Oper Res 2015; 15(1): 53–63.

15.

Fast

Hicks

IV.

A branch decomposition algorithm for the p-median problem. INFORMS J Comput 2017; 29(3): 474–488.

16.

Gollowitzer

Gouveia

Laporte

et al . A comparison of several models for the hamiltonian p-median problem. Networks 2014; 63(4): 350–363.

17.

Hajiaghayi

et al . A constant factor approximation algorithm for fault-tolerant k-median. ACM T Algorithms 2016; 12(3): 1–20.

18.

Drezner

Brimberg

Mladenovic

et al . New heuristic algorithms for solving the planar p-median problem. Comput Oper Res 2015; 62: 296–304.

19.

Dominguez

Munoz

PJ.

An efficient neural network algorithm for the p-median problem. Lect Notes Comput Sc 2002; 2527: 460–469.

20.

Levanova

Loresh

MA.

Algorithms of ant system and simulated annealing for the p-median problem. Automat Rem Contr+ 2004; 65(3): 431–438.

21.

Lee

Younis

Lee

Connectivity restoration in a partitioned wireless sensor network with assured fault tolerance. Ad Hoc Netw 2015; 24: 1–19.

22.

Lloyd

Xue

Relay node placement in wireless sensor networks. IEEE T Comput 2007; 56(1): 134–138.

23.

Senel

Younis

Akkaya

. A robust relay node placement heuristic for structurally damaged wireless sensor networks. In: IEEE 34th conference on local computer networks, Zurich, 20–23 October 2009, pp.633–640. New York: IEEE.

24.

Lee

Younis

Optimized relay placement to federate segments in wireless sensor networks. IEEE J Sel Area Comm 2010; 28(5): 742–752.

25.

Zenghua

Nan

Zhibin

et al . Long-distance IEEE 802.11 wireless mesh networks. Chin J Comput 2012; 35: 1209–1222.

26.

Jintao

Peng

et al . End-to-end delay analysis in wide-area heterogeneous wireless network for industrial monitoring and control applications. Sci China: Inform Sci 2015; 10: 1249–1262.

27.

Chow

Tsitsiklis

JN.

The complexity of dynamic programming. J Complexity 1989; 5(4): 466–488.

28.

Brazil

Graham

Thomas

et al . On the history of the Euclidean Steiner tree problem. Arch Hist Exact Sci 2014; 68(3): 327–354.