Joint Access Configuration and Beamforming for Cell-Free Massive MIMO Systems with Dynamic TDD

We address the trade-off between system throughput and user equipment (UE) fairness in dynamic time division duplex (TDD) cell-free (CF)-massive multiple-input multiple-output (mMIMO) systems, developing to that end a joint access point (AP) access configuration, power allocation, and beamforming design scheme. Unlike state-of-the-art (SotA) methods, which mostly focus on sum-rate maximization or max-min worst case optimization, our design approach is to maximize the geometric mean of the UEs’ throughput performance under a constrained transmit power, which is shown to yield an excellent compromise between system throughput and user fairness. The direct reformulation of the resulting optimization problem is, however, hard to solve due not only to the non-convexity of the objective function, but also to the binary constraint induced by the AP access configuration. We thus present an efficient (i.e., polynomial time complexity) solution for the problem via a combination of fractional programming (FP) and convex–concave procedure (CCP), assisted by a negative entropy regularizer that promotes a binary solution. Numerical simulations are offered to evaluate the throughput performance and fairness index of the proposed algorithm in comparison not only to conventional TDD-based solutions, but also to recent dynamic TDD SotA designs, which illustrate the effectiveness of the proposed approach over existing methodologies both in throughput and fairness.


I. INTRODUCTION
As recently reported in [1], total wireless traffic continues to grow exponentially worldwide as a result of strong demands for high-speed and massive connectivity. In order to satisfy growing demands, the recently deployed fifth generation (5G) wireless networks rely on massive multiple-input multiple-output (mMIMO) technologies that exploit the vast spatial degrees of freedom offered by large antenna arrays [2].
Although mMIMO technologies have proven to be effective both theoretically and experimentally [3], [4], challenges such as high spatial correlation due to co-located antennas [5] and inter-cell interference at cell edges [6] still remain. One of the promising solutions to this problem is the cellfree (CF) mMIMO architecture [7]- [10], in which a large number of access points (APs), each equipped with a few antennas, are distributed over the service area and connected to one or multiple central processing units (CPUs) via wired fronthaul links. Thanks to the antenna distribution strategy of CF-mMIMO systems, harmful effects due to spatial correlation that limits effective spatial degrees of freedom can be significantly reduced, enabling the ideal performance of the original mMIMO concept without the cell-edge effect to be reached.
The CF architecture is, however, not without its own challenges, one of which is the design of beamforming schemes for distributed settings. To cite a few contributions addressing this matter, an optimized precoding strategy with power allocation for the downlink CF-mMIMO systems was proposed in [11], which aimed to maximize the minimum throughput among user equipments (UEs) based on a therederived lower bound of the achievable throughput of the system. A theoretical performance analysis of uplink CF-mMIMO systems with either zero forcing or conjugate beamformers was also offered in [12], where power allocation strategies based on a tight approximation of the capacity were proposed. More recently, a joint design of analog beam selection and precoding was proposed in [13], with the aim of reducing power consumption while enhancing performance; while the overall energy efficiency in downlink CF-mMIMO systems was considered in [14], in which a energy-efficient joint power allocation and precoding scheme was designed using optimization and game theories.
Another challenge of CF-mMIMO systems is optimizing the allocation of UE to APs. Related to the latter problem, a joint power allocation and user-AP connection design scheme was proposed in [15], aiming to maximize system throughput in time division duplex (TDD) mode. A similar attempt to optimize the AP access configuration was investigated in [16] for energy-efficient green CF-mMIMO systems, which proposed a dynamic AP ON/OFF strategy depending on traffic demand. In turn, a joint power allocation and usergrouping method under quality of service (QoS) constraints for the downlink of TDD CF-mMIMO systems was designed in [17] via generalized benders decomposition, seeking to minimize total power consumption. An AP selection method using the effective channel gain without instantaneous channel information was considered in [18].
Notwithstanding the aforementioned advances, the scalability of CF-mMIMO systems has been recognized as its main bottleneck [8]. In light of that, a new fronthaul topology enabling APs to share a serial fronthaul link with a peruser limited fronthaul bandwidth was considered in [19]. Similarly a joint power allocation and AP selection problem for uplink TDD CF-mMIMO systems was tackled in [20].
Most existing contributions in this area, including the aforementioned ones, assume a TDD-based transmission protocol with the assumption of channel duality between uplink (UL) and downlink (DL) and distinct time slots scheduled for UL and DL. As argued in [21], however, the overhead is associated with separating UL and DL communications in different time slots. Thus, It may impose inevitable delays in TDD-based systems, leading to inferior overall performance.
The concept of dynamic TDD has been proposed to overcome the latter issue in cellular mMIMO systems [22]. Dynamic TDD is characterized by the flexible allocation of UL/DL directions, based on the communication request from users in the cell. Under the CF architecture, the conventional dynamic TDD concept may be effective in reducing the possible number of standby UEs, due to the lack of a cellular structure. In light of the above, a joint UL and DL transmission scheme referred to as network-assisted full duplex (NAFD), which is similar to dynamic TDD, has been recently proposed [23] to further improve resource allocation efficiency in CF-mMIMO systems.
In the NAFD approach, APs operate in conventional half duplex (HD) mode, but each is allocated to either UL or DL transmission, such that from a network point-of-view, the system can be considered as full duplex (FD) [24]. The performance of such systems was analyzed in [25], while the AP allocation algorithm aiming to maximize the sum throughput while adopting the conventional linear minimum mean square error (LMMSE) beamforming strategy was proposed in [23]. We remark that the aforementioned works should be distinguished from FD CF-mMIMO systems, where each AP is capable of FD transmission. Contributing to the general line of work outlined above, we propose in this article a novel joint AP access configuration, power allocation, and beamforming algorithm for dynamic TDD-based CF-mMIMO systems, which aims at a reasonable optimized solution that balances system throughput and user fairness. To this end, we consider the maximization of the geometric mean of the UEs' throughput so as to maximize the total throughput utility while reducing the variance among the UE's achievable throughputs, which is motivated by recent studies such as those in [26], [27]. The resulting maximization problem involving the geometric mean objective function is highly intractable, due not only to the non-convexity of the signal-to-interference-and-noise ratio (SINR) expression, but also to the binary (discrete) feasibility set of DL/UL AP access configuration. This challenge is tackled here via a combination of a negative entropy based regularizer, fractional programming, and convex-concave procedure (CCP). The resulting algorithm is based on an iteration of convex quadratically constrained program (QCP) sub-problems. In addition, we consider both cases where each UE is equipped with a single or multiple antenna(s), tailoring the proposed algorithm for both scenarios. To the best of our knowledge, no similar joint AP access configuration, power allocation, and beamforming mechanism with the aim of addressing the performance-fairness trade-off in CF-mMIMO systems, has been proposed in the literature yet.
The remainder of the article is organized as follows. The channel and system models are described in Section II. The resource allocation and beamforming design scheme for dynamic TDD CF-mMIMO, and an initialization scheme for resource allocation method are introduced in Section III. Numerical results illustrating the efficacy of the proposed method and the associated computational analysis are offered Section IV. Finally, conclusions are given in Section V.
Throughout the text, the following notation is applied. The set of complex numbers is denoted by C. The complex conjugate is denoted by (·) * , transpose by (·) T , Hermitian transpose by (·) H , inverse of (·) −1 , and the unit matrix of N × N by I N . Also, the set difference of the set X and Y is This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. denoted by X \ Y Finally, · p denotes the p-norm.

II. SYSTEM AND CHANNEL MODEL
Consider a CF-mMIMO system comprising of L multiple APs randomly distributed over a certain coverage area, each equipped with M antennas. All APs assumed to be connected via a wired fronthaul link to a common CPU as illustrated in Fig. 1. Also, we employing a dynamic TDD scheme to jointly serve K distinct randomly distributed UEs, each of which operate in either UL or DL mode during a given channel realization. In order to focus on the achievable performance of the system, it is assumed that perfect channel state information is available at the CPU. We emphasize, however, that an extension of the approach to imperfect channel state information setups is in principle possible [28]- [30], which left for future work.
Hereafter, we shall denote the sets of all UEs by K = {1, . . . , K}, and the sets of UEs operating in DL and UL modes respectively by K dl ⊂ K and K ul = K \ K dl . Similarly, the sets of APs and of APs in DL and UL are respectively denoted by L = {1, . . . , L}, L dl ⊂ L, and L ul = L \ L dl . We emphasize that all devices are assumed to be operate in HD mode, implying that there is no overlap between K dl and K ul , and similarly between L dl and L ul .
In what follows, two distinct system setups depending on the antenna configuration at UEs will be considerednamely, single-antenna or multiple-antennas UE scenarios -with the corresponding system and channel models introduced accordingly.

A. SINGLE ANTENNA CASE
In this subsection, we first consider the case of single-antenna UEs, describing the corresponding channel model as well as the associated received signal characterization in the considered dynamic TDD CF-mMIMO system.

1) Channel Model
With k ∈ K, k ∈ K \ k, ∈ L and ∈ L \ , (i.e., k = k and = ), it is assumed that the communication channel between the k-th UE and the -th AP, as well as the interference channel between the -th AP and the -th AP are spatially correlated, whereas the channel between the k-th UE and k -th UE is uncorrelated since UEs are independently and randomly distributed over the service area, and each has a single antenna.
Given the assumption mentioned above, the three distinct channels -i.e., the communication channel h k, between the k-th UE and the -th AP, the inter-AP interference channel H , between the -th and the -th APs, and the inter-UE interference channel h k,k between the k-th and the k -th UEs -can be respectively modeled as where the small-scale fading matrixH , ∈ C M ×M are modeled as vec H , ∼ CN (0, I M 2 ).
In the latter equations, the matrices R k, ∈ C M ×M characterizes the spatial correlation of the corresponding UE-AP communication channel; while the matrices R , ∈ C M ×M and R , ∈ C M ×M are the spatial correlation matrices of inter-AP interference channels.
With regards to the channel models described above, and in particular their corresponding correlation matrices, various approaches have been proposed for characterization over the last decades [31]- [33]. To name but a few, the local scattering model based on the Gaussian angular distribution was considered in [34], whereas a low-scattered angularsparse model was studied in [35]. An azimuth spread of departure (ASD)-dependent flexible correlation model was investigated in [36], and a highly correlated scenario with a single scattered cluster can be modeled by the one proposed in [37]. In addition, uncorrelated models have also been studied in recent related works such as [38]- [41].
Among these and various other options, we adopt in this article the approach most widely employed in CF-mMIMO systems [8], [42], [43], namely the local scattering model of [44], in which the spatial correlation matrices R k, , R , and R , are characterized by the equation [44,Eq. (2.23)]: where [·] r,c is an operator that extracts the element at the r-th row and c-th column of a given matrix, g denotes a large-scale fading coefficient, d H is the antenna spacing normalized by wavelength, respectively, and the variable of integrationφ is given byφ = ϕ + δ, with ϕ denoting a certain (deterministic) VOLUME 4, 2016 angle of arrival or departure and δ ∼ N (0, σ 2 ϕ ) modeling a random fluctuation in angle spread with standard deviation σ ϕ .

2) System Model
Given the channel model described above, in this subsection we turn our attention to the system model of the corresponding CF-mMIMO network operating in dynamic TDD mode. As one can easily infer, one of the main bottlenecks of such a system is the interference due to joint UL and DL transmission, which limits system performance if not properly handled.
In order to circumvent this issue, APs are assumed to be capable of performing coordinated beamforming, whereas single-antenna UEs can only control their transmit power. Under such assumption, the UL signal vector y ul ∈ C M ×1 received at the -th AP with ∈ L ul can be written as where d ul k ∈ C represents the data signal from the k-th UL UE, w ap k , denotes the precoding vector employed by the AP toward the k -th DL UE, and n ∼ CN (0, σ 2 ul I M ) is an additive white Gaussian noise (AWGN) vector with perelement variance σ 2 ul . In light of the latter, the aggregated estimate of the data symbol of the k-th UE at the -th AP can be written as where v ap k, ∈ C M denotes the combining vector at the -th AP to detect the intended data transmitted by the k-th UL UE.
Similarly, the DL signal y dl k ∈ C received at the k-th UE with k ∈ K dl can be expressed as where d dl k ∈ C, w ap k , ∈ C M , and n k ∼ CN (0, σ 2 dl ) denotes the DL data signal toward the k-th DL UE, the precoding vector employed by the -th AP, and the AWGN affecting the k-th UE, respectively.
Since each UE possesses a single antenna, an estimate of the intended signal can be simply written as 1 Assuming that the total transmit power at each DL AP is limited to p dl,max , i.e., k∈K dl w ap k, 2 2 ≤ p dl,max , and that the power of each data symbol from UL UEs is normalized, the SINR at UL and DL can be respectively written as 1 One can consider the case where each UE is capable of multiplying a scalar quantity with the received signal y dl k . Although we omit such case here, the latter will be discussed in the subsequent section, where the more general case of UEs with multiple-antennas and beamforming capability at the receiver is addressed. where

B. MULTIPLE-ANTENNAS CASE
Next, we turn our attention to the scenario with multiple antenna UEs, generalizing the above system model accordingly.
Throughout the section, the number of antennas equipped at each UE is assumed to be N .

1) Channel Model
Since the inter-APs channel model is independent of the number of antennas at UEs, suffice it to consider the channel models associated with UEs, namely, the communication channel H k, between the k-th UE and the l-th AP, and the inter-UE interference channel H k,k between the k-th and the k -th UEs, which in the case of UEs with multiple antennas can be rewritten as where R k, ∈ C M ×M , R ,k ∈ C N ×N and R k,k ∈ C N ×N are spatial correlation matrices modeled by equation (2), while vec H k, ∼ CN (0, I M N ) and vec H k,k ∼ CN (0, I N 2 ).

2) System Model
Straightforwardly, the UL and DL received signals can be respectively written as where w ue k ∈ C N is the precoding vector employed by the k-th UE, and n k ∈ C N and n ∈ C N are the AWGN vectors at the k-th UE and the -th AP, respectively.
In turn, the estimates of the UL and DL data symbols at the k-th UE are respectively given by and where v ue k ∈ C N is the combining vector employed by the k-th UE at DL.
Consequently, the corresponding SINRs are given by where

III. PROPOSED METHOD
In this section, we propose a joint dynamic resource allocation and beamforming scheme for CF-mMIMO systems operating in the dynamic TDD mode, aimed at improving the system's total throughput, while maintaining the QoS fairness among UEs. A key ingredient of the proposed scheme is replacement of objective functions such as the total system throughput and throughput min-max, which are popular in related works [44]- [52], by the geometric mean measure, which is widely utilized in other research fields such as portfolio optimization [53]. This approach is motivated by arguments that maximizing the total throughput often leads to an unfair condition with widely varying user experiences [44], and that the maximization of minimum rate limits total system throughput [44]. In turn, it has been recently shown that maximizing the geometric mean throughput is an efficient approach to achieve the aforementioned trade-offs [26]. In view of the above, we propose an optimization-based method that maximizes the geometric mean of UL and DL data rates by jointly optimizing the beamforming weights and AP allocation.
As shall be soon clarified, the problem formulation that ensues can be tackled efficiently via Fractional programming (FP) techniques, which have been recently employed in wireless communication to deal with optimization problems involving fractional terms such as SINR and energy efficiency (EE) [54], [55]. Although techniques such as Dinkelbach's algorithm have been used in the past in association with FP methods, new frameworks that yield tractable quadratic formulations while preserving convergence guarantee with first-order optimality for sum-of-ratio problems were recently proposed [48], [54], [56], which shall be adopted here as well, due to their advantage in terms of computational efficiency. For the sake of readability without compromising completeness, an introductory description of the FP techniques used here are offered in Appendix A.

A. PROBLEM FORMULATION: GEOMETRIC MEAN MAXIMIZATION
Let us now proceed to formulating our geometric mean maximization problem aimed at the joint optimization of beamforming weights and UL/DL AP allocation, and which admits efficient solution via FP and CCP.

1) Single-Antenna UE
Consider the following weighted geometric mean maximization problem: where µ k are weighs of choice, assumed given; p dl k, is the transmit power at the -th AP towards the k-th UE, collected into the vector p dl [p dl 1,1 , · · · , p dl K,1 , · · · , p dl 1,L , · · · , p dl K,L ] T ; constraints (16b), (16c) and (16d), withη jointly enforce the selection of APs to UL and DL modes and the corresponding power limitations; constraints (16f) and (16g) enforce fairness by guaranteeing minimum numbers N ul and N dl of APs operating in UL and DL, respectively; and Γ k is the SINR corresponding to the data intended for the k-th UE, which is compactly rewritten as equation (17). As indicated by constraints (16b), (16c), and (16d), theth AP is allocated for UL transmission if and only if (iff) η = 1, or else for DL transmission iff η = 0. For the sake of notation simplicity, these binary indicator variables are collected in the vector η [η 1 , . . . , η L ] in the statement (16a) of equation (16).
The main difficulties of solving equation (16) are the facts that the objective is non-convex and that the selection variable η is binary. Fortunately this can be relaxed by the combination of FP and negative entropy regularization. As VOLUME 4, 2016 shown in Appendix B, the resulting continuous and convexified problem can be obtained as is the penalty function, λ ≥ 0 is a hyper-parameter used to adjust the strength of the penalty, and , and the SINR quantities Γ qt,1 k and Γ qt,2 k respectively defined as and as in equation (21), respectively. At this point it is worth noticing that the objective function in (18) is the difference between two concave functions, which motivates us to leverage the difference of concave (DC) programming technique. In particular, we adopt CCP to find a solution of the problem, whereby equation (18) is further modified into [57] where (·) t−1 denotes the solution obtained at the t − 1 iteration, and the first term of the objective is given by withΓ qt,2 k as given in equation (24). Finally, equation (22) can be solved by numerical convex optimization solvers such as SeDuMi and SDPT3 [58]. For Initialize: Combining vectors v ap , precorders w ap , weights λ and µ k , ∀k, incrementation parameter λ + , and maximum number of iterations t max 1 η , ∀ ← Initialize as per equation (41) 2 t ← 0 3 repeat // Update auxiliary variables : k , ∀k // Optimization of beamformer and AP allocation : 8 Update η, v ap , w ap , p dl , γ and s by solving (22) // Update hyperparameter : Then, the geometric mean maximization problem posed in equation (61), extended to multiple-antenna UEs, yields where Γ m k is as given in equation (27). The bottleneck of equation (26) is the non-convexity of equation (26a), which can be addressed by following the same approach employed in Subsection III-A1, namely, by transforming the problem using the lagrangian dual transform (LDT) and the quadratic transform (QT). For the sake of brevity, we omit detailed derivations and offer the resulting formulation, which yields subject to (16e) to (16g), (61b) to (61d), (26b) and (26c) where and Γ qt,m,2 k as given by equation (31). Equation (28) is equivalent to equation (22), extended to the case of multiple-antennas UEs. Similar to equation (22), however, the problem described by equation (28) is still not convex due to the fact that the term in equation (31) are products of different variables, which introduces variable coupling.
Finally, for fixed AP beamformers and the optimization problem to update UE beamforming counterparts is given by subject to (26b) and (26c), withΓ qt,m,1 and Γ qt,m,2 k,UE as given in equation (39). The optimization problems given in equations (32) and (36) can be solved using convex optimization solvers, leading to an alternating procedure for the joint resource allocation and beamforming design in CF-mMIMO with multipleantenna UEs, which we summarize as a pseudo code in Algorithm 2.

B. INITIALIZATION: PRE-SELECTION OF APs
Despite all the measures taken to convexize the problems whose solutions led to Algorithms 1 and 2, the fact that the underlying original problem given in equation (16) is nonconvex, and in particular the combinatorial nature of the AP allocation vector η, makes the performance of the methods described above dependent on reasonably good initializers.
We therefore introduce in this subsection a novel UL/DL pre-allocation method which can be used to initialize the AP selection vector η by taking advantage of the available channel state information (CSI) knowledge. To that end, consider for any given -th AP the k-th user such that k = argmax k∈K (E[ h k, .
Finally, let us consider a threshold ρ that characterizes the greediness to pre-allocate an AP, and initialize η by In plain words, what equation (41) implies is that: a) if a given -th AP is surrounded by an UL UE with a ρ-dominant RRSSI, then that AP is pre-allocated to UL; or else b) if conversely the AP is surrounded by a DL UE with a ρdominant RRSSI it is pre-allocated to DL; or else c) the AP is pre-allocated to UL or DL randomly.
We emphasize that this procedure is heuristic and obviously not optimal, but again, the procedure is intuitively better than a purely random pre-allocation, and only employed to initialize the vector η before the execution of Algorithm 1 or 2, from which optimized AP allocations are obtained, as shall be demonstrated in the sequel.

IV. NUMERICAL RESULTS
Before moving to a direct performance assessment of the proposed algorithms, it is useful to establish state-of-theart (SotA) benchmarks for comparison, and to analyze the complexity of the corresponding solutions, which will be pursued in the next subsections.

A. STATE-OF-THE-ART BENCHMARK SOLUTIONS
Two distinct approaches that form the basis of SotA benchmark solutions are the sum-rate maximization approach of [45]- [49], and the max-min worst-case strategy pursued in [50]- [52], both of which are briefly reviewed below for the convenience of the reader.
Since in the aforementioned articles only the case of single-antenna UEs is addressed, and since (as shall be latter shown) the proposed method already outperforms the latter in such a scenario, suffice it to consider here only benchmark approaches for the single-antenna UEs case. Under such conditions, the sum throughput maximization and minmax worst-case approaches are respectively described by the optimization problems maximize η,v ap ,w ap ,p dl k∈K subject to (16e) to (16g) and (61b), to (61d), and maximize b,η,v ap ,w ap ,p dl b (43a) subject to µ k log 2 (1 + Γ k ) ≥ b, ∀k (43b) (16e) to (16g) and (61b) to (61d), where b is the lower bound of the throughput of all UEs.
Obviously, the latter problems are not convex, due to the objective function in (42) and the constraint (43b), respectively. In existing works relying on these approaches [44]- [52], the first-order Taylor approximation is typically employed to convexize and solve the problem.
Here, however, for the sake of consistency and in order to maintain an equal footing for the subsequent comparison with our own proposed method, we instead employ the LDT and QT techniques of Lemmas 1 and 2 shown in Appendix A, respectively, without prejudice to the performance of the SotA methods [48], [54], [56].

B. IMPROVED BENCHMARK SOLUTIONS
As described in Subsection III-A, the solutions of the problems described by equations (44) and (49) still suffer from the effect of relaxing the UL/DL allocation variable η from its discrete domain η ∈ {0, 1} to its convex hull η ∈ [0, 1], which can be alleviated by the introduction of the weighted negative entropy discretizer λP(η ) into the problem.
Concisely, under such modification, the η-discretized sum throughput maximization problem of equation (42) and the max-minimum throughput problem of equation (43) subject to µ k log 2 (1 + Γ k ) ≥ b, ∀k (47b) (16e) to (16g) and (61b) to (61d), which after the convexification techniques described in Subsection III-A1, i.e., LDT, QT and CCP, respectively turn into again, with f fin k as defined in equation (23). To offer a quick glimpse on the performances of these benchmark approaches relative to one another, Figure 2 shows the cumulative distribution functions (CDFs) of the sum spectral efficiencys (SEs) achieved by both methods, with simulation parameters as specified in Section IV-D. In this and all subsequent figures, the labels "Min-rate" and "Sum-rate" refer respectively to the max-min worst-case and the sum-rate maximization methods, with the additional word "SotA" and dash-dot lines used to distinguish the methods revised in Subsection IV-A from the improved benchmark versions of this subsection.
The figure shows that indeed the SotA schemes are significantly outperformed by the improved benchmark variations, due to rounding errors in the indicator variable η . In view of these results, we will hereafter consider the improved benchmark alternatives when assessing the performance of our own proposed solution, emphasizing, however, that these benchmark schemes already incorporate contributions here introduced.

C. COMPLEXITY ANALYSIS
Before moving on to the numerical performance evaluation of the proposed Geometric Mean maximization-based method in comparison to the aforementioned benchmarks, it is of interest to analyze their computational complexities.
For starters, let us recall that the proposed method for the scenario with single-antenna UEs is fundamentally described by equation (22), while that for the case of multiple-antenna VOLUME 4, 2016 UEs is described by the pair of equations (32) and (36), which correspond to the optimization of APs and UEs, respectively. In turn, the improved sum-rate maximization and max-min rate SotA methods employed as benchmark are respectively described fundamentally by equations (48) and (49).
Next, notice that the most expensive operation associated with all these optimization problems is the computation of the corresponding QCPs ε-solution, whose canonical arithmetic complexity C can be upper-bounded by [59], [60] whereM ,Ñ , and Q m respectively denote the number of constraints in the problem, the size (i.e., vector dimension) of the real-valued multidimensional variable, and the size of the m-th constraint space; while the constant quantity digit(ε) is the order of precision of the ε-solution in terms of its distance to the optimum [60]. Leaving details to Appendix C, after transforming equation (22) dominant terms from (72c) and (72d) , where K ul and K dl are the number of UL UEs, respectively; L denotes the number of initialized APs according to equation (41); andK = log 2 (K) j=1 2 log 2 (K) −j is the number of variables (including original and auxiliary) required to reformulate the geometric mean product described in equation (22) into its QCP canonical form as in equation (74).
From equations (50) through (51c), it follows that the total complexity of Algorithm 1 can be estimated at where in the last line only the term of higher order is kept. Next, we turn our attention to Algorithm 2, which differs from Algorithm 1 in the assumption that UEs are equipped with multiple antennas, implicating on the one hand that the optimization of the APs is updated from equation (22) to equation (32), and on the other hand that UEs are capable of beamforming as described by equation (36).
Looking at the AP side first, it is evident from the comparison of equations (22) and (32) that these differ only in the objective function, which are based on equations (23) and (33), respectively. In turn, these objective functions are dependent on the quantities Γ qt,1 k andΓ qt,2 k described in equations (19) and (24) in the single-antenna UE case; and onΓ qt,m,1 k,AP andΓ qt,m,2 k,AP ) as per equations (34) and (35) in the multiple-antenna UE case, respectively. Comparing therefore equations (34) and (35) to equations (19) and (24) directly, it can be seen that the only increase in complexity in solving problem (32) as opposed to (22) is the slight increase in the number of complex multiplications corresponding to the upgrade of channel vectors h k to channel matrices H k , which is negligible compared to the costs accounted for in equation (52).
In other words, and in summary, we conclude that the computational complexity of the first part of Algorithm 2 corresponding to the optimization of APs as per equation (32) is of the same order of that of Algorithm 1, already estimated in equation (52).
In turn, the computational cost of the second part of Algorithm 2 corresponding to UE beamforming as per equation (36) can be evaluated following a similar procedure used previously, namely, by transforming equation (36) into a QCP canonical described by equation (74) in Appendix C, with the number of constraints, the real-valued multivariable dimension, and size of constraint space respectively given bỹ the total complexity of Algorithm 2 can be estimated at Finally, we consider the complexities of the total throughput maximization and min-max worst-case approaches, which are summarized fundamentally by equations (48) and (49), respectively.
To that end, first recognize that equation (48) differs from (32) only in the evaluation of a sum (as opposed to the geometric mean) in the objective function. Since this distinction has negligible impact in computational cost, it is concluded that the throughput maximization and the proposed method described in Algorithm 2 have the same complexity.
Similarly, moving the rate functions f fin k to a set of constraints, as done in the min-max worst-case approach under equation (49b), only slightly decreases the number of constraints of the QCP canonical form corresponding to equation (49), such that the min-max worst-case approach has actually a slightly lower computational complexity than the proposed summarized in Algorithm 2, which however results in a negligible overall change in the total cost.
Consequently, we will consider hereafter that the SotA benchmark methods are on par with the proposed method in terms of computational complexity.

D. COMPUTER SIMULATIONS 1) Simulation Set-up
Having established the complexity analysis of the proposed geometric-mean-based scheme for joint AP access configuration and beamforming, we proceed to numerically evaluate the performance of the proposed method under various conditions, comparing it with the SotA benchmarks.
In what follows it is assumed that APs are distributed in a regular grid formation within the effective service area, while UEs are randomly and uniformly distributed over the area. The remaining simulation parameters are summarized in Table 1, with the quantity d used in the modeling of largescale fading denoting the distance between the transmitter and receiver, and the quantity z representing shadowing with a standard deviation of 10 [dB], as considered in related CF-mMIMO literature such as [8], [44]. Also following [8], [44], it is assumed that APs have a height of 10 m above the ground surface.
In order to focus on the assessment on the improvement provided by the proposed method itself, perfect CSI knowledge is assumed at transmitters and receivers. This assumption have been widely used in other CF mMIMO literatures [61]- [64]. The convex optimization problems (22), (32), and (36) in Algorithms 1 and 2 are solved via CVX [58], using the SDPT3 backend and with the maximum number of iterations set to t max = 30.
In each realization of simulations performed, the weight applied to the negative entropy function is initialized to λ = 0 and, after the fifth iteration, increased at each iteration by the increment λ + , which is numerically optimized so as to maximize the system throughput of each tested method, including the improved benchmark SotA approaches described in Subsection IV-B. Finally, in order to avoid numerical instability, the solution obtained at the previous iteration is utilized as the final output of the algorithm if the solver outputs an invalid solution in the middle of the loop.

2) Single-antenna UEs Case
Let us start by evaluating the performance of the proposed and SotA methods in the scenario where each UE has a single antenna. Figure 3 shows the CDFs of the sum SEs for the different approaches compared, and different ratios of UEs operating in UL mode, as indicated by the quantity q.
The performance of the "MMSE" method in conventional TDD mode is also included as a reference, with UL and DL minimum mean square error (MMSE) beamformers respectively obtained from v ap (55b) where DL MMSE beamformer [65]- [67] is allocated by same power for all UEs as p dl k = p dl,max /(qK). In all figures, the legend "Geo-Mean" refers to the performance of Algorithm 1, while "Sum-Rate" and "Min-Rate" refer to improved sum-rate maximization and max-min worst case benchmark approaches via equations (48) and (49), respectively. The legend "DTDD" corresponds to a system that operates in dynamic TDD regardless of the AP access configuration, but with beamforming design employed, whereas "TDD" corresponds to a TDD-based system with time resource partition relative to the access request ratio q. The additional token "w/o Init" added to some of the legends, refer to systems in which the initialized pre-allocation of APs according to equation (41) is not performed. The curves corresponding to such alternatives are plotted with dash lines for further visibility. Finally, in all methods compared, the beamforming quantities are initialized by the maximum ratio (conjugate) beamforming method [8].
From Figure 3 -which covers a predominantly DL scenario (q = 0.25), a predominantly UL scenario (q = 0.75), and a balanced UL/DL scenario (q = 0.5) -it is found that dynamic TDD yields significantly higher SEs compared TDD mode, both with the proposed and SotA methods. This indicates that thanks to joint AP allocation and beamforming optimization, dynamic TDD takes better advantage of time resources to serve users, and of spatial resources to control harmful interference, ultimately resulting in higher rates. It is also observed that among all methods compared, the min-max worst-case approach is the most sensitive to initialization, with large performance gains resulting from the preallocation of Subsection III-B observed under that particular method. The min-max worst-case scheme is also found to be the most sensitive to the UL and DL demands, with large differences in overall SE observed between predominantly DL, balanced, and predominantly UL scenarios.
In turn, it is found that both the sum-rate maximization and the proposed geometric mean maximization approaches are most reliable and best-performing methods, with the latter slightly outperforming the former especially under predominantly DL conditions. At this point we might once again emphasize, however, that the sum-rate maximization method benchmark whose results are shown in Figure 3 already include the improvements described in Subsection IV-B, such that the gains of the proposed method against SotA as-is, e.g.
We also remark that both the min-max worst-case and sum-rate maximization methods impose a trade-off between performance and QoS fairness, sacrificing either of the two important communication requirements. In contrast, as shall be shown latter, the proposed algorithm is designed to tackle this issue avoiding the sacrifice of fairness while improving upon the SE performance of SotA alternatives.

3) Convergence
Before evaluating the fairness of the methods compared, let us briefly study the average convergence behavior of the three distinct approaches in terms of the number of APs allocated to UL, under the same balanced conditions of Figure 3(b).
To this end, we plot in Figure 4 the values of η obtained with the improved benchmark SotA and proposed methods, as a function of the number of iterations. The results indicate that the AP access configuration reached by the proposed scheme differs significantly from those of the improved SotA methods. It is also noticeable, in fact, that the proposed approach is the only one that arrives at a configuration in which the same number of APs are allocated to UL and DL, as intuitively expected for the balanced UL/DL case considered. Given that fundamentally only the objective function of the proposed method differs from those of the improved SotA alternatives, and since the same negative entropy based discretizer is utilized in all methods compared, this finding confirms the efficacy of the geometric mean rate, as opposed to the other more traditional figures of merit, as the objective of joint allocation and beamforming optimization techniques for CF-mMIMO systems.
Finally, we compare in Figure 5 the convergence behavior of the proposed and improved benchmark SotA algorithms in terms of the achieved average SE computed with the continuously-relaxed η at each iteration. Notice that the SE performance shown in Figure 5 is not the actual SE with discretized η ∈ {0, 1}. At the last iteration (t = 30), however, we fix η 's to either 0 or 1, which results in slight decrement in performance due to rounding at t = 30. It can be observed from the figure that the SE continuously increases as the number of iterations grows, except for the point t = 5 when the penalty parameter starts to add up.

4) Fairness
Next, we study the fairness among UEs achieved by the proposed and benchmarking joint AP allocation and beamforming schemes considered. To that end, plots of the CDFs of the minimum SE among UEs achieved with each system are shown in Figure 6. The results show that among the dynamic TDD techniques compared, the sum-rate maximization method is the one that exhibits the poorest minimum SE performance, regardless of the UL/DL ratio q, which is as expected since this approach tends to allocate spatial degreesof-freedom to the UEs with higher SINRs, in detriment of UEs with lower SINRs, leading to a maximization of total throughput at the expense of unfair conditions for UEs. Interestingly, the proposed geometric mean approach is found to be the best among all those compared in terms of the minimum SE performance, outperforming even the max-min worst case scheme, despite the fact that the latter is designed precisely to maximize the rate of the weakest user. Although this result can seem counter-intuitive at first, it can be justified as follows.
Firstly, in an average sense, the max-min worst case method is significantly more sensitive to the initial state of the AP allocation variable vector η than its counterparts, as previously observed in Figure 3 and now further evidenced by the results of Figure 6. Secondly, in a per-realization sense, the max-min worst case approach is also found to be, among the methods compared, the one that exhibits the largest drop at 5th iteration, in which the incremental penalty λ + is added so as to promote a binary selection for the AP access configuration as seen in Figure 4. These results together indicate that the max-min worst case is the method that proves least capable of handling the binary nature of the AP allocation variable.
In contrast, the proposed geometric mean approach is effective in addressing the non-convexity and the binary constraint on the AP access configuration, resulting in fairer solutions than the max-min worst case method.
These results notwithstanding, in order to further evaluate the fairness of the three techniques, we compare their Jain's fairness indices [68], defined by Jain's Fairness This index lies between 1/K to 1. Thus, when the index be- Non-surprisingly, here the max-min worst case algorithm is shown to deliver the most fair results among the dynamic TDD schemes compared, which however comes at the expense significantly lower rates, as shown previously. This is also why the TDD mode results are the ones with the highest Jaines fairness indices, since under such method there is no inter-AP/UE interference, although the throughput itself is inferior. All in all, it is concluded from Figures 3 to 7 that the proposed geometric mean method yields the best excellent compromise of performance and fairness among those compared.
To elaborate, the proposed geometric mean method outperforms the greedy sum-rate maximization scheme in achieving the largest total rate, as shown in Figure 3, while exhibiting Jain's fairness indices not far from the latter, as seen in Figure 7; outperforms the fairness-oriented max-min worst case method in ensuring a higher rate to poorest UE, as seen in Figure 6; and achieves the best balance in the allocation of APs to UL/DL, as seen in Figure 4.

5) Multi-antenna UEs Case
Having demonstrated the effectiveness of the proposed geometric mean approach compared to SotA benchmarks, we now turn our attention to the gains obtained by increasing the number of antennas at UEs, under the correlated channel model of [44], described by equation (2), with the angular standard deviation is σ ϕ = 20 [deg] and the antenna spacing is d H = 1/2.
To that end, we compare in Figure 8 the SE performance of Algorithm 2 against that of the sum-rate maximization approach, for various cases with different number N of antennas at each UE. For N = 1, the performance of Algorithm 1 is also shown as a further reference. In addition, we clarify that since the maximum ratio transmission can not be directly utilized as a beamforming vector in multipleantennas scenarios, a Gaussian initialization is adopted to initialize the beamforming vectors in this figure.
It can be observed from the figure that the SE performances of both the proposed geometric-mean and sum-rate maximization approaches increase with N , with a significant advantage of the proposed method over the SotA, both in terms of total and minimum SE. Here, it is worth emphasizing that due to the geographical nature of the local scattering model [44], in which the spatial distribution of UEs and APs directly determine the parameters d H and ϕ, strongly influencing the correlation matrices R k, , the multiple-input multiple-output (MIMO) channels to which Algorithm 2 is subjected to tend to be highly correlated. This, in turn, indicate that a large portion of the gains obtained by Algorithm 2 is due not to the added channel diversity, but rather to the effective exploitation of the added degrees of freedom, by the proposed joint AP allocation and beamforming scheme, to aligning harmful interference while increasing the intended signal power.   It is also worth-noting that for N = 1, the performance of Algorithm 2 is essentially identical to that of Algorithm 1, as expected, with only a slight difference in terms of sum SE observed, which is attributed to the scalar weight "precoding" included in Algorithm 2, but not in Algorithm 1.
Finally, it is also noticeable that the performance of sumrate maximization approach degrades greatly compared to Algorithm 2, especially when N is small. Given that in such scenarios (N → 1) UEs are more limited, this indicates that Algorithm 2 is indeed capable of finding a better solution to the allocation of APs than the sum-rate maximization approach.

V. CONCLUSION
We studied a dynamic-TDD aided CF-mMIMO system in which UL and DL UEs are simultaneously served by the system over shared wireless resources, proposing a novel joint AP access configuration and beamforming scheme to optimize both the system throughput and user fairness of such a system, both for the case of UEs with single and multiple antennas. The proposed method is based on geometric-mean figure of metric, incorporating also a novel mechanism to enforce the discreteness of the AP allocation variables η via the introduction of a regularized negative entropy function in the objective. In order to relax the non-convexity of the newly formulated problem, FPs techniques were employed; and a simple heuristic to initialize η is also introduced. It is shown via software simulations that the proposed method is not only capable of balancing rate performance and fairness, but in fact outperforms the conventional sum-rate maximization and max-min worst case methods. Having verified the effectiveness of the proposed approach in dynamic-TDD aided CF-mMIMO systems, open problems such as imperfect CSI, AP cooperation under limited fronthaul links, and hyperparameter tuning are left as some of the possible directions for future work.

APPENDIX A. INTRODUCTION OF FP
In this section, we introduce FP techniques that enable a convex approximation of a function involving fractions and logarithmic terms. With these remarks made, and for the sake of clarity of exposition, it will prove useful to compactly summarize the following approximation techniques, respectively referred to as the LDT and the QT and proposed in [48], [56].  where α m (x) and B m (x) are given in Lemma 1.
Applying QT to the above sum-of-ratio maximization problem, we obtain where S = [s 1 , . . . , s M ] denotes a newly introduced auxiliary variable that possesses an optimal solution given by s m = B −1 m (x) α m (x). Together, Lemmas 1 and 2 yield an iterative procedure to solve sum-of-logarithms problems on the auxiliary variables s and γ, whose convergence was studied in [48], [56].

APPENDIX B. CONVEXIFICATION VIA FRACTIONAL PROGRAMMING AND NEGATIVE ENTROPY
In this section, we discuss in detail what was omitted in Section III-A1, namely, the conversion of binary variable η to continuous numbers and the convexification of the geometric mean objective function.
Despite this reformulation and the removal of the combinatorial issue on η, equation (61) is still intractable because of the non-convexity of the objective function on the remaining variables. To resolve this challenge, Lemma 1 can be applied to equation (61) yielding maximize η,v ap ,w ap ,p dl ,γ K k=1 f k v ap , w ap , p dl , γ 1 K + L =1 λP(η ) (62) subject to (16e) to (16g) and (61b) to (61d).
Equation (62) can be seen as a weighted SINR maximization problem, wherein the binary AP selection is aided by the negative entropy penalty function. This problem is, however, still not straightforwardly solvable due to the convex-overconvex formulation of the SINR expression, which however can be easily convexified via the QT. Thus, it applies Lemma 2, the SINR expression in equation (62) can be reformulated to convex function as equation (18).

APPENDIX C. CANONICAL FORM OF QCP FORMULATION
In this section, we describe how problem (22) is put into the QCP canonical from, so as to enable the calculation of its computational complexity. To that end, first consider the canonical form of a real-valued conic QCP, which is expressed as [60] A m x + b m 2 ≤ c T m x + d m , ∀m ∈ {1, . . . ,M }, (65c) x ∈ RÑ , b m ∈ R Qm .
Next, consider the following minimization problem equivalent to the maximization problem of equation (22) (67) must be expressed in terms of real variables, with a distinction for the UEs operating in UL and DL, respectively.