extfieldsofstudy:
Computer Science
Mathematics
Engineering

text:
Hybrid Radio/Free-Space Optical Design for Next Generation Backhaul Systems

The deluge of date rate in today's networks poses a cost burden on the backhaul network design. Designing cost efficient backhaul solutions becomes an interesting, yet challenging, problem. Traditional technologies for backhaul networks include either radio-frequency backhauls (RF) or optical fibers (OF). While RF is a cost-effective solution as compared to OF, it supports lower data rate requirements. Another promising backhaul solution that may combine both a high data rate and a relatively low cost is the free-space optics (FSO). FSO, however, is sensitive to nature conditions (e.g. rain, fog, line-of-sight, etc.). A more reliable alternative is therefore to combine RF and FSO solutions through a hybrid structure called hybrid RF/FSO. Consider a backhaul network, where the base-stations (BS) can be connected to each other either via OF or via hybrid RF/FSO backhaul links. The paper addresses the problem of minimizing the cost of backhaul planning under connectivity and data rates constraints, so as to choose the appropriate cost-effective backhaul type between BSs (i.e., either OF or hybrid RF/FSO). The paper solves the problem using graph theory techniques by introducing the corresponding planning graph. It shows that under a specified realistic assumption about the cost of OF and hybrid RF/FSO links, the problem is equivalent to a maximum weight clique problem, which can can be solved with moderate complexity. Simulation results show that our proposed solution shows a close-to-optimal performance, especially for practical prices of the hybrid RF/FSO.

I. INTRODUCTION
Cellular networks, flooded by a large demand for mobile data services, are expected to undergo fundamental transformation. In order to significantly increase the data capacity, coverage performance, and energy efficiency, the next generation mobile networks (5G) [1] are expected to move from the traditional single, high-powered base-station (BS) to the deployments of multiple overlaying access points of diverse sizes (i.e. microcell, picocell, femotcell, etc.) using different radio access technologies. The resulting network system architecture is refereed to as heterogeneous networks (HetNets). Inter and intra-cell interference represents a major cause of the degradation of the HetNets performance. A successful interference mitigation architecture is the so-called heterogeneous cloud radio access network architecture obtained by connecting BSs from different tiers to a cloud (central processor) through backhaul links. The cloud largely improves the performance of HetNets; however, it necessitates a considerable amount of backhaul communications in order to share the data streams between all BSs across the network. Giving that the links are capacity limited, upgrading the backhaul and increasing its capacity to support the tremendous amount of data is a necessity [2].
Traditional technologies for the backhaul network design include copper, microwave radio links (RF), and optical fibers (OF). The price of copper lines, the most widely used backhaul technology [3], is growing linearly with the provided capacity. For high data rates, copper lines become expensive. Microwave radio technologies are also used for the backhaul network design. However, they represents only 6% of the total used transport media [3] due to the initial investment in the licensed part of the spectrum [4]. Moreover, low frequency radio is limited in terms of data-rates due to interference problems and high frequency radio are limited in the transmission coverage area.
Optical fiber (OF) backhaul links provide high data-rates over long distances. However, they are expensive to be deployed and require considerable initial investment. Recently, the free-space optics technology (FSO) becomes an interesting substitute [5] for the next generation cellular backhaul networks design. Generally, an FSO link refers to a laser beam between a pair of photo-detector transceivers using the freespace as medium of transportation. Giving that its wavelength is in the micrometer range, which is an unlicensed band, FSO links not only are free to use but they are also immune to electromagnetic interference generated by the radio-frequency (RF) links. The high bandwidth and interference immunity features make an FSO link up to 25 times more effective than a RF link in terms of capacity [6]. In addition, FSO represents a cost effective solution compared to OF.
Unlike OF links that are always reliable, FSO links are less reliable since they are sensitive to weather conditions, such as fog, snow, and rain [7]. Therefore, reliability should be taken into account when designing FSO-based backhaul networks. In order to cope with the varying reliability and combine the advantages of RF (reliability) and FSO (capacity), the hybrid RF/FSO technology has been proposed. Hybrid RF/FSO transmits, when possible, simultaneously on both the RF and FSO links. In harsh weather conditions, in which the FSO link is affected, the data are transmitted solely on the RF link [8]. Moreover, hybrid RF/FSO transceivers can be quickly deployed over several kilometers [9] and can also be easily combined with OF [10]- [12]. All these benefits make hybrid RF/FSO a suitable complementary option for upgrading the existing backhaul network [7], [13].
In the past few years, hybrid RF/FSO has attracted a large amount of research. Most of the existing work focused on the determination of the factors affecting the FSO link performance and finding solutions to improve the quality [13], [14]. However, fundamental problems of hybrid RF/FSO architecture optimization for the backhaul network topology design are only at their beginning.
The authors in [15] designed an efficient and scalable algorithm to optimize a given physical layer objective for 2 and 3 optical transceivers per node with a minimum number of links. Kashyap et al. [6] designed a routing algorithm for hybrid RF/FSO networks that backs up the traffic to the FSO routes when it could not be carried by the RF links. In [16], the authors considered an hybrid RF/FSO network in which the RF and FSO links operate at different data rates. They derived an upper bound for the capacity per node that is asymptotically achievable for random networks. Rak et al. [8] introduced a linear integer programming model to determine routing in hybrid RF/FSO network in which the FSO link availability is varying with the weather conditions.
Numerous mixed integer programming model have been proposed to formulate the problem of backhaul network design using hybrid RF/FSO technology. In particular, Son et al. [17] presented an algebraic connectivity-based formulation for the design of the backbone of wireless mesh networks with FSO links and solved it using a greedy approach that iteratively inserts nodes to maximize the algebraic connectivity. The authors in [18] proposed to maximize the network throughput by installing as many FSO links as possible under the constraint that the number of FSO links in a node is bounded. In [19], Ahdi et al. introduced a mixed integer programming model to find the optimal placement of FSO links in order to upgrade an existing RF backhaul network. Similarly, reference [20] propose to upgrade an existing RF backhaul network with FSO links using the minimum number of FSO links to guarantee a target network throughput when RF links are non available due to interference. The authors in [11] considered the upgrade of a pre-deployed OF backhaul network using FSO links and mirrors for nodes not in line-of-sight of each other. For two link-disjoint paths networks, they formulated the problem as a mixed integer programming and extended the study in [21] to K link-disjoint paths. In [12], the same group of authors studied the impact of the parameter K on the optimal solution.
In this paper, we consider the problem of minimizing the network deployment cost under target rate constraints, so as to determine which type of connections exists between each pair of BSs, i.e., OF, hybrid RF/FSO, or none. A major concern is to guarantee network connectivity, which is achieved by connecting, each pair of nodes in the network, possibly via multiple hops. While the deployment cost of hybrid RF/FSO links depends mainly on the cost of the hybrid RF/FSO transceivers, the deployment cost of OF links depends mostly on the distance between the two end nodes. On the other hand, OF links always satisfy the data rate constraint. The performance of hybrid RF/FSO, however, degrades with the distance and the number of installed links. The paper solves the problem using graph theory techniques by introducing the corresponding planning graph. The paper main contribution is to provide a close to optimal explicit solution to the problem. The paper shows that under a specified realistic assumption about the cost of OF and hybrid RF/FSO links, the problem can be reformulated as a maximum weight clique problem, which can be globally solved using efficient algorithms [22].
The rest of this paper is organizes as follows: Section II presents the considered system model and the problem formulation. Section III illustrates the proposed solution. Before concluding in Section V, simulation results are presented in Section IV.

II. SYSTEM MODEL AND PROBLEM FORMULATION A. System Model and Parameters
In this paper, we consider a backhaul network connecting a set B = {b 1 , · · · , b M } of M base stations. We assume that all nodes (interchangeably denoting base-stations in this paper) have line-of-sight connections. Note that, with minor modifications, the results claimed in this paper apply to a network with pre-deployed OF links (scenario of [11], [12], [21]). However, due to space limitation, we will only focus on networks with no existing links. Each node can be connected to any other node with either an OF or a hybrid RF/FSO connection, as in Figure 1 which shows a network containing 5 base-stations.
Let d(., .) be the distance operator between any two nodes. In other words, d(b, b ) is the distance between the nodes b and b , ∀(b, b ) ∈ B 2 . Let π (O) (x) be the cost of an OF link and π (h) (x) the cost of a hybrid RF/FSO link as a function of distance x. We have π (O) (0) = π (h) (0) = 0. Both functions π (O) (x) and π (h) (x) are increasing functions of the distance x. Moreover, since hybrid RF/FSO is a cost effective solution, we assume that the minimum distance between two base-stations. Figure 2 represents an example of such cost of OF and hybrid RF/FSO links with respect to the distance.

B. Problem Formulation
This paper considers the problems of minimizing the network deployment cost under the following constraints: 1) Connections between nodes can be either OF or hybrid RF/FSO. 2) Each node has a data rate that exceeds the target data rate. 3) Each node can communicate with any other node through single or multiple hop links. In other words, the graph is connected.
Let X ij , 1 ≤ i, j ≤ M be a binary variable indicating if base stations b i and b j are connected with an OF connection. Similarly, let Y ij , 1 ≤ i, j ≤ M indicates if they are connected with an hybrid RF/FSO link.
The objective function can be written as: Hence, the objective function can be written as Moreover, since X ij = X ji , Y ij = Y ji , and d(b i , b j ) = d(b j , b i ), then the objective function can be written as: Define C = [c ij ] as the adjacency matrix as follows: Since only one type of connections exists between the same BSs then X ij Y ij = 0, hence c ij is a binary variable (i.e. c ij ∈ {0, 1}).
From a graph theory perspective [23], [24], the graph connectivity constraint is expressed as a function the Laplacian matrix L defined as The connectivity condition of the matrix can be written using the algebraic formulation proposed in [23], [24] as λ 2 > 0.
The rate constraint can be written as for all BS i as follows: where R (O) and R (h) are the rates of the OF and the hybrid RF/FSO links, respectively, and R t the target rate. For simplicity, in this paper, we consider the normalized rates. Hence the constraint can be written as: Since X ij and Y ij are binary variable, then the constraint can be rewritten, for all BS i, using the normalized rates as follows: where R (O) (x) and R (h) (x) are the normalized data rates of an OF and hybrid RF/FSO links, respectively, as a function of the distance x. By normalized rates, we refer to the minimum between the actual rate, divided by the target rate, and unity. In other words, if the rate is R, the target rate is R t and the normalized rate R n , then the relationship linking the three quantities can be expressed as R n = min(R/R t , 1). Since the OF provides high data rates, without loss of generality, we assume that is a decreasing function of the distance. Figure 3 represents an example of such normalized data rates of OF and hybrid RF/FSO links with respect to the distance. The problem of minimizing the cost of the backhaul network planning can be formulated as: where the optimization is over both binary variables X ij and Y ij . The optimization problem (8) is equivalent to a weighted Steiner tree problem which is NP-hard [25], [26]. We refer to the solution of this problem as the optimal planning (OPT-PL). In rest of this paper, we propose efficient method to solve the problem, under the assumption that the hybrid RF/FSO connection between two nodes that are far away from each other is always more expensive than the cost of the OF connections between one of the nodes and its closest neighbour. The rationale for such assumption is that the solution for such problem becomes mathematically tractable as shown in the next subsection.
The notations used in this paper is summarized in Table I.

III. HEURISTIC SOLUTION
As highlighted above, the original optimization problem (8) is an NP-hard problem. The difficulty in solving the problem rate of hybrid RF/FSO link of x unit of distance lies particularly in the structure of constraint (13) and in simultaneously optimizing (8) over both binary variables X ij and Y ij . In this section, we present an efficient heuristic to solve the problem under the assumption that the hybrid RF/FSO connection between two nodes that are far away from each other is always more expensive than the OF connections between one of the nodes and its closest neighbor. The heuristic is based on first finding the solution to the problem when only OF links can be used. Afterwards, it solves an approximate of the backhaul network planning problem via relating problem (8) to solution reached by the planning problem when only OF links are allowed.

A. Optimal Planning Using Optical Fibre Only
For simplicity, in this section we consider that the cost of OF links is linear with the distance. This assumption is not only realistic but also comes without loss of generality. Indeed, the general case be obtained by simply replacing the distance operator d(., .) by the corresponding cost of the OF link π (O) (d(., .)) in the problem formulation and the algorithm. The following lemma introduces the reduced problem when only OF links are allowed. Lemma 1. The problem of backhaul design with minimum cost when only OF links are allowed is the following: Algorithm 1 Optimal planning using only OF links Require: B, d(., .).
Proof: The proof of this lemma can be found in Appendix A To solve the aforementioned problem, we propose to cluster the BSs according to the minimal distance. First, a cluster is created for all BS in the system. Afterwards, the two clusters at minimum distance of each other are merged into a single cluster, and connected. The distance between two clusters is defined as the minimum distance between all BS in each cluster. When two clusters are merged, the two BSs, in each cluster, at minimum distance are connected through an OF link. The process is repeated until only one cluster remains in the system. The steps of the algorithm are summarized in Algorithm 1. The following theorem characterizes the solution produced by Algorithm 1 with respect to the problem defined in Lemma 1: Theorem 1. The solution reached by Algorithm 1 is the optimal solution to the problem proposed in Lemma 1. We refer to this solution as the optimal OF only planning (OF-O-PL).
Proof: The proof of this theorem can be founded in Appendix B.

B. Problem Approximation
In this section, we approximate the backhaul network planning problem (8) under the assumption that an hybrid RF/FSO connection between two nodes that are not neighbours is more expensive than an OF links between one of the nodes and its closest neighbour. We first define b i * as be the closest node to base-station b i as follows: The set of neighbours N i of base-station b i is defined as the set of base-station that are closest to base-station b i , and that satisfy the connectivity condition. Mathematically, the condition can be written as: where X ij , 1 ≤ i = j ≤ M is the optimal solution to the OF only planning problem (15).
The results presented in this paper do not depend on the definition of the set of neighbours N i of node b i as long as N i ⊂ N i . Intuitively, as the set N i gets bigger and bigger, the approximation of the solution is more and more tight. For N i = B \ {b i }, the proposed algorithm reduces to an exhaustive search.
The assumption that two nodes that are far away from each others (i.e. not neighbors) connected with hybrid RF/FSO link generate a cost greater that the cost of the same nodes connected with OF links with their closest neighbors can be written as: Based on the above assumption, The following lemma approximates the optimization problem (8) under the assumption (21). Lemma 2. The problem of backhaul network cost minimization design using OF and hybrid RF/FSO connections can be approximated by the following problem: Proof: To proof this lemma, we first proof that any solution to (22), is a feasible solution to (8). To show that a feasible solution to (22), is a feasible solution to (8), we just show that constraint (27) is included in constraint (13) since it is the only constraint changing from one formulation to the other. Constraint (27) ensures that, for all connections X ij = 1 that are generated by Algorithm 1, then a similar connections (OF or hybrid RF/FSO link) between nodes b i and b j should exist. For connections X ij = 0, the constraint is always satisfied and connection may or may not exist. From Theorem 1, Algorithm 1 produces a connected graph. In other words, λ 2 > 0. Therefore, constraint (27) is included in constraint (13) which conclude that a feasible solution to (22), is a feasible solution to (8). In Theorem 2, we show that the optimal solution to (22) is the optimal solution to (8) in many scenarios (but not all). Therefore, the approximation of problem (8) by the problem (22) is tight.

C. Proposed Solution
This section proposes the solution for the approximate problem (22). The solution is based first on constructing the network planning graph, and then on formulating the problem (22) as a graph theory problem that can be optimally solved with moderate complexity.
1) Planning Graph: In this section, we introduce the undirected planning graph G(V, E), where V is the set of vertices and E the set of edges. Before stating the vertices construction and the edge connection, we first introduce the clustering C i for each node b i , 1 ≤ i ≤ M in the network as follows: Note that the size of a clustering C i is always strictly bounded by 2 |Ni|+1 . For each cluster α i ∈ C i , a vertex v ij , 1 ≤ j < 2 |Ni|+1 is generated. Two distinct vertices v ij and v kl are connected with an edge in E if the two following conditions are satisfied: 1) C1: i = k: The vertices represents different nodes in the network.
The vertices are non conflicting. 2) Proposed Algorithm: The following theorem characterize the solution of the approximated backhaul network planning problem (22).
ij , Y * ij ), 1 ≤ i, j ≤ M be the optimal solution to the planning problem (22) then we have The proof of this theorem can be found in Appendix C.
The following theorem links the solution of problem (22) to the planning graph.
Theorem 3. The solution of the approximation of the backhaul network design using hybrid RF/FSO can be formulated as a maximum weight clique, among the cliques of size M , in the planning graph, in which the weight of each vertex v ij is the weight of the corresponding cluster α i defined in (30).
Proof: The proof of this theorem can be found in Appendix D.

IV. SIMULATION RESULTS
This section shows the performance of the proposed solution to the backhaul network planning problem using hybrid RF/FSO technology. The base stations are randomly placed on a 5 Km long square. The cost of an OF link is taken to be 13.2$ per meter and the cost of an hybrid RF/FSO link is taken to be independent of the distance. Given the prices offered by the different constructors (fSONA, LightPointe, and RedLine), 2 type of price are identified π (F ) = 10000 and 20000. The price π (F ) = 40000 is proposed as a cut of price for which hybrid RF/FSO do not represent any advantage. The normalized data rate of an hybrid RF/FSO links is taken to be constant over a distance x after which it decays exponentially. For illustration purposes, the distance x is assumed to be 3Km unless indicated otherwise. The numbers of base-stations, price of the hybrid RF/FSO transceivers, and the distance x vary in the simulations so as to study the methods performance for various scenarios.
The optimal planning solution (solution of (8)) denoted by "OPT-PL", the planning solution using only OF links (solution of (15)) denoted by "OF-O-PL" and our proposed solution (solution of (22)) denoted by "FSO-PL" are simulated in this section. Figure 4 plots the cost of the network versus the number of BSs, for different costs of the hybrid RF/FSO transceivers. We clearly see that the degradation of our proposed solution against the optimal solution becomes less severe when first the number of base stations increases, and secondly when the hybrid RF/FSO transceivers becomes more expensive. The increase in performance in the first case can be explained by the fact that the connectivity opportunities of nodes increase as the number of base-stations increases, due to the increase in the neighbours sets N i . The gain in performance when the price of the hybrid RF/FSO transceivers increases can be explained by the fact that our assumption (21) becomes more valid as the price of the hybrid RF/FSO transceivers increases. Figure 5 and Figure 6 illustrate the cost of the network and the ratio of the OF link, respectively, against the cost of the hybrid RF/FSO transceivers. As shown in Figure 4, the performance of our proposed algorithm is more and more close to the one of the optimal planning as the cost of the hybrid RF/FSO transceivers increases. From Figure 6, we clearly see that if the hybrid RF/FSO transceivers are expensive enough, both the optimal and our proposed solution contain only OF links. In fact for expensive hybrid RF/FSO transceivers, the OF links offers a noticeable rate advantage which explain their use. It is worthy mentioning that for a cost π (h) ≥ 30000 in Figure 6, even though the link nature used in the optimal solution and our proposed solution are different, the total cost of the network is almost the same ( Figure 5 for M = 7).
Finally, to quantify the performance of the proposed algorithms with respect to the distance x, Figure 7 and  hybrid RF/FSO transceivers. Again, we clearly see that for expensive hybrid RF/FSO transceivers both the optimal solution and our proposed solution use exclusively OF links. Figure 7 shows that, even for cheap hybrid RF/FSO transceivers, our proposed solution performs as good as the optimal solution for a distance x ≥ 4. Figure 8 depicts that for a small x, our proposed solution uses more OF links than the optimal solution. Whereas for a x ≥ 4 the ratio is almost the same. This can be explained by our choice of neighbours N i . The connectivity opportunities of our proposed solution are less than the one of the optimal solution. Hence, for small x, to satisfy the rate constraint our proposed solution connects to the neighbours using OF links since the hybrid RF/FSO links do not satisfy the constraint. The optimal solution connects to more nodes (outside the neighbours sets) to satisfy the rate constraint.

V. CONCLUSION
In this paper, we consider the problem of backhaul network design using the OF and hybrid RF/FSO technologies. We first formulate the planning problem under connectivity and rate constraints. We, then, solve the problem optimally when only OF links are allowed. Using the solution of the OF deployment, we formulate an approximation of the general planning problem and show that under a realistic assumption about the relative cost of the OF links and the hybrid RF/FSO transceivers, the solution can be expressed as a maximum weight clique in the planning graph. Simulation results show that our approach shows a close-to-optimal performance, especially for practical prices of the hybrid RF/FSO.

APPENDIX A PROOF OF LEMMA 1
The problem formulation when only OF links are allowed can be simply obtained by setting Y ij = 0, 1 ≤ i, j ≤ M in (8). Constraints (10) and (11) become redundant and the problem can be written as: Note that since π (O) (x) is a linear function with respect to the distance, then the objective function can reduce to . Furthermore, since the graph is connected (i.e. λ 2 > 0), then each node is at least connected with another node. Hence, we have (12) become redundant. Therefore, the problem can be simplified to:

APPENDIX B PROOF OF THEOREM 1
To proof this theorem, we first proof that Algorithm 1 produces a feasible solution to the problem. Afterwards, we proof that any graph that can be reduced, using Algorithm 2, to a single cluster includes the graph outputted by Algorithm 1. Finally, we show that any solution that cannot be reduced to a single cluster is not optimal.
Algorithm 2 can be seen as a complement of Algorithm 1. As for Algorithm 1, in Algorithm 2 begins by generating a cluster of each BS in the system. Afterwards, two clusters at minimum distance of each other and whose BSs ta minimum distance are connected are merged in a single cluster. The process is repeated until no further connection can be found. The steps of the algorithm are summarized in Algorithm 2.
Let Z = {Z 1 , · · · , Z |Z| } be the clustering at any step of Algorithm 1. First note that Z is a partition of B. We proof by induction that nodes inside any cluster Z i , 1 ≤ i ≤ |Z| are connected. Clearly, for a cluster Z i with |Z i | = 1 (the cluster contains a single node), all nodes inside the cluster are connected. Assume that all clusters Z i of size |Z i | ≤ n are connected. From the last step of Algorithm 1 clusters of size n + 1 can be generated only by merging two clusters Z j and Z k with |Z j |, |Z k | < n and |Z j | + |Z k | = n + 1.

Algorithm 2 Clustering Algorithm
Since by construction such clusters are connected (X jk = X kj = 1 with b j ∈ Z j and b k ∈ Z k ), then the resulting cluster Z i from merging Z j and Z k is also connected. Therefore, all nodes within any arbitrary cluster Z i , 1 ≤ i ≤ |Z| are connected. Finally, since Z contains a single cluster at the end of Algorithm 1 (i.e. |Z| = 1) and it is a partition of B (contains all the nodes in the network), then all the nodes are connected. Hence the outputted solution satisfy constraint (17). By construction, we can easily see that the connections are binary and symmetric. In other words, the outputted solution satisfy constraints (16) and (18) which conclude that it is a feasible solution.
Let X ij , 1 ≤ i, j ≤ M be a feasible solution and let Z be the outputted clustering by Algorithm 2 when inputted X ij . We can clearly see that Z is a partition of B. Therefore, if |Z| = 1, then Z = Z since there exist only a unique partition containing a single element (the set B itself). This concludes that the any graph that can be reduced, using Algorithm 2, to a single cluster includes all the connections that are created in the graph outputted by Algorithm 1. Since π (O) is a strictly positive function and that the graph outputted by Algorithm 1 have the minimum number of connections among all the graphs that can be reduced to a single cluster using Algorithm 2, then the solution of Algorithm 1 is the best solution among the solutions that can be reduced to a single cluster using Algorithm 2. Now assume that |Z| = 1. We can clearly see that |Z| ≥ 3. Otherwise, if there exist only two clusters (i.e. |Z| = 2) and since the solution is feasible, then they are connected. Due to the fact they are only two, then two cases can be distinguished: • The BSs at minimum distance of each others are connected and hence they can be reduced to a single cluster. Therefore, clustering |Z| = 2 cannot be outputted by Algorithm 2. • The BSs at minimum distance of each others are not connected. Then the solution having the same connections except for the connection between the two cluster being replaced with the connection of BSs at minimum distance of each others produces a feasible solution at a lower cost. Therefore, the initial solution is not optimal.
Let Z = {Z 1 , · · · , Z |Z| }, |Z| ≥ 3 be the outputted clustering by Algorithm 2. Define the reduced clustering as We can clearly see that |Z| = 3 with none of the clusters connected and at minimum distance of each other.
Lemma 3. For any three points in the plane, there must exist two points at minimum distance of each others.
Proof: Let a, b, c be the three points in the plane and assume that there do not exist two points at minimum distance of each other. The only possible configuration (up to a permutation of the points) is that a at minimum distance from b and b is not. Hence b is at minimum distance from c which is at its turn at minimum distance from a. These conditions yield d(a, b) < d(a, c), d(b, c) < d(b, a) and d(c, a) < d(c, b). Since the distance operator is symmetric then, d(c, a) < d(c, b) = d(b, c) < d(b, a) = d(a, b) < d(a, c). In other words, d(a, c) < d(a, c), which is impossible. Therefore for any three points in the plane, there must exist two points at minimum distance of each others.
From Lemma 3, there must exist two clusters at minimum distance of each other. Since they have not been reduced to a single cluster, then they are not connected. For simplicity, assumeZ 1 andZ 2 are at minimum distance of each other and since the graph is connected thenZ 1 andZ 3 are connected and similarly forZ 2 andZ 3 . Moreover, it can be easily concluded that all nodes inside the clusters are connected. Otherwise, assumeZ 3 can be split into two non connected clustersZ and Z withZ 1 connected only toZ andZ 2 connected only toZ . Then, sinceZ 1 andZ 2 are not connected, the whole graph is not connected and the solution is not feasible.
The clustering connectingZ 1 withZ 2 andZ 1 withZ 3 (or Z 2 withZ 3 ) produces also a feasible solution at a lower cost since the sum of the distances is minimized (d(Z 1 ,Z 2 ) < min(d(Z 1 ,Z 3 ), d(Z 2 ,Z 3 ))). Therefore,Z is not optimal and by extension Z is also not optimal. Finally, we can conclude that the optimal solution is the solution containing a single cluster. Therefore, the solution outputted by Algorithm 1 is the optimal solution to the problem proposed in Lemma 1.

APPENDIX C PROOF OF THEOREM 2
In this theorem, we show that the optimal solution Remark 2. Note that if any feasible solution X ij , Y ij to the general problem (8) that verify Theorem 2, then the solution is feasible to (22). In other words, if X ij , Y ij feasible to (8) and This can be easily concluded given the construction of N i , ≤ i ≤ M as the minimum set of nodes that can generate a connected graph. Since all (x, y) such that X xy = 1 are at the edge of at least one of the N i , then the only connected solution that satisfies Theorem 2 is feasible to (22). In that scenario, the optimal solution of (8) and (22) are the same.
Assume ∃(x, y) such that X * xy + Y * xy = 1 and (x, y) / ∈ N y × N x . Two scenarios can be distinguished: For the first scenario, consider the reduced network {b x , b y ,b k }. Clearly, we have b x * , b y * ∈b k . Define the following planning: in which the connection between b x and b y is replaced by two connections between b x and b x * and between b y and b y * . We can clearly see that the network is connected. Moreover, since b x and b y are connected with an OF link, then the data rate constraint is satisfied for both nodes. Therefore, X ij ,Ỹ ij , 1 ≤ i, j ≤ M represents a feasible solution. Moreover, The difference in cost between the optimal planning X * ij , Y * ij and the planningX ij ,Ỹ ij is lower bounded by: From assumption (21), the difference is positive. This concludes that X * i,j , Y * i,j , 1 ≤ i, j ≤ M is not the optimal solution.
Remark 3. For scenario 1, X * i,j , Y * i,j , 1 ≤ i, j ≤ M can be the optimal solution to the original problem (8). Hence, for this configuration, the optimal solution of (8) and (22) are the same.
For scenario 2, let the network be reduced to {b x , b y ,b k ,b l } with b x connected tob k , which is a connected subgraph, b y connectedb l , which is a connected subgraph, andb k andb l are not connected. Since X * i,j , Y * i,j , 1 ≤ i, j ≤ M is a feasible solution to (22), then it satisfies constraint (27). In other words, X * k,l +Y * k,l = 1, ∀ k, l such that X kl = 1. Note that X xy = 0. Otherwise, by construction of the neighbours sets, we have (x, y) ∈ N y × N x . Define the planningX ij ,Ỹ ij such that The planningX ij ,Ỹ ij satisfy (27). However, we can clearly see that the graph is not connected. Therefore, X * i,j , Y * i,j , 1 ≤ i, j ≤ M is not a feasible solution. This concludes that scenario 2 is not feasible. Finally, we conclude that the optimal solution X * i,j , Y * i,j , 1 ≤ i, j ≤ M to (22) should satisfy X * xy + Y * xy = 1 only if (x, y) ∈ N y × N x . Remark 4. Scenario 2 can be a feasible scenario if X * i,j , Y * i,j , 1 ≤ i, j ≤ M is the optimal solution to the original problem (8). In that case, two scenarios can be distinguished: • b x * ∈ N y and b y * ∈ N x . In that caseX ij ,Ỹ ij presented for scenario 1 produces a feasible solution with lower cost. Therefore, the optimal solution of (8) and (22) are the same. • b x * / ∈ N y or b y * / ∈ N x . In this configuration, no conclusion can be reached about the optimal solution of (8) and (22) is an upper bound of the minimum of (8).
Therefore, the approximation of the problem (8) by the problem (22) is a tight approximation.

APPENDIX D PROOF OF THEOREM 3
To proof this theorem, we first prove that there is a one to one mapping between the set of feasible solution of a modified version of problem (22) and the set of cliques of degree M in the planning graph G(V, E). To conclude the proof, we show that the weight of the clique is equivalent to the merit function of the optimization problem (22).
From Theorem 2, we have X ij = 0 and Y ij = 0, ∀ (b i , b j ) / ∈ N j × N i . Hence the objective function (22) and constraint (26) can be replaced by: Similarly, since X ij = 0 and Y ij = 0, ∀ (b i , b j ) / ∈ N j ×N i , then (X ij , Y ij ) = (X ji , Y ji ) is always verified for (b i , b j ) / ∈ N j × N i . Hence the constraints (23), (24), and (25) can be replaced by: By definition of the set N i , we have X ij = 0, ∀ j / ∈ N i . Therefore, constraint (27) may be written as: 3) The problem (22) can be reformulated as: } be the new variable. Using the variables α i and the definition of the sets C i , the problem can be written as: (D.5b) Let C be the set of cliques of degree M in the planning graph and let F be the set of feasible solutions to the optimization problem (D.5). We first proof that any clique C = {v 1 , · · · , v M } ∈ C satisfy constraints (D.5a), and (D.5b). Then, we prove the converse. In other words, for element in F, there exists a clique in C.
Let C = {v 1 , · · · , v M } ∈ C. Assume ∃k, i, j such that v i , v j ∈ C k . Since all the vertices in a clique are connected, then from the connectivity condition C1, vertices v i and v j are not connected. Hence k, i, j such that v i , v j ∈ C k . Given that the clique contain M elements, then constraint (D.5a) is satisfied. The connectivity condition C2 ensures that (X ij , Y ij ) = (X ji , Y ji ) for all vertices. Therefore, C is a feasible solution to (D.5). Similarly, let {c 1 , · · · , c M } be a feasible solution to (D.5), then clearly the vertices corresponding to each cluster are connected. Finally, there a one to one mapping between C and F. Moreover, the weight of the clique is w(C) = M i=1 w(α i ) which conclude that the solution of (D.5) is the maximum weight clique, among the clique of size M , in the planning graph.


MIMO Precoding for Networked Control Systems with Energy Harvesting Sensors

In this paper, we consider a MIMO networked control system with an energy harvesting sensor, where an unstable MIMO dynamic system is connected to a controller via a MIMO fading channel. We focus on the energy harvesting and MIMO precoding design at the sensor so as to stabilize the unstable MIMO dynamic plant subject to the energy availability constraint at the sensor. Using the Lyapunov optimization approach, we propose a closed-form dynamic energy harvesting and dynamic MIMO precoding solution, which has an event-driven control structure. Furthermore, the MIMO precoding solution is shown to have an eigenvalue water-filling structure, where the water level depends on the state estimation covariance, energy queue and the channel state, and the sea bed level depends on the state estimation covariance. The proposed scheme is also compared with various baselines and we show that significant performance gains can be achieved.

I. INTRODUCTION
Networked control systems (NCSs) have become quite popular recently due to their growing applications in industrial automation, smart transportation, remote robotic control, etc.. A typical NCS consists of a multiple-input multiple-output (MIMO) dynamic plant 1 , a multiple antenna wireless sensor 2 and a controller. These are connected over a communication network and form a closed-loop control, as illustrated in Figure  1. Compared to traditional wireless sensors that are powered by a non-rechargeable battery, we consider an energy harvesting sensor, which is equipped with an energy harvesting device (e.g., a solar panel or a micro wind turbine) so that the sensor can harvest energy from the surrounding environment. Such NCSs with energy harvesting sensors have many advantages. For example, they do not require the replacement of batteries (for traditional sensors) and the maintenance process is simplified for NCSs in dangerous environments (e.g., chemical plants).
For an NCS, a potentially unstable plant is stabilized using a sensor to measure the state and feeds back to a controller, which in turn drives a stabilizing control signal to the plant. 1 MIMO dynamic plant refers to the dynamic plant with vector state evolution and vector control inputs [1]. 2 In some typical applications (such as 2.4GHz ISM WLAN and IEEE 802. 15.4), the number of available channels for NCS system operation can be very limited [1]. As such, multi-antenna sensor will be useful to enhance the spectral efficiency despite the limited bandwidth in the system. Unlike conventional feedback control design where the feedback path is assumed to be perfect, the sensor in the NCS observes the plant state and transmits the observation to a remote controller over a wireless MIMO fading channel. The performance of the NCS is closely related to the communication resource control at the energy harvesting sensor over the MIMO wireless channel. For an energy harvesting sensor, the communication resource is related to the energy available in the battery. Due to the random nature of the renewable energy source, it is difficult to predict the future evolution of the energy arrivals, and hence there is a tradeoff of using the available energy to support the current transmission and saving the energy for future good transmission opportunities.
In fact, NCS is a multidisciplinary subject involving both control theory and communication theory. The main objective of the control theory is to control the plant such that its state evolves in a desired manner. The plant is usually modeled as a linear dynamic system which is described by a set of first order coupled linear difference equations representing the evolution of the state variables. The representation of a plant with linear system provides a convenient and compact way to model and analyze the plant. Stability is an important characteristic of control system. When a plant is unstable, the state of the plant may be unbounded even though the input to the plant is bounded. General uncertainties and external disturbances may cause the dynamic plant to be unstable and this may incur costly physical damage. For instance, an unstable aircraft may crash and an unstable chemical plant may explode. Therefore, it is extremely important stabilize the unstable dynamic plant.
The theory of Lyapunov drift has a long history in the field of stochastic control to analyze the system stability in the field of control and communications. The authors of [2] first applied the Lyapunov drift theory to develop a general algorithm which stabilizes a multi-hop packet radio network. The negative Lyapunov drift terms play a central role when applying the Lyapunov drift theory to analyze the stability of dynamic systems. Intuitively, the negative Lyapunov drift is a stabilizing force that pulls the system state back to the equilibrium point. In [3], it is shown that negative Lyapunov drift ensures network stability because whenever the data queue length vector leaves a certain bounded region, the negative drift eventually drives it back to that region. In [4], the negative Lyapunov drift is utilized to analyze the stability of Markov chains.
From the communication side, it is important to under-stand how to optimize the communication resources (such as MIMO precoding) targeted for control applications. In particular, MIMO precoder optimization has been well studied in wireless communications. A MIMO precoder is essentially a multimode beamformer which splits the transmit signal into spatial eigenbeams and assigns higher power along the beams when the channel is strong [5]. In [6], [7], the authors obtain the MIMO linear precoding solution to minimize the MMSE [6] or maximize the SINR [7] of the MIMO wireless systems. However, these precoding solutions are not tailored for NCS applications because the optimization objectives (MMSE or SINR) are merely physical layer metrics in wireless communications, which may not be directly related to the performance metric of the NCS applications (such as stability of the plant or plant state estimation errors). Another related development in wireless communications is the device-to-device (D2D) or machine-to-machine (M2M) communications where a communication device can communicate and exchange information with a peer device autonomously. One important application of D2D or M2M systems (in addition to delivering content) is to support real-time industrial control where the devices may act as sensors and/or actuators. As such, these new application scenarios embrace both wireless communications and control in the context of networked control systems. In this paper, we are interested in studying the design of the MIMO precoder in the communication subsystem to support NCS applications.
Note that communication resource allocation for an NCS with an energy harvesting sensor is quite challenging. This is because such a problem embraces the information theory (to model the physical layer over the wireless channels), the queuing theory (to model the energy queue dynamics) and the control theory (to model the plant dynamics under imperfect state feedback control). There are some works on the dynamic resource control design for NCSs with an energy harvesting sensor. In [8], the authors study the mean square average state estimation error minimization under renewable energy constraints. To obtain the optimal communication power control policies, the associated stochastic optimization problems are solved using the numerical value iteration algorithm, which induces huge complexity, and suffers from slow convergence and lack of insights [9]. In this paper, we propose a lowcomplexity closed-form dynamic MIMO precoding solution at the sensor powered by renewable energy to stabilize the unstable MIMO dynamic plant. The following summarizes the key contributions in the paper.
• Direct Analog State Transmission with Peak Power Constraint: Unlike existing approaches in NCS [10], [11], [12] where the plant state is first quantized and transmitted over a simplified digital channel, we propose a novel analog state transmission where the sensor simply transmits a spatially rotated state measurement (rotated by the MIMO precoder) to the remote controller without quantization or coding. The consideration of renewable energy source also poses a unique challenge. For instance, the sensor needs to equip with a limiter, which is a nonlinear module to attenuate peaks of signals to satisfy the peak power constraint 3 . Such design will simplify the datapath design of the sensor and the remote controller and directly take advantage of the MIMO communication channels. • Closed-form Dynamic Event-Triggered MIMO Precoder Policy: Event-triggered sensor control has been proposed in various existing NCS applications [12], [13]. However, the existing solutions either has no closedform triggering solution [12], [13] or the solution is not truly dynamic [12]. In addition, in all these existing solutions, there is no consideration of the multiantenna MIMO communication channels. In this paper, we derive a closed-form fully dynamic solution which adapts to the complete system state (MIMO channel state for transmission opportunity, sensor energy state for energy availability and plant state estimation error for the urgency of the state transmission). These agile adaptivity are very important for superb performance in the NCS application. • Closed-form Performance Characterizations: In this paper, we also derive closed-form requirement of the renewable energy arrival rate and the battery capacity requirement to attain stability of the MIMO NCS. Such results give important guideline for the dimensioning of the resource needed at the sensor. Furthermore, we have derived closed form MSE and study how the performance depends on key system parameters. Such closed-form derivation is challenging due to the dynamic MIMO precoding policy as well as the coupled system state evolutions where the dynamic evolution of the three system states (the energy state, the plant state and the MIMO channel state) are tightly coupled together in a very complicated manner. Notations: Uppercase and lowercase boldface denote matrices and vectors, respectively. The operators (·) T , (·) † , (·) H , Tr (·), | · |, Re {·}, 1{·} are the transpose, element-wise conjugate, conjugate transpose, trace, cardinality, real part, and indicator function, respectively; ||A|| and ||a|| denote spectrum norm of matrix A and Euclidean norm of vector a, respectively; and diag(a) means diagonal matrix with diagonal elements being a; A ij denotes the element in the i-th row and j-th column of matrix A; µ i (A) denotes the i-th largest eigenvalue of matrix A. R m×n (C m×n ) represents the set of m × n dimensional real (complex) matrices.
II. SYSTEM MODEL In this section, we introduce the model of the MIMO NCS with an energy harvesting sensor, plant dynamics, MIMO channel model, energy queue model, as well as the information structures at the sensor and the controller.
A. MIMO Networked Control System with an Energy Harvesting Sensor Figure 1 shows a typical networked control system (NCS) with a MIMO plant (potentially unstable), a multi-antenna sensor (with N s transmit antennas) with energy harvesting capability and a multi-antenna remote controller (with N c receive antennas). The sensor and the controller are geographically separated and connected through a wireless MIMO fading channel. We consider a time-slotted system with slot duration τ . The sensor has perfect observation of the MIMO plant state x (n) at every time slot. Due to the consideration of renewable energy source at the sensor, the observed plant state x (n) is passed through an energy limiter before the transmission so as to satisfy an instantaneous peak power constraint determined by the instantaneous available energy. The output of the limiter is q (n). Unlike conventional digital approaches in NCS, we consider direct analog transmission of the measured system state at the sensor. Specifically, the analog state measurement is spatially rotated by a MIMO precoder F(n) and transmits to the remote controller over the multi-antenna wireless fading channel. We assume the perfect channel state information can be obtained at the sensor by causal signaling feedback from the controller as illustrated in Figure 1. The received signal at the controller is y(n), which is passed to the state estimator to obtain a state estimatex (n). Based onx (n), the remote controller further generates the control action u (n). The controller is physically attached to the actuator of the plant so that the control action u (n) generated at the controller can be applied to the plant directly. The actuator, which is colocated with the plant, then applies u (n) for plant actuation. The goal of the MIMO NCS is to stabilize the potentially unstable MIMO plant with limited wireless communication resources.

F(n)q(n)
Signaling link

x(n)
Energy harvesting sensor

B. Stochastic Dynamic MIMO Plant Model
We consider a discrete-time stochastic MIMO plant system with state dynamics: x (n + 1) = Ax (n) + Bu (n) + w (n), n ≥ 0, x (0) = x 0 , where x (n) ∈ R K×1 is the plant state process, u (n) ∈ R D×1 is the plant control action, A ∈ R K×K , B ∈ R K×D , and w (n) ∈ R K×1 is the plant noise. We assume the plant noise w (n) ∈ R K×1 is zero mean with covariance E w (i) (w (j)) T = δ i,j W for ∀i, j, where δ i,j = 0 if i = j and δ i,j = 1 otherwise. We assume that the plant noise covariance is finite 4 . The MIMO plant system (A, B) is assumed to be controllable with A containing possibly unstable eigenvalues 5 .

C. MIMO Wireless Channel Model
We model the wireless communication channel between the multi-antenna sensor and the controller as a wireless MIMO fading channel. Using multiple-antenna techniques, the N santenna sensor can deliver K parallel data streams to the N cantenna controller through spatial multiplexing. We assume K ≤ min {N s , N c }. At the n-th time slot, the received signal y (n) at the controller is given by where H (n) ∈ C Nc×Ns is the MIMO channel fading matrix, F (n) ∈ C Ns×K is the MIMO precoding matrix, q(n) is the output of the energy limiter 6 , and z (n) ∼ CN (0, I Ns ) is the additive complex Gaussian noise. We have the following assumption on H (n) :

D. Energy Harvesting and Energy Queue Model at the Sensor
We assume the sensor is solely powered by renewable energy sources (such as a solar panel). Due to the random nature of the energy source, a battery or ultra-capacitor is needed to store the harvested energy at the sensor. Let α (n) be the amount of harvestable energy at the sensor at time slot n. We have the following assumption on harvestable energy process α (n).
Assumption 2: (Renewable Energy Model) The harvestable energy α (n) is i.i.d. across different time slots according to some distribution.
The sensor has an energy storage (or battery) so let E (n) be the amount of energy left in the storage device at time n. We assume that the sensor is causal in the sense that new energy arrivals are observed after the control actions are performed at each time slot. Hence, the energy queue dynamics at the sensor is given by where [x] + = max {x, 0}, F (n) q (n) 2 τ measures the energy consumed at time slot n for delivering q (n) using precoding action F (n) over the MIMO wireless fading channel, and θ is the sensor battery capacity. Furthermore, at any time slot n, the precoding control action must satisfy the following energy availability constraint: Challenge 1: Energy Availability Constraint and Saturation.
The energy availability constraint (3) means that the energy consumption at each time slot cannot exceed the current available energy in the energy buffer at the sensor. Such a constraint greatly complicates the design of the dynamic MIMO precoder at the sensor. At any time slot, the dynamic MIMO precoder has to strike a balance between how much energy to consume, the good transmission opportunities induced by the fading channel and the urgency induced by the state estimation errors at the controller. Furthermore, there is an effective peak transmission power constraint at the sensor at each time slot. In order to satisfy this constraint, the sensor needs to be equipped with a limiter. Since the input process x (n) is non-stationary with unbounded support, there is nonzero chance that the input x (n) exceeds the limiter range, causing saturation. This causes non-linearity in the feedback loop and substantially complicates the optimization problem. In the next section, we introduce the design of the limiter, and we will elaborate how we tackle the issue of saturation in Section III-B.

E. Energy Limiter at the Sensor
An energy limiter is a non-linear functional module, and the amplitude relationship between its input and output is illustrated in Figure 2. Specifically, an energy limiter is parameterized by its dynamic range L (n) and saturation value M . If the limiter input signal's amplitude does not exceed the dynamic range L(n), then the input signals are allowed to pass with a scaling of M L(n) . On the hand, if the input signal's amplitude is greater than the dynamic range L(n), then its amplitude is attenuated to the value M . The input of the energy limiter is the plant state x(n) and the output is q(n), and the specific structure of the energy limiter is given by: where The design of the limiter's dynamic range L (n) is very critical to the system performance. It's clear that, due to the nonlinearity of the energy limiter, the magnitude information of x (n) will be partially lost when the limiter saturates, i.e., x (n) > L (n), which will greatly deteriorate the state estimation quality at the remote controller. We shall address how to design the limiter dynamic range L (n) in Section III.

Magnitude of Plant State Input
The sensor then transmits the magnitude-limited plant state q (n) to the remote controller over a wireless MIMO fading channel, and based on the limiter structure in (4), we have q T (n) q (n) ≤ M 2 . Therefore, we restrict the energy availability constraint in (3) as follows 7 :

III. DYNAMIC MIMO PRECODING VIA LYAPUNOV OPTIMIZATION
In this section, we shall address the design of the energy limiter, and formulate the dynamic MIMO precoding problem to stabilize the unstable MIMO dynamic plant using Lyapunov optimization. We first have the following definition on the stability of the dynamic MIMO plant [15]: where the expectation is taken with respect to the randomness of the plant noise, the channel noise, the MIMO channel state, and the energy state.

A. Virtual State Estimation Covariance Process
The nonlinear structure of the limiter introduces the nonlinear state estimation, and thus the evolution of the plant state estimation mean square error (MSE) is very difficult to be characterized explicitly. To get around this, we first establish a tight analytical bound on the state estimation MSE. Letx (n) be the minimum mean square error (MMSE) MIMO plant state estimate at the controller. We have the following lemma, which gives an upper bound of the MSE of the plant state estimation where Σ (n) is a virtual state estimation covariance process with the following dynamics: 7 It means that if the constraint in (5) is satisfied, then the constraint in (3) is satisfied. For analytical tractability, we restrict the energy availability constraint (3) to constraint (5) so that the dynamics of state estimation error covariance is independent of the state x (n) and is therefore analytically tractable. Since our proposed MIMO precoding F(n) has an event-driven structure (as shown in Theorem 1) and the sensor can be in dormant mode most of the time to save energy, the difference in terms of average energy consumption between constraint (5) and constraint (3) is quite small. and Proof: Please see Appendix A. The upper bound of the MMSE in (6) can be achieved by using a low complexity state estimation 8 We consider a control law of the form u(n) = −ΨA x (n) 9 for the dynamic MIMO plant, where the feedback gain ΨA is chosen such that A − BΨA is a Hurwitz matrix. Note that the certainty equivalent controller is given by and P ∈ S K + and R ∈ S D + are the weighting matrices for the plant state deviation cost and plant control cost of the LQG control associated with the certainty equivalent controller [17], [18]. Hence, one possible way to design such Ψ is let From Lemma 1, if we can achieve stability in the virtual state estimation error covariance process, i.e., limsup 2 ) will also be bounded, which in turn leads to the bounded state process as in Definition 1. This is formally stated in the following Lemma. Lemma 2: (Connection between Stability of Σ (n) and Stability of x (n)) Under the control law Proof: Please see Appendix B. Note that the evolution of the virtual state estimation MSE process Σ (n) depends on the dynamic MIMO precoder F (n) according to (7). As a result, we will focus on the design of dynamic MIMO precoding F (n) in the MIMO NCS to achieve stability of Σ (n).

B. ε-Saturation Energy Limiter Design
The magnitude information of the plant state is partially lost when the limiter saturates. In conventional literature, bounded noise support is assumed, and hence saturation can be avoided with proper dynamic range design. In this work, we considered unbounded noise support, and hence saturation cannot be completely eliminated. Instead, we need to adapt the dynamic range L (n) of the limiter to maintain a small probability of 8 Note that due to the non-linearity of the energy limiter, it is difficult to compute the MMSE state estimatorx (n) . 9 Note that this form of control law is quite general and also covers the certainty equivalent control which has been widely used in adaptive control systems [16]. saturation ε. We have the following lemma on the design of the dynamic range of the limiter.
Lemma 3: (ε-Saturation Limiter Dynamic Range Adaptation) Suppose the dynamic range L (n) of the limiter evolves according to the following: where is a constant (Q and T are any positive definite symmetric ma- Then holds for any time slot n. Proof: Please see Appendix C. Remark 1: (Computation of L (n) at the Sensor) Since the input state process x (n) is a non-stationary random process, in order to maintain a small saturation probability ε, the dynamic range L (n) of the limiter needs to be adaptive as in (8) to keep track of the instantaneous covariance of the input state process x (n). Note that the limiter's dynamic range L (n) is a function of the virtual state estimation covariance Σ (n). The sensor can obtain Σ (n) via the dynamics (7) based on the local information only without explicit signaling feedback, and hence it is implementation friendly. Due to the local availability of Σ (n) and L (n) at the sensor, the MIMO precoding solution in Theorem 1 can be implemented at the sensor in a decentralized way.

C. Dynamic MIMO Linear Precoding via Lyapunov Optimization
We focus on deriving a dynamic MIMO precoding policy to achieve stability of the virtual covariance process Σ (n) using Lyapunov techniques [19], [20]. Intuitively, a negative term in the Lyapunov drift is a stabilizing force that pulls the system state back to the equilibrium point. As a result, we shall design the dynamic MIMO precoder to maximize the negative Lyapunov drift to achieve stability under the renewable energy resource constraints.
We define a Lyapunov function as follows [20]: The associated Lyapunov drift is given by where the expectation E is w.r.t. the randomness of the channel state and the energy state under a given MIMO precoding rule F(n). We have the following lemma on the Lyapunov drift.
Lemma 4: (Lyapunov Drift) Given a MIMO precoding rule F(n), the Lyapunov drift can be upper bounded as follows: Proof: Please see Appendix D. The first expectation on the R.H.S of (13) is taken w.r.t. the randomness of energy arrival, the second expectation on the R.H.S of (13) is taken w.r.t. the randomness of both the channel state and the energy state under given Σ (n) and MIMO precoding rule F(n). Note that the negative Lyapunov drift plays a central role when applying the Lyapunov drift theory to analyze the stability of dynamic systems. Intuitively, the negative Lyapunov drift is a stabilizing force that pulls the system state back to the equilibrium point. The drift term in equation (13) that contributes to positive Lyapunov drift is bracketed as a "bad" term for stabilization of Σ (n), because positive Lyapunov drift leads to instability of Σ (n). On the contrary, the drift term that contributes to negative Lyapunov drift is bracketed as a "good" term for stabilization of Σ (n), because negative Lyapunov drift is the stabilizing force for the stability of Σ (n). Hence, we focus on deriving a dynamic MIMO precoder to minimize the drift (13). This is equivalent to considering the following optimization problem: Problem 1: (MIMO Precoding via Lyapunov Optimization) For given realizations E(n) = E, Σ (n) = Σ, and H(n) = H, the MIMO precoding F (n) is given by the solution of the following problem: where L(Σ) is given in (8).
While Problem 1 is equivalent to a convex problem, it is still very challenging to obtain closed-form solution. This is due to the tight coupling in the objective function (involving a mixture of real, trace and inverse).

IV. EVENT-DRIVEN ENERGY HARVESTING AND MIMO PRECODING SOLUTION
In this section, we shall give the MIMO precoding solution to Problem 1 and discuss its structural properties.

A. Drift Minimizing MIMO Precoder
The underlying structure of the objective function of Problem 1 provides some key insights into the MIMO precoding solution structure. Specifically, if Re F H (n)H H HF(n) , Σ −1 and F H (n)F(n) can be simultaneously diagnalized, then the objective function is substantially simplified because the matrix operations will only involve diagonal matrices. This may be feasible if F(n) can be expressed as a product of three matrices, where the leftmost matrix is the Hermitian of the right singular matrix of H and the rightmost matrix is the Hermitian of the left right singular matrix of Σ. Inspired by this observation, let the singular value decomposition of the matrix of H be H = VΠU H , where V ∈ C Nc×Nc and U ∈ C Ns×Ns are unitary matrices, and Π is a rectangular diagonal matrix with diagonal elements in a descending order. Denote Π K as the leading principal minor of order K of Π. Let the eigenvalue decomposition of Σ −1 be Σ = SΛS T , where S ∈ C K×K is an orthogonal matrix, and Λ ∈ R K×K is diagonal with diagonal elements in a descending order. Then the drift-minimizing solution of Problem 1 is given as follows.
Theorem 1: (Drift-Minimizing Solution) The drift-minimizing MIMO precoding is as follows: +EI is not positive definite, then where β is given by (14). Proof: Please see Appendix D.

Remark 2: (Interpretation on the Structure of the Driftminimizing MIMO Precoding Solution)
Event-driven Structure: The MIMO precoding solution in Theorem 1 also has an event-driven structure, in the sense that the sensor either transmits or shuts down aperiodically depending only on whether the dynamic threshold θI − + EI is positive definite or not.
• Large Π K (good channel condition) or large E (sufficient energy in the storage) leads to active mode of the sensor, which means it is better that the sensor be active when there are good transmission opportunities. • Large MIMO plant state estimation error implies large Σ , as shown in Lemma 1, which also leads to the active mode of the sensor. This is because a large MIMO plant state estimation error means high transmission urgency and hence it is better that the sensor be active. Figure 3 illustrates a sample path of the state estimation error and transition between the active and dormant modes. It can be observed that the state estimation error x (n) −x (n) 2 increases during the dormant modes and is reset to a low value during the active modes. As such, the MIMO precoding solution has an event-driven structure with aperiodic reset of state estimation error.  Note that the virtual state estimation covariance Σ can be decomposed as Σ = SΛS H with S and Λ being unitary and diagonal, respectively. Hence, each diagonal element of Λ corresponds to a subsystem with Λ ii as the state estimation error of the i-th subsystem.
Dynamic Spatial Channel Activation Structure: The MIMO precoding solution F * (n) dynamically activates the spatial channels. Specifically, whether the i-th spatial channel (the spatial channel corresponds to the i-th subsystem) is activated or not depends only on whether the dynamic threshold + − β is positive or not. Eigenvalue Water-filling Structure: The MIMO precoding solution F * (n) also has an eigenvalue water-filling structure.
Specifically, for the i-th subsystem, AA T ([θ−E] + +β)τ is the dynamic water level which adapts to the virtual state estimation covariance Σ, energy storage E, and the i-th spatial channel (Π K ) ii , and (Λ) −1 ii is the dynamic seabed level. Large energy storage (i.e., large E) or better channel condition of the i-th spatial channel (i.e., large (Π K ) ii ) will lead to a high water level which will allow the sensor to allocate more transmission energy for the i-th subsystem. High transmission urgency (i.e., large (Λ) ii ) will lead to a low seabed level, which induces more energy to be allocated to the subsystem.

Remark 3: (Comparison with Traditional MIMO Precoding Solutions)
The traditional waterfilling MIMO precoding solution [21] that maximizing the physical layer capacity has the form for all i = 1, . . . , K, λ j are the Lagrange multiplier that enforce energy causality and µ j are the Lagrange multipliers that enforce no-energy overflow conditions. Similarly, the waterfilling solution under energy harvesting constraints also only adapts to the energy state and channel state. In contrast, the proposed MIMO precoding solution is truly dynamic and adaptive to the MIMO channel state, sensor energy state and plant state transmission urgency. It is noticed that the traditional waterfilling solution and the waterfilling solution with energy harvesting constraints both fail to exploit plant state transmission urgency and therefore, they have inferior performance as shown in Section VI.

B. Application Example
We consider the following simple toy example to illustrate the structural properties of the MIMO precoding solution.
Example 1: (Noisy Decoupled Plant System) We consider an aggregation of two decoupled SISO plant subsystems as illustrated in Figure 4. The dynamics of the two subsystems are given by: x 1 (n + 1) = 1.6x 1 (n) + u 1 (n) + w 1 (n) (Subsystem 1) where w 1 (n) and w 2 (n) are independent Gaussian random noise with zero mean and covariance 1. Let the control rule be u 1 (n) = −0.8x 1 (n) and u 2 (n) = −0.55x 2 (n). Let θ = 36, τ = 1 and the parameters of the ε-saturation limiter be ε = 0.1, M = 1 and Q = I. The energy-harvesting sensor observes the limiter output q (n) = (q 1 (n) , q 2 (n)), applies a precoder F = • When the Sensor has to be Active? Based on the event-driven structure in Remark 3, the sensor will be activated if either 2 is negative. Therefore, better channel condition (large h 1 or h 2 ), large energy storage E or increased transmission urgency (large σ 1 or σ 2 ) will lead to the active mode of the sensor.
• Which Spatial Channel to Turn On? Suppose h 1 = 4 and σ 1 = 70, Figure 5 illustrates the number of spatial channels activated at different regions of the state space (h 2 , σ 2 ) for different E. It can be observed that the region where both spatial channels are turned on enlarges as the available energy E increases. This is reasonable because large E means a good transmission opportunity which allows more spatial channels to be activated.  V. STABILITY ANALYSIS In this section, we derive the sufficient condition for NCS stability under the drift-minimizing MIMO linear precoding solution in Theorem 1, and analyze the MIMO plant state estimation performance.

A. Sufficient Condition for NCS Stability
Instability measures of the unstable dynamic plants provide key information directly related to stabilizability and performance limitations. One of the key measures that has been proposed in the literature for measuring the instability of the unstable dynamic plants is defined as the product of the modulus of the unstable eigenvalues of matrix A [23], [24], i.e., In this section, we shall establish a sufficient condition for NCS stability under the proposed MIMO precoding policy in terms of the instability measure M (A), the energy arrival α (n) as well as the battery capacity θ. The sufficient condition for NCS stability is summarized in Theorem 2. Theorem 2: (Sufficient Condition for NCS Stability) If the following condition is satisfied: and π is the unordered singular value of H, δ = 2 ε (1 + A − BΨA Θ) BΨ A and Θ is a constant given in Lemma 3, then the dynamic MIMO plant is stable.
Proof: Please see Appendix E. The sufficient condition in Theorem 2 delivers some key design insights into the dynamic MIMO plant system. Specifically, we have the following system design insights. . This means the more unstable the MIMO plant is, the larger the capacity of the energy storage devices is required. • Energy arrival requirement: The sufficient condition also implies that in order to achieve stability, there is a requirement on the minimum average energy arrival rate, i.e., , and larger average energy arrival rate (stronger energy harvesting capability) is required for the more unstable MIMO plant. The sufficient condition (16) is not difficult to check because it is expressed in terms of key system parameters from individual components of the system. For example, In (16), δ, K, M (A) and M AA T are the parameters from the MIMO dynamic plant. E 1 α is the statistical parameters for renewable energy source, which can be obtained from offline measurements. θ is the battery size of the sensor. ε is the sensor limiter saturation probability. Pr (π < ξ) and E π −1 π ≥ ξ are the statistical parameters for the MIMO fading channel, which again can be obtained from offline measurements.
Remark 4: In [25], [26], the authors show that to achieve stability there is a minimum data rate requirement in terms of the instability measure (15). In our scenario, if there is no sufficient energy storage, the sensor tends to be in sleep mode and there is no data transmission between the sensor and the controller. Hence, the available transmission energy at the sensor will affect the communication data rate. This intuitively indicates that to achieve stability there shall exist some requirement on the available transmission energy at the sensor. This effect is shown in the sufficient condition (16), where there is a requirement on the minimum average energy arrival rate in terms of instability measure (15).

B. State Estimation MSE Performance
We are interested in analyzing the achievable state estimation MSE using the proposed MIMO linear precoding policy. This is summarized in the theorem below.
Theorem 3: (MSE of State Estimation Error) If the sufficient condition for NCS stability (16) is satisfied, then the MSE satisfies: where Proof: Please see Appendix G. The state estimation performance bound (17) reveals the fact that large average energy arrival rate (small E 1 α ) and good MIMO channel quality (small E π −1 π ≥ ξ * and Pr (π < ξ * )) will result in a better state estimation performance.

VI. NUMERICAL RESULTS
In this section we compare the performance of the proposed MIMO linear precoding scheme with the following baselines via numerical simulations.
• Baseline 1 (MIMO Water-filling Precoding Maximizing Capacity) [27], [21]: The MIMO precoding action is given by The simulation is Monte Carlo simulation and for a given pair of average energy arrival E [α] and battery capacity θ, we simulate 5000 sample paths of plant state evolution and estimation, and each sample path contains 300 time slots.
The system parameter configurations, namely the harvestable energy process α (n), battery capacity θ, MIMO channel H, limiter saturation probability ε, and limiter parameter δ of Baseline 1-3 are the same as the proposed scheme, and the system parameter configurations satisfy the stability condition (16) in Theorem 2. The Baseline 1-3 and the proposed scheme also satisfy the peak power restriction and the energy harvesting constraints.
A. State Estimation MSE Versus Battery Capacity θ Figure 6 illustrates the normalized MSE of the state estimation error versus the sensor battery capacity θ under the average energy arrival E [α] = 40J. Note that Baseline 1 is the MIMO precoding that adapts to the CSI and energy state only and maximizes the physical layer capacity. Baseline 2 covers a lot of existing NCS literature where a static transmission scheme with the sensor periodically being active is adopted [8]. Baseline 3 is the MIMO precoding that maximizes the receiver MSE. Baseline 4 and 5 are supplied with constant power E [α] to compare the effect of energy harvesting constraints.
It can be observed that there is a significant performance gain of the proposed scheme compared with all the baselines. This is because the proposed scheme is plant state and CSI aware and is tailored for NCS applications. Baseline 1 is not tailored for NCS applications in the sense that Baseline 1 merely optimizes the physical layer throughput in wireless communications. Baseline 2 has the worst performance because it neither fully exploits the CSI nor adapts to the plant state. Baseline 3 has better MSE performance compared with the capacity maximizing solution (Baseline 1  VII. CONCLUSION In this paper, we consider an NCS with an energy harvesting sensor that delivers the MIMO dynamic plant state information to a controller over a MIMO fading channel. We derived a low-complexity closed-form dynamic MIMO precoder at the sensor via Lyapunov optimization to stabilize the unstable MIMO dynamic plant. The solution has an event-driven structure and an eigenvalue water-filling structure. We also analyzed the sufficient condition for NCS stability under the proposed MIMO precoder solution. Compared with the various existing MIMO linear precoding schemes, the proposed MIMO precoder solution has substantial state estimation performance gain.

B. Proof of Lemma 2
Note that x (n + 1) = Ax (n) − BΨA x (n) + w (n), it follows: Since A − BΨA is a Hurwitz matrix, then Tr (W) . Combining with Lemma 1 and note that Similarly, if A − BΨA = 0, then Therefore, Lemma 2 is proved.

C. Proof of Lemma 3
By the standard Lyapunov stability theory for discrete-time linear systems, there exist positive definite symmetric matrices Q and T such that It follows [31]: , it follows: although the one step drift x T (n + 1) Qx (n + 1) − x T (n) Qx (n) is positive, we have • If x (n) > Θ BΨA x (n) − x (n) + w (n) , the one step drift x T (n + 1) Qx (n + 1)−x T (n) Qx (n) is negative, hence x (n + 1) will stay in the ball: x (n + 1) : Hence holds for any time slot n.

Note
that ≤ ε holds, then taking expectation on both sides w.r.t. the randomness of I c (n − 1) will yield E 1 x (n) 2 > L 2 (n) = Pr x (n) 2 > L 2 (n) ≤ ε. Therefore, in order to guarantee the saturation probability of the limiter is at most ε, it is sufficient to design the limiter range L (n) such that E 1 x (n) 2 > L 2 (n) I c (n − 1) ≤ ε holds for all n.
. Let X(n) = M L(Σ) S∆ 1 G (n). Then Problem 1 is transformed into the following problem: AA T X(n) 2X T (n)X(n) + Σ −1 −2 Note that for non-convex problems, the KKT condition is necessary and in order to find the global optimum, we need to test all the solutions that satisfy the KKT conditions. Therefore, without loss of optimality, we restrict the feasible solution set of the optimization problem to be the solutions that satisfy the KKT conditions (37)-(40). By exploiting the structural properties of the solutions that satisfy the of the KKT conditions, the objective function of the optimization problem can be further simplified. Let X(n) = S T X(n)S, the KKT condition (37)-(40) are equivalent to the following conditions (41)-(44).

F. Proof of Theorem 2
Intuitively, the negative Lyapunov drift is related to the stability of the MIMO plant. Therefore, the sufficient condition for stability should be the sufficient condition to guarantee negative Lyapunov drift. Substituting the MIMO precoding solution in Theorem 1 into (36), it follows: whereπ = π Tr(Π −1 Taking limits of the inequality (57) as N → ∞ yields Note that from Lemma 1, we have limsup n=0 E[Tr (Σ (n))], hence Theorem 3 is proved.