On the Achievable Max-Min User Rates in Multi-Carrier Centralized NOMA-VLC Networks

Visible light communications (VLC) is gaining interest as one of the enablers of short-distance, high-data-rate applications, in future beyond 5G networks. Moreover, non-orthogonal multiple-access (NOMA)-enabled schemes have recently emerged as a promising multiple-access scheme for these networks that would allow realization of the target spectral efficiency and user fairness requirements. The integration of NOMA in the widely adopted orthogonal frequency-division multiplexing (OFDM)-based VLC networks would require an optimal resource allocation for the pair or the cluster of users sharing the same subcarrier(s). In this paper, the max-min rate of a multi-cell indoor centralized VLC network is maximized through optimizing user pairing, subcarrier allocation, and power allocation. The joint complex optimization problem is tackled using a low-complexity solution. At first, the user pairing is assumed to follow the divide-and-next-largest-difference user-pairing algorithm (D-NLUPA) that can ensure fairness among the different clusters. Then, subcarrier allocation and power allocation are solved iteratively through both the Simulated Annealing (SA) meta-heuristic algorithm and the bisection method. The obtained results quantify the achievable max-min user rates for the different relevant variants of NOMA-enabled schemes and shed new light on both the performance and design of multi-user multi-carrier NOMA-enabled centralized VLC networks.


Introduction
Wireless optical communications, particularly visible light communications (VLC), has emerged as a bandwidth-abundant, secure, and cost-effective communications technology. It complements the existing radio frequency (RF) systems or even replaces them for some indoor and outdoor applications, such as conference and exhibition halls, office rooms, trains and airplane cabins, and so forth. It can also be deployed outdoors in vehicle-tovehicle (V2V) and vehicle-to-everything (V2X) applications, as well as applications which are short-range, and which have a very high data rate in beyond-fifth generation (B5G) networks. VLC is based on the principle of modulating light from light emitting diodes (LEDs) without any adverse effects on the human eye and at the required illumination levels to transmit data. Clearly, this provides an excellent opportunity to utilize the existing illumination infrastructure for very high-speed and secure wireless communications [1].
In a conventional distributed VLC architecture, a set of L transmit LEDs communicate with a set of N users in a way that each light emitting diode (LED) acts as an access point (AP) that serves its own subset of the N users. Besides, all the access points (APs) contain a base-band unit followed by an optical front-end and connect to each other through the data backbone, as well as the electrical grid. While in centralized light access network adopted in this work, as it is attractive for the max-min user rate requirement [9]. The reason is that by re-distributing the rate gains between the clusters and achieving cluster fairness, the rate gap between the users' rates will reduce and this reduces the search space, and thus the algorithm complexity, for reaching the max-min user rate point for all the users in the network.
The subcarrier allocation problem is typically a non-deterministic polynomial-time (NP)-hard optimization problem in both the number of users and the number of channels to be allocated [12]. The optimal resource allocation problem is known to be an NP-hard problem even for the single carrier case [12]. The NP class contains problems that have the characteristic that only non-deterministic algorithms can solve them in polynomialtime. However, we can verify a given solution to these problems in polynomial time. The non-deterministic algorithms solve a problem by evaluating possible guesses of the solution. A problem is NP-hard if it is at least as hard as the hardest problem in NP, and all NP problems are reducible to it. In NP-hard problems, the search space size is so large that an exhaustive search is infeasible. Therefore, we determine approximate solutions to NP-hard problems. Evolutionary algorithms (EAs) are a popular approach to solve NPhard problems [13]. Simulated Annealing (SA) is a well-known EA that has demonstrated excellent performance in solving optimization problems in various applications, including wireless communications [13][14][15][16]. In this work, we use the SA algorithm to allocate the subcarriers to user-pairs. For the power allocation, the bisection method is commonly used to optimally solve the power allocation problem for max-min user rate optimization in the NOMA literature [17,18].

Related Work
When it comes to rate-optimal NOMA-enabled VLC networks, the majority of works have considered maximizing the system sum-rate metric [19][20][21][22][23][24][25][26], and this is detailed in Table VIII of the survey by Maraqa et al. [7]. In particular, early works on rate-optimal NOMA-enabled VLC networks have considered a single-LED, single-carrier setup as in [19][20][21][22][23]. Then, the research evolved to consider a single-LED multi-carrier setup as in [24,25]. Later, a few works have considered a multi-LED, multi-carrier setup as in [26][27][28]. Among those, only two works have considered the optimization of another common fairness measure, which is the max-min achievable user rate alongside maximizing the system sum-rate metric [27,28]. This is referred to as the reason that the optimization of the max-min user rate tends to reduce the total throughput of the users, as compared to sumrate maximization; however, it guarantees fairness among the users. This is an appealing feature for the operation of ultra-massive machine-type communication (umMTC) and internet-of-everything (IoE) in B5G systems that need reliable connectivity with generally low data-rate requirements [29]. The work in [27] has considered the maximization of the minimum rate for NOMA-enabled VLC networks, where an approximation of the objective function was utilized to convexify the optimization problem, and a gradient projection (GP) algorithm is used to adjust the power allocation coefficients. The maximization of the minimum user rate for the cell-edge users in NOMA-enabled VLC networks in [28], using hybrid NOMA and/or linear zero-forcing (ZF) pre-coding schemes, was approached by numerical convex optimization of the standard determinant maximization. However, both aforementioned works have considered the optimization of both the power allocation and user grouping. Neither of them has analyzed the effect of subcarrier allocation on the max-min user rate performance of NOMA-enabled VLC networks, which is the main objective of this work.
Our contributions can thus be summarized as follows: • The optimization of the achievable max-min user rates for NOMA-enabled centralized VLC is investigated through formulating a joint problem for the user pairing, the subcarrier allocation, and the power allocation. Then, a low complexity solution is proposed. • The development of Simulated Annealing (SA)-assisted algorithm for tackling the subcarrier allocation in the max-min user rate optimization problem. The obtained results are further verified using the Tabu-search (TS) algorithm. • The implementation of both of the NOMA-imposed schemes, where all the users are grouped into pairs, and the NOMA-not-imposed schemes; and besides these, the investigation of the effect of the different network parameters on the achievable max-min user rate.
The remainder of the paper is organized as follows. In Section 2, we provide the system and channel models for the utilized C-LiAN architecture. The formulation of the proposed max-min user rate optimization problem is analyzed in Section 3. In Section 4, the details of the proposed heuristic-based solution for the formulated max-min user rate optimization are explained. Finally, simulation results, paper conclusions, and future research directions are given in Section 5, Section 6 and Section 7, respectively. Table 1 summarizes the notations used in this paper.
The ratio between the electrical signal power and the optical transmit power κ The optical to electrical conversion efficiency of the photo-diodes (PDs) The power of equivalent AWGN, where Z o denotes the noise power spectral density and B L denotes the base-band modulation bandwidth The achievable rate of an arbitrary user (j-th user) of an arbitrary LED in the network (i.e., R j can be a strong user or a weak user in an arbitrary user pair) A binary variable to denote that a user is served by a LED l i and a subcarrier k S k ∈ {0, 1} A binary variable to denote that a user is served by a subcarrier k S j l i ∈ {0, 1} A binary variable to denote the user u j is served by LED l i Γ(l i ) A set of all user-pair combinations for an arbitrary LED l i K l i The maximum number of subcarriers for an arbitrary LED l i

Heuristic-based Solution's Variables
The most suitable LED for a user u j ∈ N f 1 (l i ) A function that returns the number of users assigned to the LED l i f 2 (u j ) A function that denotes the distance of the user u j from the LED to which it is currently allocated A function that denotes the maximum distance of u j from any LED N l i , N l i A vector that store the users of LED l i after binding, the number of users in N l i

System and Channel Models
A multi-user NOMA-enabled C-LiAN architecture [2] is shown in Figure 1. In this paper, we assume that the ceiling of the room contains uniformly distributed LEDs denoted by L = {l 0 , l 1 , ..., l L−1 }, where L is the number of LEDs. Each LED is assigned up to K subcarriers that can be shared by several users, where the subcarriers of each LED are denoted as K = {s 0 , s 1 , .., s K−1 }. In addition, the users are randomly located in the room, and denoted by U = {u 0 , u 1 , ..., u N−1 }; where the total number of users in the network is N. In multi-user systems, the transmission resources described in time, wavelength, and/or space are divided into resource units. A key challenge is the allocation of the resource units in a way that key performance metrics, such as user fairness and spectral efficiency, are fulfilled. This motivates allowing the served users or a subset of them to share the subcarriers or a subset of them. Hence, in each signaling interval, we assume that each pair of users can occupy one or more subcarrier(s) in only one LED. In this work, the line-of-sight (LoS) optical channel model is adopted, since in a typical indoor environment, the strongest diffuse component is usually much weaker than the LoS component [30]. Hence, the channel gain between the LED l 0 and user u s is given by [16]: where ϕ l 0 ,u s and ψ l 0 ,u s denote, respectively, the angle of irradiance and angle of incidence between the LED l 0 and user u s . As shown in Figure 2, d l 0 ,u s is the distance between the LED l 0 and user u s . Ψ 1/2 is the FoV semi-angle of the user u s . m is the order of Lambertian emission and is equal to −1 log 2 (cos(φ 1/2 )) , where φ 1/2 is the semi-angle of the LED l 0 . A p is the area of the PD for user u s , T s (ψ l 0 ,u s ) is the optical filter gain, and χ denotes the refractive index.
We assume that served users by LED l 0 are ordered according to their channel gains in ascending order as |h l 0 ,u 0 | ≤ |h l 0 ,u 1 | ≤ ... ≤ |h l 0 ,u N 0 −1 |; where we assume here that the total number of served users by LED l 0 is N 0 . In general, if the user u 0 (i.e., the user that has the weakest instantaneous channel gain) of LED l 0 occupies an arbitrary subcarrier (for example, the k-th subcarrier), then the signal-to-noise-interference ratio (SINR) of the user u 0 can be expressed as: where a 0 represents the power allocation factors for the user u 0 that is attached to LED l 0 . P o,k is the LED optical transmit power per subcarrier and is equal to P o /(K − 2) according to the DCO-OFDM principle, where P o denotes the optical transmit power at the output of the LED. The ratio between the electrical signal power and the optical transmit power can be expressed as ι = P o √ P e , where P e denotes the electrical signal power. κ denotes the optical-to-electrical conversion efficiency of the PDs, σ 2 k is the power of the equivalent additive Gaussian noise (practical indoor VLC networks are affected by both thermal noise due to the receiver pre-amplifier, and ambient light shot noise due to the possible sun-light through the windows and/or non-VLC indoor light sources. The total noise can be modeled as signal-independent Gaussian noise [31,32]), whose variance is the sum of the variances of these two noise components and equals to Z o B L /K; where Z o denotes the noise power spectral density. In the simulation section of this work, the value of Z o accounts for both thermal noise and light shot noise. B L denotes the base-band modulation bandwidth. In Equation (2), the first term in the denominator, represents the interference induced from the users that utilize the k-th subcarrier in the interfering LEDs (i.e., inter-LED interference) where S l i ,k ∈ 0, 1 is a binary variable that denotes the user has received an interfering signal from a LED l i on a subcarrier k. Likewise, the second term in the denominator, o,k , represents the interference that remains after successive interference cancellation (SIC) decoding (i.e., inter-user interference) where a j represents the power allocation factors for the rest of the users that is attached to LED l 0 . It is worthy to note that we assume in the above equation that multiple users can be served by one subcarrier of LED l 0 (i.e., multiple users per cluster). For the rest of the paper, we consider the case where two users are served by the same subcarrier(s) (i.e., two users per cluster). Studying the system performance while considering multiple users per cluster is left for future research. For the two-user pairing case, we use the divide-and-next-largest-difference userpairing algorithm (D-NLUPA) scheme for pairing the users, as mentioned in Section 1, where the strong user (i.e., user u s ) is close to the LED with strong channel gain and the weak user (i.e., user u w ) is farther away from the LED with weak channel gain so that the strong user can successfully decode and subtract the weak user signal before decoding its own signal (i.e., performing SIC). In Figure 1, for LED l 0 , we can express the SINR of the strong user, the SINR of the weak, and the achievable rates of the strong and the weak users using a well-known lower bound for the capacity, given in [33,34], as where S k ∈ 0, 1 is a binary variable to denote that the user is served by a subcarrier k, and it is assumed that the two-user pair (user u s and user u w ) occupies one or more subcarrier(s) of the LED l 0 and the power allocation factors satisfy the condition a s + a w = 1. Noting that the analysis in Equations (3) and (4) represents the rate equations for one user-pair of LED l 0 , the rate equations for the rest of the user-pairs of LED l 0 and for the user-pairs of any other LED in the network can be analyzed similarly.

The Max-Min User Rate Optimization Problem
In this section, we formulate a joint optimization problem for user pairing, subcarrier allocation, and power allocation that maximizes the minimum achievable rate of the served users. Let R j denote the achievable rate of an arbitrary user (i.e, j-th user) of an arbitrary LED (i.e, LED l i ) in the network (i.e., R j can be a strong user or a weak user in a user pair, as described in Equations (3) and (4)). Then, the max-min optimization problem can be expressed as: subject to where Γ(l i ) is a set of all user-pair combinations of an arbitrary LED l i . S j k ∈ {0, 1} and S j l i ∈ {0, 1} are binary variables to denote that the user u j is served by a subcarrier k and the user u j is served by LED l i , respectively. By Equation (6), we ensure that any user attached to LED l i is served by at least one subcarrier, noting that the weak user and the strong user in a user-pair are served by the same number of subcarriers according to the NOMA principle. By Equation (7), the maximum number of subcarriers per LED is limited to K 2 − 1 to be aligned with the DCO-OFDM principle. By Equation (8), we ensure that each user attached to LED l i is served only by LED l i . By Equation (9), we limit the values of the power allocation coefficients in a user-pair. The above optimization problem is of a combinatorial nature, for the user pairing and the subcarrier allocation, and has a non-linear objective function, and hence a non-convex NP-hard optimization problem [12,35]. The optimal values of the power allocation coefficients can be determined using the bisection method [17].

The Heuristic-Based Solution for the Max-Min User Rate Optimization Problem
The optimization problem in Equations (5)-(9) is an involved optimization problem, as stated before. Subsequently, the search space to handle the problem jointly is too large. Noting that performing an exhaustive search to solve this problem jointly is infeasible from a practical point of view, to reduce the complexity of this optimization problem, the user pairing is assumed to follow the D-NLUPA algorithm, as mentioned before, where that pairing algorithm is attractive for the max-min user rate requirement [9]. Hence, the optimization problem in Equations (5)-(9) reduces to a joint subcarrier allocation and power allocation problem, which is still an NP-hard problem because of the remainder of the subcarrier allocation problem [12,35,36], though with less complexity as compared to the original optimization problem. Consequently, the proposed heuristic-based solution, that is discussed in this section, tackles the original optimization problem in three steps, as in Algorithm 1:

1.
Binding of users to LEDs.

3.
Optimizing subcarrier(s) allocation to user-pairs in each LED and power allocation within each pair (i.e., subcarrier allocation and power allocation).
As stated before, the max-min optimization problem in Equations (5)-(9) involves both subcarrier and power allocation; however, since each of the NOMA users in each pair will share the same assigned subcarriers with a fixed power per subcarrier, then the subcarrier allocation and the power allocation problems may be done in an iterative manner through both the SA algorithm and the bisection method. This is discussed in detail later in this section.

Algorithm 1:
Overview of the proposed heuristic-based solution.
1 Binding of users to LEDs: Assign users to different LEDs in such a way that a user cannot be assigned to more than one LED; 2 Determining the user-pairs for each LED: Apply the D-NLUPA algorithm to determine the user-pairs per LED. Store the user-pairs in Γ(l i ); 3 Optimizing subcarrier(s) allocation to user-pairs in each LED and power allocation within each pair: This task allocates the subcarriers to the user-pairs to maximize the minimum achievable rate of the served users. We apply the SA algorithm with an innovative objective function that allocates the users-pairs to the subcarriers while ensuring an almost uniform data-rates for all the served users.

Binding of Users to LEDs
Although a user can receive data using any number of subcarriers, all subcarriers should belong to a single LED. Therefore, it is necessary to bind users to LEDs. In the VLC literature, there are different approaches to attach users to LEDs [28]. In this work, we adopted the strategy to bind the users to the LEDs that have the channels of maximum strength with the users. We denote the most suitable LED for a user u j ∈ N with λ(u j ), and its value can be determined using the following equation: The imposition is that the NOMA scheme, termed in this work as the NOMA-imposedscheme (i.e., two users should share the available subcarrier(s)), requires that the number of users allocated to each LED should be even. Otherwise, it is not possible to impose NOMA fully, and as a result, some subcarriers are allocated to only one user (this scheme is termed in this work as the "NOMA-not-imposed scheme"). When the number of users allocated to any LED is odd, we can formulate a linear constrained optimization problem and solve it to ensure that each LED has an even number of users or zero users. Next, we describe the optimization problem. In Equation (11), |.| denotes a complement operation, and this operation returns one only when |λ(u j ) − i| = 0. The function f 1 (l i ) returns the number of users assigned to the LED l i .
The distance of the users from their allocated LEDs can be given by: In Equations (12) and (13), f 2 (u j ) denotes the distance of the user u j from the LED to which it is currently allocated, and f 3 (u j ) denotes the maximum distance of u j from any LED. The linear constrained optimization problem is shown below: subject to The above linear-constrained optimization problem can be solved using many methods, and in our work, we employed an iterative greedy algorithm [37] to solve it. The solution is represented by {λ(u 0 ), λ(u 1 ), ..., λ(u N−1 )}, where λ(u j ) ∈ {l 0 , l 1 , ..., l L−1 }. Algorithm 2 summarizes the aforementioned approach followed for binding users to LEDs. The pseudo-code outlines an iterative method whose termination criterion could be the maximum number of iterations (e.g., 1000 iterations). The iterative method makes a random change in the current binding of users and only lets the new solution replace the current one if the new binding is better than or equal to the current one using Equations (14) and (15). The simulation results, shown in Section 5.1, illustrate that the above-described method is efficient in preventing the LEDs from having an odd number of users with a relatively low number of iterations.

Determining of the User-Pairs for Each LED
We used the D-NLUPA method [9] to determine the user-pairs for each LED. The D-NLUPA method is efficient in finding user-pairs, and it creates pairs in which the first user has a strong channel (or strongest among users not already part of any pair) with a user of a relatively weaker channel. Algorithm 3 shows the D-NLUPA algorithm. The first step is to sort the users with respect to their channel gains. The second step is a loop that creates pairs of strong and weak users and stores them in Γ(l i ). The actions on lines 6-7 only execute if the NOMA-not-imposed scheme is deployed and the number of users under any LED is odd. We should place the last user (i.e., the user having minimum channel strength) into a pair whose weak user is null. Insert P j in Γ(l i ); 5 end 6 if NOMA-not-imposed scheme==TRUE && N l i is an odd number then 7 Create a pair in which the user at index N l i − 1 is the strong user and a null value for the weak user, and insert it into Γ(l i ); 8 end 9 return Γ(l i );

Optimizing Subcarrier(s) Allocation to User-Pairs in Each LED and Power Allocation within Each Pair
In this subsection, we discuss the problem of allocating the subcarriers to user-pairs, as well as a solution method. As we mentioned in the previous subsection, for each LED, the user-pairs allocated to it are represented by Γ(l i ), and each LED has up to K subcarriers. The allocation aims to maximize, as well as keep uniformness in the data-rates of the users. We also propose a novel objective function to achieve the dual goals through applying a single function. In the following, we first discuss our proposed objective function, and then discuss a solution method using the SA algorithm. We represent the solution of subcarrier allocation using the following notation: x 0,0 , x 0,0 , ..., x 0,K−1 x 1,0 , x 1,1 , ..., x 1,K−1 . . . , . . . , . . .
where x i,k is an integer, and k ∈ {0, ..., K − 1}, i ∈ {0, ..., L − 1}, and x i,k denote the index of the user-pair of Γ(l i ) assigned to the k-th subcarrier of the LED l i . The objective of the subcarrier allocation is to maximize the minimum achievable rate of the served users through finding a decision matrix (X), and is given by maximize O(X).
We propose two new constraints into the above objective function (O(X)) to help the iterative heuristics to obtain a solution in which the data-rates of all users have the following features: 1.
The data-rates of all users are closer to each other.

2.
The users should not have a zero data-rate. 3.
The data-rate should be as maximal as possible considering the above two conditions.
To include the above enhancements, we propose the objective function as follows: In the above equations, Equation (20) ensures that the percentage difference between the maximum and minimum data-rates of any two users should be lesser than a predefined constant C ∈ [0, 1]. Equation (21) ensures that no user suffers a service outage.
Finally, we convert the above constrained problem into an unconstrained optimization problem using the penalty method [38,39], as follows: The new objective function can be given by: maximize O (X).
In the above equation, P 1 and P 2 denote the penalty factors. The function O(X) denotes a data-rate value of a user-pair, and hence it could be zero or positive, and the rates can vary between different user pairs. Therefore, P 1 is set for the avoidance of solutions with zero minimum data-rates (i.e., f cons > 0), P 2 is set for the avoidance of solutions in which the difference among the data-rates of the users is more than the given value (i.e., f diff > 0). We prioritize the two constraints by setting P 1 >> P 2 . The function (O (X)) is a non-linear integer programming (NLIP) problem, which is NP-hard in general [40]. In this work, we adopted the SA algorithm, which is a popular method of solving NP-hard optimization problems [14][15][16]. Now, we briefly describe the SA algorithm for finding an allocation of subcarriers to user-pairs in order to maximize the objective function Equation (28). Algorithm 4 shows an overview of the SA algorithm. The input parameters are as follows: (i) T 0 is the initial temperature of the SA algorithm, and its value should be high; (ii) α is the cooling-rate and its value should be between (0, 1), usually, its value is kept very close to 1, such as 0.97, 0.99; (iii) M indicates the number of iterations in the Metropolis function; (iv) β indicates the increase in the value of M. The first four lines in the SA algorithm initialize the variable and the current solution (X current ). The 'while' loop contains the main algorithm, and the termination criterion could be the maximum time or iterations. The 'while' loop in each iteration calls the Metropolis function, which is responsible for exploring the current solution's neighborhood to improve the current solution (X current ).
Algorithm 5 shows the Metropolis function. The Neighbor function in line 2 creates a new solution by mutating the current solution. In this work, we employ the random mutation that consists of the following two steps: (i) Randomly choose a LED and a subcarrier of it; and (ii) change the user-pair already allocated to the selected subcarrier to another user-pair attached to the same LED. The new solution (X new ) could replace the current solution based on comparing its objective function values with that of the current solution. The RANDOM in the pseudo-code indicates a random number between (0, 1). The readers can refer to [13] for more details on the SA algorithm. For the optimal power allocation, the bisection method has been utilized, which can be found in [17].

Algorithm 4:
Overview of the SA algorithm.
Input: T 0 , α, β, M Output: X 1 Initialize X to a random solution 2 T = T 0 , X current = X, X best = X 3 c current = O (X) 4 c best = c current 5 while stopping criterion is not reached do 6 Call Metropolis(X current , c current , X best , c best , T, M) 7 Replace M by β × M, and T by α × T 8 end 9 return X best

The Complexity Analysis of the Proposed Heuristic-Based Solution
In this section, we briefly discuss the time complexity of the algorithms used to accomplish different tasks of the proposed heuristic-based solution. The binding method, discussed in Section 4.1, requires the determination of the maximum value in an array which can be implemented using max-heap or linear determination [37]. The time complexity of the operation to find the maximum value is O(L), and since the operation should be repeated for each user, therefore, the total complexity of the binding method is O(NL). A brief analysis of the time complexity of the algorithm to implement the D-NLUPA method, discussed in Section 4.2, is as follows. The sorting can be accomplished using heapsort [37] that has a complexity of O(NlogN), and the creation of user-pairs is a linear operation that has complexity equal to O(N). The time complexity of this step is O (LNlogN) because the heap-sort should be performed for each LED. The computation for each LED is independent of the others. Therefore, we can perform parallel computations, and in that case, the complexity reduces to O(NlogN). We optimized the allocation of user-pairs to the subcarriers using the SA algorithm. The complexity of this algorithm is usually expressed in terms of the number of iterations needed for it to converge to its best solution [41]. In the Results and Discussions section, we analyzed the convergence aspects of the SA algorithm in detail.

Results and Discussions
In this section, the obtained results for the achievable max-min user rates in a C-LiAN architecture are presented. The user pairing, the sub-carrier allocation, and the power allocation were carried out using the D-NLUPA, the proposed SA algorithm, and the bisection method, respectively, as detailed before. To assess the performance of the proposed heuristic-based solution, the Monte-Carlo simulation that is averaged over 100 different users' location realizations is used; specifically, each point in the performance curves is an average of implementations of 100 different users' location realizations. The proposed heuristic-based solution was implemented in R and C++ using R version 4.0.2, RStudio version 1.3.1, and Rcpp. We ran the simulations on a desktop computer with an Intel 2.6 GHz processor and 128 GB of memory. In the simulations, an empty room was considered, with a dimension of 5 × 5 × 3, unless otherwise stated. The room was equipped with four LEDs, unless otherwise stated, mounted in the ceiling of the room with an equal distance between adjacent LEDs, in a square lattice topology, similar to [16,42]. The users were distributed in the room according to the uniform random distribution [43,44]. The LoS optical channel model in Equation (1) was considered, where the electrical signal power of each LED was set to 35 dBm unless otherwise stated, the default values for half the viewing angle of the LEDs and the FoV of the PDs were set to 60 • and 85 • , respectively, unless otherwise stated. The PDs faced upward towards the ceiling with an area of 1 cm 2 [33]. The number of subcarriers in each LED was set to either 16 or 32 subcarriers, and the number of users was set to 20 users, unless otherwise stated. All network parameters used in the simulation are listed in Table 2. The SA algorithm parameters were chosen as, β = 1.0005, α = 0.995, M = 50, and T 0 = 1.0. The penalty factors P 1 and P 2 , in (27), were set to 1 × 10 5 and 10, respectively.

Validation and Convergence of the Proposed Heuristic-Based Solution
In this subsection, to validate the results that the SA algorithm produces for the subcarrier allocation problem, the subcarrier allocation problem is solved using another well-known meta-heuristic, which is the TS algorithm that is also efficient in solving non-linear optimization problems [49]. For validation purposes, both the SA and the TS algorithms are simulated under the same common parameters. It is worth mentioning that there are two specific parameters to the TS algorithm: (i) The "Tabu-search list" that includes the most recently visited solutions and set to be 10, and (ii) the "Tabu-searchcandidate list" that contains the examined subset of neighborhood solutions and set to be 4. In addition, we provide the convergence curves of the SA algorithm for both the NOMA-imposed and the NOMA-not-imposed schemes. Figure 3 represents a comparison figure between the Simulated Annealing and the Tabu-search algorithms for both the NOMA-imposed and the NOMA-not-imposed schemes, while changing the total number of users in the network. From this Figure, it can be noticed that the obtained results of both algorithms are similar. This can verify that the design, the chosen parameters, and the obtained results for the Simulated Annealing algorithm are suitable and accurate.
. Figure 4a,b shows the optimization convergence curves (i.e., objective function versus iterations) of the SA algorithm, for the parameter value just mentioned, for both the NOMAimposed scheme and NOMA-not-imposed scheme, respectively. The curves indicate that the search process successfully skips through several locally optimal solutions and converges to a good-quality solution. The curve contains negative values due to the violation of constraints. The equations of the objective function Equations (27) and (28) show that the violation of constraints causes it to return negative values. It should be noted that the number of iterations along the x-axis is equivalent to the number of function evaluations because we computed the objective function once in every iteration. The number of function evaluations is a metric used in evolutionary algorithms (EAs) to denote the convergence and time relationship of EAs [41,50].   In Figure 5, the convergence analysis of the user-binding iterative greedy algorithm is provided. As discussed in Section 4.1, this algorithm takes place only when the NOMA-imposed scheme is in operation to make sure that the number of users under each LED in the network is even for forming NOMA user-pairs. In Figure 5a, box-plot type is used to show the results; therefore, for the reader's convenience, it is necessary to briefly describe the main elements of the box-plot: (i) The median of the number of iterations needed for convergence (y-axis values) is shown by the red line in the middle of the boxes, (ii) the ends of the boxes show the lower (Q1) and the upper (Q3) percentile of the y-axis values, (iii) the small horizontal lines above and beneath the boxes, called the whiskers, show the lowest and the highest value of the number of iterations needed for convergence excluding the outlier points, (iv) the outlier points are denoted here in a red plus sign and refer to the unexpected values. In Figure 5a, the number of iterations needed for convergence considering a different number of users in the network, [10,20,30,40] users, is demonstrated. One can observe from the figure that for all considered cases, the algorithm successfully prevents the LEDs from having an odd number of users with a relatively low number of iteration. Comparing the number of iterations needed for convergence in the SA algorithm, provided in Figure 4, and the number of iterations needed here, we can see that this algorithm converges with a much lower number of iterations. This is related to the reason that the search space for the sub-carrier allocation is much larger than the search space for user-binding. In Figure 5b we provide a detailed look into the number of iterations needed for convergence for the test case of N = 10. In this Figure, the convergence is verified using a histogram plot, for 100 users' location realizations, where the number of iterations needed for convergence is shown in the x-axis versus the number of users' location realizations depicted in the y-axis. For example, when the x-axis value equals zero, this indicates that around 17 (out of 100) user realizations have an even number of users under each LED from the beginning (the iterative greedy algorithm is not needed). Next, when the x-axis value equals one, this indicates that 5 (out of 100) user realizations need only one iteration in the iterative greedy algorithm to converge, and so on. This figure consolidates that the greedy algorithm always converges with a small number of iterations for all the considered 100 uniform random user realizations.  Figure 6a,b illustrates the achievable max-min user rate for the different total number of users, [10,20,30,40] users, and different total number of LEDs, [4,9], in the network. The LEDs are deployed in a square lattice topology, that is, 2 × 2 and 3 × 3. The first observation here, and in the other max-min user rate curves in the other figures (i.e., Figures 7-10), is that the NOMA-imposed scheme tends to achieve lower rates as compared to the NOMAnot-imposed scheme. The reason for this is that the NOMA scheme's imposition might pair some users with their second-best LED in order to impose the NOMA scheme fully. The advantage of implementing the NOMA-imposed scheme comes from its ability to utilize the resources more efficiently (i.e., serving the users with fewer resources) than the NOMA-not-imposed scheme that might serve one user in a cluster. The second observation here is that as the number of users in the network increases, the achievable max-min user rate decreases. This is expected, as serving more users by the same resources increases the inter-LED interference in the network and results in decreasing the achievable max-min user rates. The third observation here is that the achievable max-min user rates with 32 subcarriers available per LED are less than the achievable max-min user rates with 16 subcarriers. This is because (i) the effect of the subcarrier allocation with 16 subcarriers is more apparent on the achievable max-min user rates compared to the network with 32 subcarriers; and (ii) there is a trade-off between the high utilization of subcarriers in the LEDs and the interference in the network, in which, if we allocate a high number of subcarriers to users, then the interference in the network will be high, and subsequently, the achievable max-min user rate would decrease. The proposed subcarrier allocation with the SA algorithm takes into account this trade-off and eventually, for our setup, the simulations show that the SINR and the user rates are better with 16 subcarriers compared to 32 subcarriers. Notably, this observation is valid for small to medium values of LEDs' electrical signal power (i.e., P e < 45 dBm); however, this trend will change for large P e values, as discussed later in the discussion of Figure 9. Finally, one can observe from Figure 6b that the achievable max-min user rate decreases, again, as the number of subcarriers increases, leading to a more interference-limited scenario. Additionally, the performance of both the NOMA-imposed and the NOMA-not-imposed schemes become closer to each other compared to Figure 6a, which is due to the fact that for a larger number of active LEDs in the room, each user can receive strong channels from multiple LEDs, as LEDs become closer to each other. Subsequently, in the NOMA-imposed scheme, the binding of a user to another adjacent LED will not change its channel gain, and hence its data rate, much.  (b) For L = 9. Figure 6. The achievable max-min user rate of the NOMA-imposed scheme and the NOMA-not-imposed scheme obtained while changing the total number of users in the network. (P e = 35 dBm, φ 1/2 = 60 • , Ψ 1/2 = 85 • ). Figures 7 and 8 show the achievable max-min user rate performance while changing the semi-angle at half illumination of the LEDs and changing the FoV of the PDs, respectively, for both the NOMA-imposed scheme and the NOMA-not-imposed scheme. In both figures, as the viewing semi-angle of the LED and the FoV of the PDs increase, a slight decrease in the max-min user rate is observed. This is because, when the viewing semiangle of the LED, as in Figure 7, is increased, the coverage region of that LED increases, and hence the possibility of binding users that are far from the LED increases. On the other hand, in Figure 8, as the FoV of the PDs increases, the possibility of binding a user with a LED that is far from the user increases. Besides, with a large LED viewing angle, the signal intensity that the user receives from the LED is reduced. Likewise, with a large user FoV, the user becomes more susceptible to interference from neighbor LEDs. A similar trend (i.e., as the FoV increases the achievable rate decreases) is reported in [26] for the one-to-many case considering the sum-rate performance metric (in this paper, we consider the max-min user rate metric) in a multi-user VLC network.    In Figure 9, we present the achievable max-min user rate performance while changing the LEDs' electrical signal power (i.e., P e = [30,35,40,45,50,55] dBm) for both the NOMA-imposed scheme and the NOMA-not-imposed scheme. In this Figure, on the x-axis, we show the electrical signal power P e but the users' SINR and rate equations (i.e., Equations (3) and (4)) are written in terms of the LED optical transmit power P o . Therefore, we transform the electrical signal power values to optical transmit power values through the following relation ι = P o √ P e to substitute these values in the simulations. From Figure 9, one can observe that for small to medium P e (i.e., P e < 45 dBm) both NOMAimposed and NOMA-not-imposed schemes with 16 subcarriers have a better achievable max-min user rate performance compared to the equivalent schemes with 32 subcarriers, while the opposite is true for large P e values (i.e., P e ≥ 45 dBm). The reason is that as the power of the interfering signals gets larger, the achievable max-min user rate tends to have a larger decrease for a smaller number of the allocated subcarriers. A similar trend in RF channels, for the secrecy user rate metric, was reported in [18].

The Performance of the Proposed NOMA Schemes
Finally, in Figure 10, we illustrate the achievable max-min user rate performance while changing the room height. It is intuitive to see that when the room height increases, the channel gains of the users get weaker and subsequently, the achievable max-min user rate decreases. Besides, with a large room height, the achievable max-min user rates become comparable for both the NOMA-imposed scheme or the NOMA-not-imposed scheme and with 32 subcarriers per LED or with 16 subcarriers per LED. This is in view of the fact that the users become prone to weak signal power and weak interference power due to the comparable weaker channel gains.

Conclusions
As both VLC technology and NOMA-enabled schemes are envisioned to be among the enablers of a high data rate and low-latency future wireless networks, the achievable maxmin user rates of NOMA-enabled centralized multi-carrier VLC networks were investigated in this paper. The D-NLUPA and the bisection method were utilized for user pairing and the power allocation per pair, respectively, and a simulated annealing-based algorithm was developed for the subcarrier allocation. The obtained results, for two variants of the proposed NOMA-enabled scheme, have quantified the maximum achievable max-min user rates and shed lights on the effect of the different network parameters, such as: (i) The number of served users, (ii) the number of LEDs in the room, (iii) the semi-angle at half illumination of the LEDs, (iv) the FoV of the PDs, (v) the LED power, and (vi) the indoor space dimensions on the achievable max-min user rates.

Extensions and Future Work
The current work in this paper can be extended in the following different directions: • By exploiting illuminating LED-arrays, one can enable the utilization of multiple-input multiple-output (MIMO) in indoor VLC networks to extend the network coverage, and further increase the system capacity [51]. Investigating the max-min user rate optimization for indoor MIMO-VLC networks can be considered as a possible direction of future research. However, the performance gains may be limited due to the effect of the peak-to-average power ratio (PAPR) problem [52]. • An important practical consideration in indoor VLC networks is user mobility. The Random Way-Point model (RWP) is the most commonly used one for user mobility in indoor VLC literature [53]. In indoor multi-user centralized VLC networks, there are different solutions worth studying which can be adopted to accommodate user mobility: (i) Horizontal handover while adopting fractional frequency reuse (FFR) scheme or the use of red, green, and blue (RGB) LEDs, or allowing for a coordinated multipoint (CoMP) transmission scheme between different LEDs, (ii) vertical handover that involves RF/VLC network or WiFi/VLC network or power line communication (PLC)/VLC network, (iii) cell-zooming strategies that dynamically adjust the coverage areas of the LEDs based on user mobility profiles, and (iv) utilizing algorithms that can accommodate for user mobility by determining solutions within the coherence time of the channel [16]. • A consequence of user mobility in indoor VLC networks is LoS link blockage [53]. Thus, some novel solutions need to be adopted-for example, a multi-directional receiver or omni-directional receiver where PDs are embedded at different sides or all sides, respectively, of a smartphone. Another possible solution for the LoS link blockage can be considered by utilizing intelligent reflecting surfaces (IRSs) inside the indoor environment. Investigating the max-min user rate optimization with such solutions can be an interesting direction of future research.