{'fraction_non_alphanumeric': 0.03739284102872612, 'fraction_numerical': 0.011486689727778613, 'mean_word_length': 4.765553869499241, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cooperative Localization in Visible Light Networks: Theoretical Limits and Distributed Algorithms

Light emitting diode (LED) based visible light positioning (VLP) networks can provide accurate location information in indoor environments. In this manuscript, we propose to employ cooperative localization for visible light networks by designing a VLP system configuration that involves multiple LED transmitters with known locations (e.g., on the ceiling) and visible light communication (VLC) units equipped with both LEDs and photodetectors (PDs) for the purpose of cooperation. First, we derive the Cram\'er-Rao lower bound (CRLB) and the maximum likelihood estimator (MLE) for the localization of VLC units in the proposed cooperative scenario. To tackle the nonconvex structure of the MLE, we adopt a set-theoretic approach by formulating the problem of cooperative localization as a quasiconvex feasibility problem, where the aim is to find a point inside the intersection of convex constraint sets constructed as the sublevel sets of quasiconvex functions resulting from the Lambertian formula. Next, we devise two feasibility-seeking algorithms based on iterative gradient projections to solve the feasibility problem. Both algorithms are amenable to distributed implementation, thereby avoiding high-complexity centralized approaches. Capitalizing on the concept of quasi-Fej\'er convergent sequences, we carry out a formal convergence analysis to prove that the proposed algorithms converge to a solution of the feasibility problem in the consistent case. Numerical examples illustrate the improvements in localization performance achieved via cooperation among VLC units and evidence the convergence of the proposed algorithms to true VLC unit locations in both the consistent and inconsistent cases.

Abstract-Light emitting diode (LED) based visible light positioning (VLP) networks can provide accurate location information in indoor environments. In this manuscript, we propose to employ cooperative localization for visible light networks by designing a VLP system configuration that involves multiple LED transmitters with known locations (e.g., on the ceiling) and visible light communication (VLC) units equipped with both LEDs and photodetectors (PDs) for the purpose of cooperation. First, we derive the Cramér-Rao lower bound (CRLB) and the maximum likelihood estimator (MLE) for the localization of VLC units in the proposed cooperative scenario. To tackle the nonconvex structure of the MLE, we adopt a set-theoretic approach by formulating the problem of cooperative localization as a quasiconvex feasibility problem, where the aim is to find a point inside the intersection of convex constraint sets constructed as the sublevel sets of quasiconvex functions resulting from the Lambertian formula. Next, we devise two feasibilityseeking algorithms based on iterative gradient projections to solve the feasibility problem. Both algorithms are amenable to distributed implementation, thereby avoiding high-complexity centralized approaches. Capitalizing on the concept of quasi-Fejér convergent sequences, we carry out a formal convergence analysis to prove that the proposed algorithms converge to a solution of the feasibility problem in the consistent case. Numerical examples illustrate the improvements in localization performance achieved via cooperation among VLC units and evidence the convergence of the proposed algorithms to true VLC unit locations in both the consistent and inconsistent cases.

I. INTRODUCTION A. Background and Motivation
Accurate wireless positioning plays a decisive role in various location-aware applications, including patient monitoring, inventory tracking, robotic control, and intelligent transport systems [1]- [4]. In the last two decades, radio frequency (RF) based techniques have commonly been employed for wireless indoor positioning [5], [6]. Recently, light emitting diode (LED) based visible light positioning (VLP) networks have emerged as an appealing alternative to RF-based solutions, providing high-accuracy and low-cost indoor localization services [7]. While visible light networks can be harnessed for enhancing localization performance in indoor scenarios [8], they also offer illumination and high speed data communications simultaneously without incurring additional installation costs via the use of existing LED infrastructure [9]. In indoor VLP networks, various position-dependent parameters such as received signal strength (RSS) [10], [11], time-of-arrival (TOA) [12], [13], timedifference-of-arrival (TDOA) [14] and angle-of-arrival (AOA) [15] can be employed for position estimation. In order to quantify performance bounds for such systems, several accuracy metrics are considered, including the Cramér-Rao lower bound (CRLB) [10], [12], [13], [16] and the Ziv-Zakai Bound (ZZB) [17].
Part of this work was presented at IEEE International Black Sea Conference on Communications and Networking (BlackSeaCom), June 2017, www.ee.bilkent.edu.tr/∼gezici/coopVLP.pdf noncooperative approach, position estimation is performed by utilizing only the measurements between anchor nodes (which have known locations) and agent nodes (the locations of which are to be estimated) [18], [19]. On the other hand, cooperative systems also incorporate the measurements among agent nodes into the localization process to achieve improved performance [19]. Benefits of cooperation among agent nodes are more pronounced specifically for sparse networks where agents cannot obtain measurements from a sufficient number of anchors for reliable positioning [20]. There exists an extensive body of research regarding the investigation of cooperation techniques and the development of efficient algorithms for cooperative localization in RF-based networks (see [18]- [20] and references therein). In terms of implementation of algorithms, centralized approaches attempt to solve the localization problem via the optimization of a global cost function at a central unit to which all measurements are delivered. Among various centralized methods, maximum likelihood (ML) and nonlinear least squares (NLS) estimators are the most widely used ones, both leading to nonconvex and difficult-to-solve optimization problems, which are usually approximated through convex relaxation approaches such as semidefinite programming (SDP) [21]- [23], secondorder cone programming (SOCP) [24], [25], and convex underestimators [26]. In distributed algorithms, computations related to position estimation are executed locally at individual nodes, thereby reinforcing scalability and robustness to data congestion [19]. Set-theoretic estimation [27]- [30], factor graphs [19], and multidimensional scaling (MDS) [31] constitute common tools employed for cooperative distributed localization in the literature.
Despite the ubiquitous use of cooperation techniques in RFbased wireless localization networks, no studies in the literature have considered the use of cooperation in VLP networks. In this study, we extend the cooperative paradigm to visible light domain. More specifically, we set forth a cooperative localization framework for VLP networks whereby LED transmitters on the ceiling function as anchors with known locations and visible light communication (VLC) units with unknown locations are equipped with LEDs and photodetectors (PDs) for the purpose of communications with both LEDs on the ceiling and other VLC units. The proposed network facilitates the definition of arbitrary connectivity sets between the LEDs on the ceiling and the VLC units, and also among the VLC units, which can provide significant performance enhancements over the traditional noncooperative approach employed in the VLP literature. Based on the noncooperative (i.e., between LEDs on the ceiling and VLC units) and cooperative (i.e., among VLC units) RSS measurements, we first derive the CRLB and the MLE for the localization of VLC units. Since the MLE poses a challenging nonconvex optimization problem, we follow a set-theoretic estimation approach and formulate the problem of cooperative localization as a quasiconvex feasibility problem (QFP) [32], where feasible constraint sets correspond to sublevel sets of certain type of quasiconvex functions. The quasiconvexity arising in the problem formulation stems from the Lambertian formula, which characterizes the attenuation level of visible light channels. Next, we design two feasibility-seeking algorithms, having cyclic and simultaneous characteristics, which employ iterative gradient projections onto the specified constraint sets.
From the viewpoint of implementation, the proposed algorithms can be implemented in a distributed architecture that relies on computations at individual VLC units and a broadcasting mechanism to update position estimates. Moreover, we provide a formal convergence proof for the projection-based algorithms based on quasi-Féjer convergence, which enjoys decent properties to support theoretical analysis [33].

B. Literature Survey on Set-Theoretic Estimation
The applications of convex feasibility problems (CFPs) encompass a wide variety of disciplines, such as wireless localization [27]- [30], [34], compressed sensing [35], image recovery [36], image denoising [37] and intensity-modulated radiation therapy [38]. In contrary to optimization problems where the aim is to minimize the objective function while satisfying the constraints, feasibility problems seek to find a point that satisfies the constraints in the absence of an objective function [35]. Hence, the goal of a CFP is to identify a point inside the intersection of a collection of closed convex sets in a Euclidean (or, in general, Hilbert) space. In feasibility problems, a commonly pursued approach is to perform projections onto the individual constraint sets in a sequential manner, rather than projecting onto their intersection due to analytical intractability [39]. The work in [27] formulates the problem of acoustic source localization as a CFP and employs the well-known projections onto convex sets (POCS) technique for convergence to true source locations. Following a similar methodology, the noncooperative wireless positioning problem with noisy range measurements is modeled as a CFP in [28], where POCS and outer-approximation (OA) methods are utilized to derive distributed algorithms that perform well under non-line-of-sight (NLOS) conditions. In [34], a cooperative localization approach based on projections onto nonconvex boundary sets is proposed for sensor networks, and it is shown that the proposed strategy can achieve better localization performance than the centralized SDP and the distributed MDS although it may get trapped into local minima due to nonconvexity. Similarly, the work in [29] designs a POCS-based distributed positioning algorithm for cooperative networks with a convergence guarantee regardless of the consistency of the formulated CFP, i.e., whether the intersection is nonempty or not.
Although CFPs have attracted a great deal of interest in the literature, QFPs have been investigated only rarely. QFPs represent generalized versions of CFPs in that the constraint sets are constructed from the lower level sets of quasiconvex functions in QFPs whereas such functions are convex in CFPs [32]. The study in [32] explores the convergence properties of subgradient projections based iterative algorithms utilized for the solution of QFPs. It is demonstrated that the iterations converge to a solution of the QFP if the quasiconvex functions satisfy Hölder conditions and the QFP is consistent, i.e., the intersection is nonempty. In this work, we show that the Lambertian model based (originally non-quasiconvex) functions can be approximated by appropriate quasiconvex lower bounds, which convexifies the (originally nonconvex) sublevel constraint sets, thus transforming the formulated feasibility problem into a QFP.

C. Contributions
The previous work on VLP networks has addressed the problem of position estimation based mainly on the ML estimator [13], [16], [40], the least squares estimator [15], [40], triangulation [11], [41], and trilateration [14] methods. This manuscript, however, considers the problem of localization in VLP networks as a feasibility problem and introduces efficient iterative algorithms with convergence guarantees in the consistent case. In addition, the theoretical bounds derived for position estimation are significantly different from those in [16], [40] via the incorporation of terms related to cooperation, which allows for the evaluation of the effects of cooperation on the localization performance in any three dimensional cooperative VLP scenario. Furthermore, unlike the previous research on localization in RF-based wireless networks via CFP modeling [27]- [29], [34], where a common approach is to employ POCS-based iterative algorithms, we formulate the localization problem as a QFP for VLP systems, which necessitates the development of more sophisticated algorithms (e.g., gradient projections) and different techniques for studying the convergence properties of those algorithms (e.g., quasiconvexity and quasi-Fejér convergence). The main contributions of this manuscript can be summarized as follows: • For the first time in the literature, we propose to employ cooperative localization for VLP networks via a generic configuration that allows for an arbitrary construction of connectivity sets and transmitter/receiver orientations. • The CRLB for localization of VLC units is derived in the presence of cooperative measurements (Section II). The effects of cooperation on the performance of localization in VLP systems are illustrated based on the provided CRLB expression (Section VI-A). • The problem of cooperative localization in VLP systems is formulated as a quasiconvex feasibility problem, which circumvents the complexity of the nonconvex ML estimator and facilitates efficient feasibility-seeking algorithms (Section III). • We design gradient projections based low-complexity iterative algorithms to find solutions to the feasibility problem (Section IV). The proposed set-theoretic framework favors the implementation of algorithms in a distributed architecture. • We provide formal convergence proofs for the proposed algorithms in the consistent case based on the concept of quasi-Fejér convergence (Section V).

II. SYSTEM MODEL AND THEORETICAL BOUNDS A. System Model
The proposed cooperative VLP network consists of N L LED transmitters and N V VLC units, as illustrated in Fig. 1. The location of the jth LED transmitter is denoted by y j and its orientation vector is given by n T,j for j ∈ {1, . . . , N L }. The locations and the orientations of the LED transmitters are assumed to be known, which is a reasonable assumption for practical systems [40], [41]. In the proposed system, each VLC unit not only gathers signals from the LED transmitters but also communicates with other VLC units in the system for cooperation purposes. To that aim, the VLC units are equipped with both LEDs and PDs; namely, there exist L i LEDs and K i PDs at the ith VLC unit for i ∈ {1, . . . , N V }. The unknown location of the ith VLC unit is denoted by For the jth PD at the ith VLC unit, the location is given by x i + a i,j and the orientation vector is denoted by n (i) R,j , where j ∈ {1, . . . , K i }. Similarly, for the jth LED at the ith VLC unit, the location is given by x i + b i,j and the orientation vector is represented by n (i) T,j , where j ∈ {1, . . . , L i }. The displacement vectors, a i,j 's and b i,j 's, are known design parameters for the VLC units. Also, the orientation vectors for the LEDs and PDs at the VLC units are assumed to be known since they can be determined by the VLC unit design and by auxiliary sensors (e.g., gyroscope). To distinguish the LED transmitters at known locations from the LEDs at the VLC units, the former are called as the LEDs on the ceiling (as in Fig. 1) in the remainder of the text.
At a given time, each PD can communicate with a subset of all the LEDs in the system. Therefore, the following connectivity sets are defined to specify the connections between the LEDs and the PDs: (2) Namely, S (j) k represents the set of LEDs on the ceiling that are connected to the kth PD at the jth VLC unit. Similarly, S (i,j) k is the set of LEDs at the ith VLC unit that are connected to the kth PD at the jth VLC unit.
The aim is to estimate the unknown locations, x 1 , . . . , x NV , of the VLC units based on the RSS observations (measurements) at the PDs. Let P (j) l,k represent the RSS observation at the kth PD of the jth VLC unit due to the transmission from the lth LED on the ceiling. Similarly, let P (i,j) l,k denote the RSS observation at the kth PD of the jth VLC unit due to the lth LED at the ith VLC unit. Based on the Lambertian formula [12], [42], P where l,k ) is the irradiation angle at the lth LED on the ceiling (at the ith VLC unit) with respect to the kth PD at the jth VLC unit, and θ (j) l,k (θ (i,j) l,k ) is the incidence angle for the kth PD at the jth VLC unit related to the lth LED on the ceiling (at the ith VLC unit). In addition, the noise components, η (j) l,k and η (i,j) l,k , are modeled by zero-mean Gaussian random variables each with a variance of σ 2 j,k . Considering the use of a certain multiplexing scheme (e.g., time division multiplexing among the LEDs at the same VLC unit and on the ceiling, and frequency division multiplexing among the LEDs at different VLC units or on the ceiling), η (j) l,k and η (i,j) l,k are assumed to be independent for all different (j, k) pairs and for all l and i.

B. ML Estimator and CRLB
Let T denote the vector of unknown parameters (which has a size of 3N V × 1) and let P represent a vector consisting of all the measurements in (3) and (4). The elements of P can be expressed as follows: Then, the conditional probability density function (PDF) of P given x, i.e., the likelihood function, can be stated as represents the total number of LEDs that can communicate with the kth PD at the jth VLC unit; that is, N From (7), the maximum likelihood estimator (MLE) is obtained asx and the Fisher information matrix (FIM) [43] is given by Then, the CRLB is stated as wherex represents an unbiased estimator of x. From (7) and (8), the elements of the FIM in (10) can be calculated after some manipulation as (12) Based on (11) and (12), the CRLB for location estimation can be obtained for cooperative VLP systems (please see Appendix F for the partial derivatives in (12)). The obtained CRLB expression is generic for any three-dimensional configuration and covers all possible cooperation scenarios via the definitions of the connectivity sets (see (1) and (2)). To the best of authors' knowledge, such a CRLB expression has not been available in the literature for cooperative VLP systems.
Remark 1: From (12), it is noted that the first summation term in the parentheses is related to the information from the LED transmitters on the ceiling whereas the remaining terms are due to the cooperation among the VLC units. In the noncooperative case, the elements of the FIM are given by the expression in the first line of (12).
Via (11) and (12), the effects of cooperation on the accuracy of VLP systems can be quantified, as investigated in Section VI.

III. COOPERATIVE LOCALIZATION AS A QUASICONVEX
FEASIBILITY PROBLEM In this section, the problem of cooperative localization in VLP networks is investigated in the framework of convex/quasiconvex feasibility. First, the feasibility approach to the localization problem is motivated, and the problem formulation is presented. Then, the convexity analysis is carried out for the resulting constraint sets.

A. Motivation
For the localization of the VLC units, the MLE in (9) has very high computational complexity as it requires a search over a 3N V dimensional space. In addition, the formulation in (9) presents a nonconvex optimization problem; hence, convex optimization tools cannot be employed to obtain the (global) optimal solution of (9). As the number of VLC units increases, centralized approaches obtained as solutions to a given optimization problem (such as (9)) may become computationally prohibitive. Besides scalability issues, centralized methods also require all measurements gathered at the VLC units to be relayed to a central unit for joint processing, which may lead to communication bottlenecks. Therefore, low-complexity algorithms amenable to distributed implementation are needed to efficiently solve the cooperative localization problem in VLP networks. To that aim, the localization problem is cast as a feasibility problem with the purpose of finding a point in a finite dimensional Euclidean space that lies within the intersection of some constraint sets. Feasibility-seeking methods enjoy the advantage of not requiring an objective function, thereby eliminating the concerns for nonconvexity or nondifferentiability of the objective function [44]. Hence, modeling the localization problem as a feasibility problem (i) alleviates the computational burden of minimizing a (possibly nonconvex) cost function in the highly unfavorable centralized setting and (ii) facilitates the use of efficient distributed algorithms involving parallel or sequential processing at individual VLC units.

B. Problem Formulation
Considering the Lambertian formula in (3)-(6), an RSS measurement at a PD can be expressed aŝ P r = P r + η (13) where P r is the true observation (as in (5) or (6)) and η is the measurement noise. Suppose that the RSS measurement errors are negative, which yieldsP r ≤ P r . 1 Then, based on (5) and (6), the following inequality is obtained: where g : R d → R is the Lambertian function with respect to the unknown PD location x, defined as 1 In order to satisfy the negative error assumption, a constant value can always be subtracted from the actual RSS measurement [45]. Decreasing the value of an RSS measurement is equivalent to enlarging the corresponding feasible set. Although this assumption does not have a physical justification, it facilitates theoretical derivations and feasibility modeling of the localization problem. It will be justified via simulations in Section VI-B that the proposed feasibility-seeking algorithms will converge for realistic noise models (e.g., Gaussian), as well. y, n T , n R , and m are known, d is the dimension of the visible light localization network, and γ is given by γ =P r Pt 2π (m+1)A . The field-of-views (FOVs) of the LED transmitters and the PDs are taken as 90 • , which implies that (x − y) T n T ≥ 0 and (y − x) T n R ≥ 0. Under the assumption of negative measurement errors, the feasible set in which the true PD location resides is given by the following lower level set of g(x): which will hereafter be referred to as the Lambertian set. In RF wireless localization networks, such feasible sets are generally obtained as balls [27], [29], hyperplanes [46], or ellipsoids [47], all of which lead to closed-form expressions for orthogonal projection. For k ∈ {1, 2, . . . , K j } and j ∈ {1, 2, . . . , N V }, the Lambertian set corresponding to the kth PD of the jth VLC unit based on the signal received from the ℓth LED on the ceiling for ℓ ∈ S (j) k is defined as follows: ℓ,k is calculated from (3). Similarly, the Lambertian set corresponding to the kth PD of the jth VLC unit based on the signal received from the ℓth LED of the ith VLC unit for where g and γ (i,j) ℓ,k is calculated from (4). The sets defined as in (17) represent noncooperative localization as they are constructed from the RSS measurements corresponding to the LEDs on the ceiling, whereas the sets in (19) are based on the signals from the LEDs of the other VLC units and represent the cooperation among the VLC units. Assuming negatively biased RSS measurements, the problem of cooperative localization in a visible light network reduces to that of finding a point in the intersection of sets as defined in (17) and (19) for each VLC unit. If the Lambertian function in (15) is assumed to be quasiconvex 2 , then the quasiconvex feasibility problem (QFP) can be formulated as follows [32], [48]:

The feasibility problem for cooperative localization of VLC units is given by
where

C. Convexity Analysis of Lambertian Sets
The Lambertian sets as defined in (16) are not convex in general. The following lemma presents the conditions under which the Lambertian sets become convex.
Lemma 1: Consider the α-level set of g ǫ (x), which is given by where ǫ is a small positive constant to avoid non-differentiability and non-continuity of g ǫ (.) at y, as in [49,Eq. 7], k ≥ 1 and γ > 0 are real numbers, and Ω ⊂ R d is defined as Then, B α is convex for each α ∈ R. Proof : Please see Appendix A. Remark 2: Lemma 1 characterizes the type of Lambertian functions whose sublevel sets are convex. Since a function whose all sublevel sets are convex is quasiconvex [50], Lambertian functions of the form (25) are quasiconvex over the halfspace Ω in (26). It can be noted that Ω consists of those VLC unit locations which are able to obtain measurements from an LED located at y due to the receiver FOV limit of 90 • .

D. Convexification of Lambertian Sets
In this part, we utilize Lemma 1 to investigate the following two cases in which the Lambertian functions can be transformed into the form of (25) and Problem 1 becomes a QFP.
1) Case 1: Convexification via Majorization: We propose to approximate the Lambertian function g(x) in (15) by a quasiconvex minorantg(x) such thatg(x) ≤ g(x) for x ∈ Ω and L ⊆L, whereL where (28) is due to the Cauchy-Schwarz inequality and x ∈ Ω, and (29) follows from the unit norm property of the orientation vector. Then, including ǫ in the denominator, we construct the Lambertian sets as (hereafter called expanded Lambertian sets) and Ω being as in (26). According to Lemma 1, L in (30) is convex,g ǫ (x) in (31) is quasiconvex over Ω and the resulting problem of determining a point inside the intersection of such sets turns into a QFP, which can be studied through iterative projection algorithms [32], [48].
2) Case 2: Known VLC Height, Perpendicular LED: In this case, as in [10], [12], [13], [51], it is assumed that the LED transmitters on the ceiling have perpendicular orientations, i.e., T for each j ∈ {1, 2, . . . , N L }, and the height of each VLC unit is known. This assumption is valid for some practical scenarios, an example of which is a VLP network where the LEDs on the ceiling are pointing downwards and the VLC units are attached to robots that move over a two-dimensional plane [7, Fig. 3]. Assuming that the height of the LED transmitters relative to the VLC units is h and T , the Lambertian function in (15) can be rewritten as follows: Then, the Lambertian set corresponding to the function in (32) by introducing ǫ in the denominator is obtained as where Ω is given by (26) andγ = γ/h m . Note that the Lambertian set in (33) is effectively defined on R 2 since the height of the VLC unit is already known. According to Lemma 1, the set defined in (33) is convex. Therefore, in this case, the noncooperative sets as defined in (17) are originally convex. Based on the discussion above, it is concluded that in the case of a known VLC height and perpendicular LED transmitter orientations, the expanded Lambertian sets in (30) defined on R 2 must be used for the measurements among the VLC units (i.e., cooperative measurements) in order to ensure that Problem 1 is a QFP. For the general case in which the LED orientations are arbitrary and/or the heights of the VLC units are unknown, all the noncooperative and cooperative Lambertian sets must be replaced by the corresponding expanded versions in (30).
A noncooperative VLP network is illustrated in Fig. 2(a), where there exist four LED transmitters on the ceiling and two VLC units. In the network, it is assumed that the heights of the VLC units are known and the LEDs on the ceiling have perpendicular orientations so that Case 2 type convex Lambertian sets can be utilized for the measurements between the LEDs on the ceiling and the VLC units. Fig. 2(b) shows the cooperative version of the VLP network with cooperative Lambertian sets including both the nonexpanded (original) sets as in (19) and Case 1 type expanded sets as in (30). It is noted from Fig. 2(b) that incorporating cooperative Lambertian sets into the localization geometry can significantly reduce the region of intersection of the Lambertian sets.

IV. GRADIENT PROJECTIONS ALGORITHMS
In this section, we design iterative subgradient projections based algorithms to solve Problem 1. The idea of using subgradient projections is to approach a convex set defined as a lower contour set of a convex/quasiconvex function by moving in the direction that decreases the value of that function at each iteration, i.e., in the opposite direction of the subgradient of the function at the current iterate [39], [52]. First, the definition of the gradient projector is presented as follows: Definition 1: The gradient projection operator G λ f : R d → R d onto the zero-level set of a continuously differentiable function f : R d → R is given by [53] where λ is the relaxation parameter and ∇ is the gradient operator. The gradient projector can also be expressed as with f + (x) denoting the positive part, i.e., f + (x) = max{0, f (x)}. In the sequel, it is assumed that G λ f (x) = x when x is outside the region where f is quasiconvex.

A. Projection Onto Intersection of Halfspaces
Since the functions of the form (25) are continuously differentiable and quasiconvex on the halfspace Ω in (26), a special case of subgradient projections, namely, gradient projections, can be utilized to solve Problem 1, under the constraint that iterates must be inside Ω to guarantee quasiconvexity. Hence, at the start of each iteration of gradient projections, projections onto the intersection of halfspaces of the form Ω in (26) can be performed to keep the iterates inside the quasiconvex region. The procedure for projection onto the intersection of halfspaces corresponding to the jth VLC unit for j ∈ {1, 2, . . . , N V }, with the halfspaces given bỹ is provided in Algorithm 1 4 . In order to find a point inside the intersection of halfspaces, the method of alternating (cyclic) projections is employed in Algorithm 1, where the current iterate is projected onto each halfspace in a cyclic manner. Convergence properties of this method are well studied in the literature [54], [55]. Γ j is guaranteed to be nonempty since it represents the set of possible locations for the jth VLC unit at which the RSS measurements from the connected LEDs on the ceiling can be acquired. However, the intersection of the halfspaces corresponding to the LEDs of the other VLC units that are connected to the jth VLC unit may be empty due to the VLC unit locations being unknown and variable during iterations.

B. Step Size Selection
An important phase of the proposed projection algorithms is determining the relaxation parameters (i.e., step sizes) associated with the gradient projector. The step size selection procedure exploits the well-known Armijo rule, which is an inexact line search method used extensively for gradient descent methods in the literature [56], [57], [58, Section 1.2]. Algorithm 2 provides an Armijo-like procedure for step size selection given a set of Lambertian functions, the initial step size value λ, a fixed constant β ∈ (0, 1) specifying the degree of decline in the value of the function, step size shrinkage factor ξ ∈ (0, 1), and the current point x. The guarantee of existence of a step size as described in Algorithm 2 can be shown similarly to [59,Lemma 4].

C. Iterative Projection Based Algorithms
In this work, two classes of gradient projections algorithms, namely, sequential (i.e., cyclic) [52] and simultaneous (i.e., parallel) [60] projections, are considered for the QFP described in Problem 1. The proposed algorithm for cyclic projections, namely, the cooperative cyclic gradient projections (CCGP) algorithm, for cooperative localization of VLC units is provided in Algorithm 3. In the proposed cyclic projections, the current iterate, which signifies the location of the given VLC unit, is first projected onto the intersection of halfspaces corresponding to the LEDs on the ceiling via Algorithm 1. Then, the resulting point is projected onto the noncooperative Lambertian set that leads to the highest function value, i.e., the most violated constraint set [32]. Similarly, projection onto the most violated constraint set among the cooperative Lambertian sets is performed and the projections obtained by noncooperative and cooperative sets are weighted to obtain the next iterate.
The cooperative simultaneous gradient projections (CSGP) algorithm is proposed as detailed in Algorithm 4. Simultaneous projections are based on projecting the current point onto each noncooperative and cooperative Lambertian set separately and then averaging all the resulting points to obtain the next iterate. At each iteration, the parallel projection stage is preceded by projection onto the intersection of halfspaces, which aims to ensure that the current iterate resides in the region where all the Lambertian functions corresponding to the fixed anchors (i.e., the LEDs on the ceiling) are quasiconvex. It should be noted that for both cyclic and simultaneous projections, the cooperative Lambertian sets are determined by the latest estimates of the VLC unit locations [29], which are updated in the ascending order of their indices. In addition, the step sizes are updated using the Armijo rule in Algorithm 2.

V. CONVERGENCE ANALYSIS
In this section, the convergence analysis of the proposed algorithms in Algorithm 3 and Algorithm 4 is performed in the consistent case. To that aim, it is assumed that for each j ∈ {1, 2, . . . , N V }, the intersection of the noncooperative and cooperative Lambertian sets in (21) is nonempty; that is, Λ j ∩ Υ j = ∅, where Λ j and Υ j are given by (22) and (23), respectively. In the following, we present the definitions of quasiconvexity and quasi-Fejér convergence, which will be used for the convergence proofs.

Algorithm 3 Cooperative Cyclic Gradient Projections (CCGP)
Initialization: Choose an arbitrary initial point x

Iterative
Step: Given the nth iterate x Projection Onto Intersection of Halfspaces Γ j by Algorithm 1:

Algorithm 4 Cooperative Simultaneous Gradient Projections (CSGP)
Initialization: Choose an arbitrary initial point x

Iterative
Step: Given the nth iterate x Projection Onto Intersection of Halfspaces Γ j by Algorithm 1: Parallel Projection Onto Lambertian Sets: wheren = n for i > j,n = n + 1 for i < j and the weights satisfy ℓ,k ≥ 0, ∀i, ℓ, k. end for Stopping Criterion: Relaxation Parameters: Initialize λ (0) j = λ 0 and update using Algorithm 2 as λ for j ∈ {1, 2, . . . , N V }, where F j and F j are given by (78) and (79) in Appendix B, respectively, and Definition 2 (Quasiconvexity [62]): Definition 3 (Quasi-Fejér Convergence [33]): A sequence {y k } ⊂ R n is quasi-Fejér convergent to a nonempty set V if for each y ∈ V , there exists a non-negative integer M and a sequence {ǫ k } ⊂ R ≥0 such that ∞ k=0 ǫ k < ∞ and For the convergence analysis, we make the following assumptions: A1. Considering any x j ∈ Λ j ∩ Υ j and x j / ∈ Λ j ∩ Υ j , the inequality g i ) holds for every iteration index n and ∀ℓ, k, i, j. A2. The sequence of path lengths taken by the iterations of the proposed algorithms are square summable, i.e., ∞ n=0 Assumption A1 is valid especially when the cooperative algorithms can be initialized at some x = (x 1 , . . . , x NV ) with x j ∈ Λ j , ∀j ∈ {1, 2, . . . , N V }. Assumption A1 implies that any point inside the intersection of the noncooperative and cooperative constraint sets is closer, in terms of the function value (whose zero-level sets are the constraint sets), to the cooperative constraint sets than any point outside the intersection. When the iterations in the cooperative case start from coarse location estimates obtained in the absence of cooperation, the corresponding cooperative sets, which are dynamically changing at each iteration, may involve the set Λ j ∩ Υ j , but exclude the points outside Λ j ∩ Υ j , which yields g On the other hand, Assumption A2 represents a realistic scenario through the Armijo rule in (40) and (41), which ensures a certain level of decline in the Lambertian functions at each iteration and generates a nonincreasing sequence of step sizes.

A. Quasi-Fejér Convergence
In the convergence analysis, the proof of convergence is based on the concept of quasi-Fejér convergent sequences, which possess nice properties that facilitate further investigation, as will be presented in Lemma 2. The following proposition establishes the quasi-Fejér convergence of the sequences generated by Algorithm 4 to the set Λ j ∩ Υ j . As the quasi-Fejér convergence of the sequences generated by the proposed algorithms is stated, the following lemma presents the properties of quasi-Fejér convergent sequences.
Lemma 2 (Theorem 4.1 in [33]): If a sequence {y k } is quasi-Fejér convergent to a nonempty set V , the following conditions hold:

B. Limiting Behavior of Step Size Sequences
In this part, we investigate the limiting behavior of the step size sequences, which are updated according to the procedure in Algorithm 2. The following two lemmas prove that the step size sequences generated by Algorithm 4 and Algorithm 3 have positive limits.
Lemma 3: Any step size sequence λ (n) j generated by Algorithm 4 has a positive limit, i.e., Proof : Please see Appendix D.
Proof : Please see Appendix E. Lemma 3 and Lemma 4 will prove to be useful for deriving the fundamental convergence properties of the proposed algorithms, as investigated next.

C. Main Convergence Results
In this part, we present the main convergence results for the proposed algorithms, i.e., convergence to a solution of Problem 1.
which implies that lim is satisfied ∀f ∈ F j ∪ F j , where the operator J f defined on R d for the set of continuously differentiable functions f : R d → R is given by For a generic Lambertian function in (25), the norm square of the gradient can be expressed as Since the sequence of iterates {x is also bounded, which implies the boundedness of ∇g ǫ (x) . Therefore, based on (59) and (60), it follows that From the Bolzano-Weierstrass Theorem [63, Section 3.4], the boundedness of {x j } ∞ n=0 has a convergent subsequence. Denote this subsequence by {x (nt) j } ∞ t=0 and its limit by x ⋆ j . From (62), it turns out that x ⋆ j ∈ Λ j ∩ Υ j . Therefore, Λ j ∩ Υ j contains a limit point of {x n=0 , which, based on Lemma 2, yields the result that {x (n) j } ∞ n=0 converges to a point inside Λ j ∩Υ j . Based on (50) and the fact that Λ j ⊂ Γ j , it follows that the sequence {x Proposition 4: Let {x (n) } ∞ n=0 be any sequence generated by

, a solution of Problem 1.
Proof 2: Applying similar steps to those in the proof of Proposition 3 and exploiting Proposition 2 and Lemma 4, the following results are obtained: Based on the most violated constraint control in (43) and (44), it is obvious thatf which implies via (63) and (64) that The rest of the proof is the same as that in Proposition 3.
VI. NUMERICAL RESULTS In this section, numerical examples are provided to investigate the theoretical bounds on cooperative localization in VLP networks and to evaluate the performance of the proposed projection-based algorithms. The VLP network parameters are determined in a similar manner to the work in [12] and [13]. The area of each PD is set to 1 cm 2 and the Lambertian order of all the LEDs is selected as m = 1. In addition, the noise variances are calculated using [64,Eq. 6]. The parameters for noise variance calculation are set to be the same as those used in [64] (see Table I in [64]).
The VLP network considered in the simulations is illustrated in Fig. 3. A room of size 10m×10m×5m is considered, where there exist N L = 4 LED transmitters on the ceiling which are located at y 1 = [1 1 5] T m, y 2 = [1 9 5]

A. Theoretical Bounds
In this part, the CRLB expression derived in Section II is investigated to illustrate the effects of cooperation on the localization performance of VLP networks.

1) Performance with Respect to Transmit Power of LEDs on Ceiling:
In order to analyze the localization performance of the VLC units with respect to the transmit powers of the LEDs on the ceiling (equivalently, anchors), individual CRLBs for localization of the VLC units in noncooperative and cooperative scenarios are plotted against the transmit powers of LEDs on the ceiling in Fig. 4, where the transmit powers of the VLC units are fixed to 1 W. As observed from Fig. 4, cooperation among VLC units can provide substantial improvements in localization accuracy (about 81 cm and 37 cm improvement, respectively, for VLC 1 and VLC 2 for the LED transmit power of 300 mW). We note that the improvement gained by employing cooperation is higher for VLC 1 as compared to that for VLC 2. This is an intuitive result since the localization of VLC 1 depends mostly on LED 1 and LED 2 (the other LEDs are not sufficiently close to facilitate the localization process), and incorporating cooperative measurements for VLC 1 provides an enhancement in localization performance that is much greater than that for VLC 2, which can obtain informative measurements from the LEDs on the ceiling even in the absence of cooperation as seen from the network geometry in Fig. 3. In addition, the CRLBs in the cooperative scenario converge to those in the noncooperative scenario as the transmit powers of the LEDs increase. Since the first (second) summand in the FIM expression in (12) corresponds to the noncooperative (cooperative) localization, higher transmit powers of the LEDs on the ceiling cause the first summand to be much greater than the second summand, which makes the contribution of cooperation to the FIM negligible. Hence, the effect of cooperation on localization performance becomes more significant as the transmit power decreases, which is in compliance with the results obtained for RF based cooperative localization networks [20].
2) Performance with Respect to Transmit Power of VLC Units: Secondly, the localization performance of the VLC units is investigated with respect to the transmit powers of the VLC units when the transmit powers of the LEDs on the ceiling are fixed to 1 W. Fig. 5 illustrates the CRLBs for localization of the VLC units against the transmit powers of the VLC units in the noncooperative and cooperative cases. As observed from Fig. 5, cooperation leads to a higher improvement in the performance of VLC 1, similar to Fig. 4. In addition, via the FIM expression in (12), it can be noted that the contribution of cooperation to localization performance gets higher as the transmit powers of the VLC units increase, which is also observed from Fig. 5. However, the CRLB reaches a saturation level above a certain power threshold, as opposed to Fig. 4, where the CRLB continues to decrease as the power increases. The main reason for this distinction between the effects of the transmit powers of the LEDs on the ceiling and those of the VLC units can be explained as follows: For a fixed transmit power of the VLC units, the localization error by using three anchors (i.e., three LEDs on the ceiling that are connected to the corresponding VLC unit) converges to zero as the transmit powers of the anchors increase regardless of the existence of cooperation. On the other hand, for a fixed transmit power of the LEDs on the ceiling, increasing the transmit power of the VLC unit (i.e., one of the anchors) cannot reduce the localization error below a certain level. Therefore, the saturation level represents the localization accuracy that can be attained by four anchors with three anchors leading to noisy RSS measurements and one anchor generating noise-free RSS measurements.

B. Performance of the Proposed Algorithms
In this part, the proposed algorithms in Algorithm 3 (CCGP) and Algorithm 4 (CSGP) are evaluated in terms of localization performance and convergence speed. For both algorithms, the initial step size is selected as λ 0 = 1, the step size shrinkage factor and the degree of decline in the Armijo rule in Algorithm 2 are set to ξ = 0.5 and β = 0.001, respectively. The VLC units are initialized at the positions of the closest LEDs on the ceiling which are connected to the corresponding VLC units.
Localization performances of the algorithms are presented in both the absence and the presence of cooperation and compared against those of the ML estimator in (9) and the CRLBs derived in Section II. In order to ensure convergence to the global minimum, the ML estimator is implemented using a multistart optimization algorithm with 100 initial points randomly selected from the interval [0 10] m at each axis. 5 In addition, two different measurement noise distributions, namely, Gaussian and exponential, are considered while evaluating the proposed algorithms as in [28]. The Gaussian noise is used to model the case in which the RSS measurement noise can be both positive and negative, whereas the exponentially distributed noise (subtracted from the true value) represents the scenario in which the RSS measurements are negatively biased, which leads to the feasibility modeling of the localization problem in Section III-B. Furthermore, the average residuals at each iteration are calculated to assess the convergence speed of the proposed algorithms [29]: In the simulations, two-dimensional localization is performed by assuming that the VLC units have known heights. Therefore, with the knowledge of perpendicular LED orientations, Case 2 type Lambertian sets in Section III-D2 are utilized for localization based on the measurements from the LEDs on the ceiling. The cooperation among the VLC units is modeled by Case 1 type Lambertian sets in Section III-D1.
1) Gaussian Noise: In Fig. 6, the average localization errors of the VLC units for the different algorithms are plotted against the transmit power of the LEDs on the ceiling for the case of the Gaussian measurement noise by fixing the transmit powers of the LEDs at the VLC units to 1 W. From Fig. 6, it is observed that the cooperative approach can significantly reduce the localization errors, especially in the low SNR regime (about 60 cm and 70 cm reduction for CSGP and CCGP algorithms, respectively, for 100 mW LED transmit power). In addition, both Algorithm 3 (CCGP) and Algorithm 4 (CSGP) can attain the localization error levels that asymptotically converge to zero at the same rate as that of the CRLB. Moreover, it can   be inferred from Fig. 6 that the proposed iterative methods achieve higher localization performance than the ML estimator in the low SNR regime for both the noncooperative and the cooperative scenarios. Although the ML estimator is forced to converge to the global minimum via the multi-start optimization procedure involving 100 different executions of a convex optimization solver, whose complexity may be prohibitive for practical implementations, it has lower performance than the proposed approaches, which depend on low-complexity iterative gradient projections. Hence, at low SNRs, the proposed algorithms are superior to the MLE in terms of both the localization performance and the computational complexity. Furthermore, the simultaneous projections outperforms the cyclic projections at low SNRs at the cost of a higher number of set projections, but the two approaches converge asymptotically as the SNR increases. Fig. 7(a) and Fig. 7(b) report the average residuals calculated by (68) corresponding to the proposed algorithms versus the number of iterations for 100 mW and 1 W of transmit powers of the LEDs on the ceiling, respectively. CSGP in the absence of cooperation has the fastest convergence rate and exhibits an almost monotonic convergence behavior. However, CSGP in the cooperative scenario shows relatively slow convergence in general and a locally nonmonotonic behavior when several consecutive iterations are taken into account. This is due to the cooperative Lambertian sets being involved in the simultaneous projection operations. In addition, cyclic projections tend to settle into limit cycle oscillations after few iterations, thus implying that the sequence itself does not converge to a point, but it has several subsequences that converge [39]. This behavior, called cyclic convergence, is encountered in cyclic (sequential) projections if the feasibility problem is inconsistent [39], [65]. Furthermore, by comparing Fig. 7(a) and Fig. 7(b), it is observed that the magnitude of limit cycle oscillations gets smaller as the SNR increases since the region of uncertainty becomes narrower at higher SNR values, thereby making the convergent subsequences close to each other.
2) Exponential Noise: To investigate the performance of the algorithms under exponentially distributed measurement noise, the average localization errors are plotted against the transmit power of the LEDs on the ceiling for the case of the subtractive exponential noise in Fig. 8. Similar to the case of the Gaussian noise, the proposed algorithms succeed in converging to the true VLC unit positions as the SNR increases. Since the projection based methods rely on the assumption of negatively biased measurements, they perform slightly better at low SNRs as compared to the case of the Gaussian noise. On the other hand, the MLE produces larger errors at low SNRs for the exponentially distributed noise since its derivation is based on the assumption of Gaussian noise. The average residuals in the case of the exponentially distributed noise are illustrated in Fig. 9(a) and Fig. 9(b) for two different LED power levels. In contrary to the case of Gaussian noise, cyclic projections do not fall into limit cycles and provide globally monotonic convergence results as the feasibility problem is consistent, which complies with the results presented in the literature pertaining to the study of CFPs [39]. In addition, it is observed that both the cyclic and the sequential projection methods have faster convergence for lower SNR values since it takes fewer iterations to get inside the intersection of the constraint sets, which becomes larger as the SNR decreases.
VII. CONCLUDING REMARKS In this manuscript, a cooperative VLP network has been proposed based on a generic system model consisting of LED transmitters at known locations and VLC units with multiple LEDs and PDs. First, the CRLB on the overall localization error of the VLC units has been derived to quantify the effects of cooperation on the localization accuracy of VLP networks. Then, due to the nonconvex nature of the corresponding ML expression, the problem of cooperative localization has been formulated as a QFP, which facilitates the development of lowcomplexity decentralized feasibility-seeking methods. In order to solve the feasibility problem, iterative gradient projections   based algorithms have been proposed. Furthermore, based on the notion of quasi-Fejér convergent sequences, formal convergence proofs have been provided for the proposed algorithms in the consistent case. Finally, numerical examples have been presented to illustrate the significance of cooperation in VLP networks and to investigate the performance of the proposed algorithms in terms of localization accuracy and convergence speed. It has been verified that the proposed iterative methods asymptotically converge to the true positions of VLC units at high SNR and exhibit superior performance over the ML estimator at low SNRs in terms of both implementation complexity and localization accuracy. An important research direction for future studies is to explore the convergence properties of Algorithm 3 and Algorithm 4 when the proposed QFP is inconsistent. In the inconsistent case, simultaneous projection algorithms tend to converge to a minimizer of a proximity function that specifies the distance to constraint sets [39], [66]. For the implicit CFP (ICFP) considered in TOA-based wireless network localization, the POCS based simultaneous algorithm is shown to converge to the minimizer of a convex function, which is the sum of squares of the distances to the constraint sets [29]. Therefore, finding proximity functions characterizing the behavior of simultaneous projections (e.g., Algorithm 4) for the inconsistent QFPs [32] would be a significant extension for the set-theoretic estimation literature.

B. Proof of Proposition 1
Since Λ j ∩ Υ j = ∅, consider any point x j ∈ Λ j ∩ Υ j . At the nth iteration, it can be assumed that x (n) j / ∈ Λ j ∩ Υ j because otherwise iterations will stop via (51) and (35), which implies quasi-Fejér convergence of {x (n) j } ∞ n=0 to Λ j ∩ Υ j based on Definition 3. Then, based on the iterative step in (51), the following is obtained: