Energy-Efficient Power Control in Cell-Free and User-Centric Massive MIMO at Millimeter Wave
Mario Alonzo Stefano Buzzi Alessio Zappone Ciro D'Elia


IEEEexample:BSTcontrol

This paper was submitted for publication on the IEEE Transactions on Green Communications and Networking on April 24, 2018 and was assigned reference number TGCN-TPS-18-0048.
It was finally accepted for publication on March 24, 2019.  




2019 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted  component of this work in other works.”






Energy-Efficient Power Control in Cell-Free and User-Centric Massive MIMO at Millimeter Wave
    Mario Alonzo, Stefano Buzzi, Senior Member, IEEE,  Alessio Zappone Senior Member,  IEEE, and Ciro D'Elia
This paper was partly presented at the 28th Annual IEEE International Symposium on Personal, Indoor and Mobile Radio Communications, Montreal (Canada), October 2017 and at the 2018 IEEE 5G World Forum, Santa Clara, CA, July 2018. The work of M. Alonzo, S. Buzzi, and C. D'Elia  has been supported by the Italian Ministry of Education and Research, under the program "Dipartimenti di Eccellenza 2018-2022". The work of A. Zappone was supported
by H2020 MSCA IF BESMART, Grant 749336.
M. Alonzo, S. Buzzi and C. D'Elia are with the Department of Electrical and Information Engineering, University of Cassino and Lazio Meridionale, I-03043 Cassino, Italy (mario.alonzo@unicas.it, buzzi@unicas.it, delia@unicas.it). S. Buzzi and C. D'Elia are also affiliated with CNIT, Parma, Italy. Alessio Zappone is with the LANEAS group of the L2S, CentraleSupelec, CNRS, UnivParisSud, Université Paris-Saclay, 91192 Gif-sur-Yvette, France. (alessio.zappone@l2s.centralesupelec.fr).
    
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

empty


In a cell-free massive MIMO architecture a very large number of distributed access points simultaneously and jointly serves a much smaller number of mobile stations; a variant of the cell-free technique is the user-centric approach, wherein each access point just serves a reduced set of mobile stations. This paper introduces and analyzes the cell-free and user-centric architectures at millimeter wave frequencies, considering a training-based channel estimation phase, and the downlink and uplink data transmission phases.  First of all, a multiuser clustered millimeter wave channel model is introduced in order to account for the correlation among the channels of nearby users; second, an uplink multiuser channel estimation scheme is described along with low-complexity hybrid analog/digital beamforming architectures. Third, the non-convex problem of power allocation for downlink global energy efficiency maximization is addressed.
Interestingly, in the proposed schemes no channel estimation is needed at the mobile stations, and the beamforming schemes used at the mobile stations are channel-independent and have a very simple structure. Numerical results show the benefits granted by the power control procedure, that the considered architectures are effective, and permit assessing the loss incurred by the use of the hybrid beamformers and by the channel estimation errors.  




Cell-free massive MIMO, user-centric approach, millimeter wave, energy efficiency, sum-rate, power control, wireless networks. 






§ INTRODUCTION

Future 5G wireless systems will heavily rely on the use of large-scale antenna arrays, a.k.a. massive MIMO, and of carrier frequencies above 10 GHz, the so called mmWave frequencies <cit.>. Indeed, on one hand, the use of massive MIMO permits serving several users on the same time-frequency resource, while, on the other one, mmWave carrier frequencies will enable the use of much larger bandwidths, at least on short-distances (up to about one hundred meters). The combined use of massive MIMO systems along with mmWave frequencies is indeed one of the key technological enablers of the envisioned wireless Gbit/s experience <cit.>. For conventional sub-6 GHz frequencies, a new communications architecture, named CF massive MIMO, has been recently introduced in <cit.>, in order to alleviate the cell-edge problem and thus increase the system performance of unlucky users that happen to be located very far from their serving AP. In the CF architecture, instead of having few base stations with massive antenna arrays, a very large number of simple APs randomly and densely deployed serve a much smaller number of MSs.
This approach has some similarities with the CoMP <cit.> and network MIMO <cit.>; these technologies exploit coordinate beamforming, scheduling, and  joint transmission using multiple distributed antennas and/or base stations, with the aim of mitigating the interference and achieving diversity gains. There are however some  key differences between these strategies and the CF massive MIMO approach discussed in this paper. In particular, we have that: (a) the CF massive MIMO strategy fully leverages the advantages and unique features of massive MIMO, such as the channel hardening effect and the use of the TDD protocol;and (b) the CF massive MIMO also makes a more limited use of the backhaul link since the channel coefficients are locally estimated at the APs and are mutually shared, and, moreover, beamformers are computed at the APs using locally available information.
The CF massive MIMO architecture can be thus thought as the scalable way of implementing distributed network MIMO deployments.  
 In the CF architecture described in <cit.>, single-antenna APs and MSs are considered, all the APs serve all the MSs, all the APs are connected through a backhaul link to a CPU, but every AP performs channel estimation locally, and channel estimates are not sent to the CPU, but are locally exploited. In particular, for the downlink communication phase, the CPU sends to the APs the data symbols to be sent to all the MSs, while for the uplink communication phase, the APs use the backhaul to send the sufficient statistics for all the MSs to the CPU, which combines them and performs uplink data decoding.  
In <cit.> a UC variant of the CF approach is introduced, wherein each APs, instead of serving all the MSs in the considered area, just serves the ones that he receives best; the results of <cit.> show that the UC approach provides savings on the required backhaul capacity and, also, provides better data-rates to the vast majority of the users. 

The paper <cit.> analyzes a CF massive MIMO system from the point of view of its energy consumption, taking into account the power consumed by the backhaul links, the number of active APs and the number of antennas at each AP. The paper also embraces the UC philosophy by considering AP-MS selection schemes aimed at maximizing the system energy efficiency. While previous papers deal with the case in which the APs and the MSs are both equipped with one antenna (with the exception of <cit.> which considers multiple antennas at the APs), in <cit.> the CF and UC  architectures are generalized to the case in which both the APs and the MSs are equipped with multiple antennas. This generalization is not trivial and indeed the system proposed in <cit.>, 
despite the use of multiple antennas at the MSs and the use of a multiplexing order larger than one, develops a beamforming scheme at the APs that does not require channel estimation at the MSs, which  adopt a channel-independent beamformer.
While <cit.> considers a simple uniform power allocation, in <cit.> power control procedures aimed at rate maximization are proposed.
The paper <cit.> investigates the use of pilot signals on the downlink, so as to enable channel estimation at the MSs; even though such an approach provides some performance improvement, it contradicts the TDD-based "Massive MIMO philosophy" wherein no channel estimation is required at the MSs. A compute-and-forward approach to CF massive MIMO was then proposed in <cit.> in order to reduce the load on the backhaul links. 

All of the above cited papers consider the case in which a conventional sub-6 GHz carrier frequency is used. On the other hand, as already stated, mmWave frequencies will play a big role in future wireless cellular systems due to availability of wide amounts of unused bandwidths <cit.>. The design of a wireless cellular system operating at mmWave poses different and new challenges with respect to conventional sub-6 GHz frequencies, due to different propagation mechanisms <cit.>, to the need of using large antenna arrays at both side of the links to counteract the increased 
path-loss <cit.>, and to the presence of hardware constraints that prevent the realization of FD precoding and postcoding beamforming structures <cit.>. Despite these tough  challenges, the attractiveness of the mmWave frequencies for cellular communications has led to intense research efforts in the last few years. 
Due to the difficulty of realizing FD beamforming architectures at mmWave with antenna arrays of large size, the CF massive MIMO architecture at these high frequencies appears particularly attractive, mainly for two reasons; first of all, it substitutes large co-located antennas with several APs 
with antenna arrays of smaller dimension, and, thus, with less hardware complexity; second, due to the limited range of mmWave communications, the distributed AP dense deployment alleviates the cell-edge problem and creates path-diversity, which is precious since signal blockages may frequently happen at these frequencies. 

Energy efficiency is another key topic that should be considered when designing modern wireless communication systems. Indeed, increasing the  bit-per-Joule energy efficiency is regarded as a key requirement of future 5G networks <cit.>. A recent survey on the most promising energy-efficient techniques for 5G has recently  appeared in <cit.>; in that paper,  radio resource allocation is identified, among other techniques, as a relevant approach to improve the energy efficiency of future wireless networks. Recent contributions on energy-efficient radio resource allocation for 5G are <cit.>.

Motivated by the above discussion, this paper considers details a CF and UC massive MIMO wireless system operating at mmWave and focuses on resource allocation strategies maximizing the global energy efficiency. To the best of authors'  knowledge, this is the first paper that considers the CF and UC massive-MIMO deployments at mmWave; preliminary results on this issue have been reported by the authors in the conference papers <cit.> and <cit.>,  wherein it has also been shown that the UC approach generally
outperforms the CF approach.

The contribution of the paper can be summarized as follows. First of all, we introduce a multiuser mmWave channel model that permits taking into account channel correlation for close users. Building upon the clustered channel model <cit.> widely used at mmWave frequencies, we extend this model to take into account the fact that if several APs and MSs are in the same area, their channels must be built using the same set of scatterers; adopting this model, users that are very close will receive beams with very close direction of arrival and so channel correlation for nearby users is intrinsically taken into account. Then, we study the UC and CF approaches at mmWave frequencies; we assume that both the APs and MSs are equipped with multiple antennas, use HY analog-digital partial ZF beamforming at the APs, while a very simple 0-1 beamforming architecture, independent of the channel estimate, is used at the MSs. 
We present simulation results for scenarios that can be representative of a lightly-loaded system, with M=80 APs and K=6 MSs transmitting on the same time-frequency slot,  and of a highly-loaded system, with M=80 APs and K=16 MSs.
Results show that over a bandwidth of 200 MHz and with a maximum transmit power of 0.1 W at each AP, taking into account channel estimation errors and using low-complexity beamforming structures, the downlink average rate per user is 1.5 Gbit/s in the lightly-loaded scenario and 400 Mbit/s in the heavily-loaded situation. Likewise,  in the lightly-loaded scenario, the average uplink rate-per-user is about 1 Gbit/s. 

This paper is organized as follows. Next section is devoted to the discussion of the system and to the description of the processing (it involves uplink training, downlink data transmission and uplink data transmission); Section <ref> concerns global energy efficiency maximization; Section <ref> concerns the sum-rate optimization under (theoretical) perfect ZF beamforming; Section <ref> contains details about the used multiuser channel model and the discussion of the numerical results. Finally, concluding remarks are given in Section <ref>.

Notation. We use lower case boldface characters to denote column
vectors and capital boldface characters to denote matrices.
We denote N × N identity matrix with I_N; (·)^T, (·)^H and (·)^* denote the transpose, conjugate transpose and complex conjugate of a matrix. Then we denote with tr(A), ∥A∥, |A| the trace, the norm and the determinant of a matrix A. x∼𝒞𝒩(μ,σ^2) indicates a Complex Gaussian random variable with mean μ and variance σ^2. ⊗ indicate the Kronecker product between two matrices. We denote with card(·) a set's cardinality. A list of used acronyms is reported in Table I.






§ SYSTEM MODEL AND COMMUNICATION PROTOCOL


We consider an area where K MSs and  M APs are randomly located.
The APs are connected by means of a backhaul network to a CPU wherein data-decoding is performed (see fig. <ref>). Communications take place on the same frequency band; downlink and uplink are separated through TDD[In TDD, using calibrated hardware, the uplink channel is the reciprocal of downlink channel.]. The communication protocol is made of three different phases: uplink training, downlink data
transmission and uplink data transmission. During the uplink training phase, the MSs send pilot sequences to the APs and each AP estimates the channels; during the second phase the APs use the channel estimates to perform pre-coding and transmit the data symbols; finally, in the third phase the MSs send uplink data symbols to the APs. 
While in the CF approach all the APs simultaneously serve all the MSs (a fully-cooperative scenario), in the UC approach each  AP serves a pre-determined number of MSs, say N, and in particular the ones that it receives best.

We assume that each AP (MS) is equipped with an antenna array with N_AP (N_MS) elements. 
The (N_AP× N_MS)-dimensional matrix H_k,m denotes the channel matrix between the k-th user and the m-th AP. Details about the channel model will be reported in Section <ref>

Additionally, we assume that each MS employs a very simple 0-1 beamforming structure; in particular, denoting by P the multiplexing order, namely the number of parallel streams sent to a given receiver, 
 the (N_MS× P)-dimensional beamformer used at the k-th MS receiver is denoted by L_k and is defined as
 L_k=I_P⊗1_N_MS/P, denoting with 1_N_MS/P an all-1 vector of length N_MS/P. Otherwise stated, we assume that the MS receive antennas are divided in P disjoint groups of N_MS/P elements, and the received data collected at the antennas of each group are simply summed together. It is APs' task, based on the uplink channel estimates and exploiting the TDD channel reciprocity, to ensure that the summed samples are, at least approximately, aligned in phase. 
Similarly when considering uplink transmission, the antennas in each group send the same signal with the same phase. 
 We describe now the three phases of the communication protocol. 




 §.§ Uplink training

During the uplink training the MSs transmit pilot sequences in order to enable channel estimation at the APs. Let 
τ_c be the length  of the channel coherence time and τ_p be the length of uplink training phase, both in discrete time samples.  Of course we must have τ_p < τ_c. We define by Φ_k ∈𝒞^P×τ_p the matrix containing on its rows the pilot sequences sent by the k-th MS. We assume that 
Φ_kΦ_k^H =I_P,  i.e. the rows of Φ_k are orthogonal, but no orthogonality is required for the pilot sequences assigned to other MSs[Of course, when KP≤τ_p it would be possible to assign to all the MSs mutually orthogonal pilot sequences. In general, the pilot sequences  contain symbols from the same constellation as the data symbols, e.g. QAM symbols or similar. In this paper, however, for the sake of simplicity, we assume that the pilot sequences are binary random sequences, and we just require that each matrix Φ_k has orthogonal rows.]. Obviously, using orthogonal pilots tout court would lead to a system immune to pilot contamination, but this would put a limit on the maximum value of the product KP that could be accommodated in the channel coherence time. 
  The received signal at the m-th AP in the τ_p signaling intervals devoted to uplink training can be cast in the following  N_AP×τ_p-dimensional matrix Y_m: 
  
    Y_m=∑_k=1^K√(p_k)S_k,mΦ_k+W_m,

  where S_k,m=H_k,mL_k, W_m is the matrix of thermal noise samples, whose entries are assumed to be i.i.d. 𝒞𝒩(0,σ_w^2) RVs. 
  In the following we briefly outline the structure of the LMMSE channel estimator.
  Defining y_m=vec(Y_m), w_m=vec(W_m), s_k,m=vec(S_k,m), we obtain the vectorized model: 
  
    y_m=∑_k=1^K√(p_k)A_ks_k,m+w_m

  with A_k=Φ_k^T⊗I. Since we need to estimate a vector, we process y_m by a matrix V_k,m^H, i.e. ŝ_k,m=V_k,m^Hy_m.  
  Then, the MSE is:
  
    𝔼[∥V_k,m^Hy_m - s_k,m∥^2]=tr(V_k,m^H𝔼[y_my_m^H]V_k,m)+
       +𝔼[∥s_k,m∥ ^2]-𝔼[2{tr(s_k,m^HV_k,m^Hy_m)}]=
       =tr(V_k,m^H(∑_l=1^Kp_kA_lA_l^H+σ^2I)V_k,m)+
       +𝔼[∥s_k,m∥ ^2]-√(p_k)tr(V_k,m^HA_k+V_k,m^TA_k^*)

Recalling that   (<cit.>)   
 ∇_Z^*tr(Z^TM)=0, ∇_Z^*tr(Z^HM)=M, and ∇_Z^*tr(Z^HMZ)=MZ, setting 
the gradient of the MSE with respect to the complex matrix V_k,m^* equal to zero and solving for   V_k,m, we find the LMMSE estimator:
  
    V_k,m^LMMSE=√(p_k)(∑_l=1^Kp_lA_lA_l^H+σ^2I)^-1A_k  .








 §.§ Downlink data transmission

After the first phase, the generic m-th AP has an estimate of the quantities S_k,m, for all k=1, …, K. 
In order to transmit data on the downlink, a ZF precoder is considered.
The precoding matrix 𝐐_k,m is designed as follows.
 First of all, the matrix 
𝒢=[S_1,1…S_K,M] is built;  then, we have

    Q_k,m=(𝒢𝒢^H)^-1S_k,m.

Finally, each precoding matrix is normalized as follows: 

    Q_k,m=Q_k,m√((Q_k,mQ_k,m^ H))  ,

∀ k=1,…,K,  ∀ m=1,…,M.

The previously described beamforming matrix  requires a FD implementation, which presumes the use of a number of RF chains equal to the number of transmit antennas. It is well-known that at mmWave frequencies hardware complexity constraints usually prevent the use of FD architectures, and thus HY beamforming structures have been proposed, where a number of RF chains N_AP^RF< N_AP is used.  
In this paper we exploit the Block Coordinate Descent algorithm <cit.> in order to decompose our beamformer in the cascade of a FD (baseband) one, represented by a (N_AP^RF× P)-dimensional matrix and of an analog one, represented by a (N_AP× N_AP^RF)-dimensional matrix whose entries have all constant norm. While the  baseband beamformer is MS-dependent, the analog beamformer  at the AP is unique and is used for transmitting to all the users. 
In particular, at the generic m-th AP, the following matrix, of dimension N_AP× KP, is 
formed[In the UC approach, to be detailed in the following, the matrix 𝐐_m^opt is
formed using the beamformers relative to the MS actually served by the m-th AP. ]

    𝐐_m^opt=[𝐐_1,m, …, 𝐐_K,m]   ,

and used as an input to the algorithm <ref> in order to obtain the N_AP× N_AP^RF matrix 𝐐_m^RF and the N_AP^RF× KP matrix 

    𝐐_m^BB=[𝐐_1,m^BB, …, 𝐐_K,m^BB]  ,
 
containing the MS-specific baseband beamformers to be used at the m-th AP. At the generic m-th AP, we will thus have as many digital beamformers as the MSs to transmit to, and only one analog beamformer, that will be used to transmit jointly to all the users. Although we are here considering, for the sake of simplicity, only one single frequency, it should be noted that with a multicarrier modulation a single analog beamformer must be used for all the subcarriers, or, if complexity permits, for each properly defined subset of contiguous subcarriers <cit.>. 




  §.§.§ The CF approach
 In this case, all the APs serve all the MSs, so the transmitted signal from the m-th AP in the n-th sample interval is: 

    s_m^CF(n)=∑_k=1^K√(η_m,k)Q_k,mx_k^DL(n),
 
where x_k^DL(n) is the data symbol intended for the k-th MS, and  η_m,k is a scalar coefficient taking into account the power used by the m-th AP to transmit towards the k.th MS.  
The k-th MS receives the following (N_MS× 1)-dimensional vector: 

    y_k^CF(n)   = ∑_m=1^MH_k,m^Hs_m^CF(n) + z_k(n)=
       =∑_m=1^M√(η_m,k)H_k,m^HQ_k,mx_k^DL(n)+
       +∑_l=1, l≠k^K∑_m=1^M√(η_m,l)H_k,m^HQ_l,mx_l^DL(n) + z_k(n),

where z_k(n) is the additive thermal noise distributed as 𝒞𝒩(0,σ_z^2).
A  soft estimate of the k-th MS data symbol is thus formed as 

    x_k^DL,CF(n)=L_k^Hy_k^CF(n).




  §.§.§ The UC approach
 In this case  the APs are assumed to serve a pre-determined, fixed number of MSs, say N; in particular, we assume that the generic m-th AP serves the N MSs whose channels have the largest Frobenious norms.  We denote by 𝒦(m) the set of MSs served by the m-th AP. Given the sets 𝒦(m), for all m = 1,…,M, we can define the set ℳ(k) of the APs that communicate with the k-th user: 

    ℳ(k)={m:k∈𝒦(m)}.

In this case the transmitted signal from the m-th AP is written as: 

    s_m^UC(n)=∑_k∈𝒦(m)√(η_m,k)Q_k,mx_k^DL(n).

The received signal at the k-th MS is expressed now as: 

    y_k^UC(n)   = ∑_m=1^MH_k,m^Hs_m^UC(n) + z_k(n)=
       =∑_m∈ℳ(k)√(η_m,k)H_k,m^HQ_k,mx_k^DL(n) +
       +∑_l=1, l≠k^K∑_m∈ℳ(l)√(η_m,l)H_k,m^HQ_l,mx_l^DL(n) + z_k(n),

As before the N_MS-dimensional vector z_k(n) represents the thermal noise at the k-th MS, and it is modeled as i.i.d. 𝒞𝒩(0,σ_z^2). Then it is possible to obtain a soft estimate of the data symbol x_k^DL(n) at k-th MS as:

    x̂_k^DL,UC(n)=L_k^Hy_k^UC(n).




Given the above equations, and assuming the use of Gaussian distributed long codewords, it is now possible to write down the achievable  rate for the k-th user in the UC approach as follows: 

    ℛ_k=Blog_2|I+R_k^-1A_k,kA_k,k^H|  ,

where

    A_k,l=∑_m∈ℳ(k)√(η_m,l)L_k^HH_k,m^HQ_l,m  .

and 

    R_k=∑_l≠ kA_k,lA_k,l^H+σ_z^2L_k^HL_k  ,

is the covariance matrix of the interference in the signal received at the k-th MS.
Now, some algebraic manipulations are needed in order to express the achievable rate ℛ_k 
in (<ref>) in a form that permits solving the optimization problems considered in Sections 
<ref> and <ref>.
By letting:

    B_k,l,m=L_k^HH_k,m^HQ_l,m⟹A_k,l=∑_m∈ℳ(k)√(η_m,l)B_k,l,m

the covariance matrix (<ref>)  can be rewritten as

    R_k   =∑_l≠ k∑_m∈ℳ(l)∑_m'∈ℳ(l)√(η_m,lη_m',l)B_k,l,mB_k,l,m'^H+
       +σ_z^2L_k^HL_k  ,

thus implying that the rate for the k-th MS can be  expressed as: 

    ℛ_k=Blog_2I+R_k^-1∑_m,m'√(η_m,kη_m',k)B_k,k,mB_k,k,m'^H

Finally, exploiting the fact that the determinant of the product of two matrices factorizes into  the product of determinants, 
we  have:

    ℛ_k   =Blog_2σ_z^2L_k^HL_k+∑_l^K∑_m∑_m'√(η_m,lη_m',l)B_k,l,mB_k,l,m'^H_g_1(η)
       -Blog_2σ_z^2L_k^HL_k+∑_l≠ k^K∑_m∑_m'√(η_m,lη_m',l)B_k,l,mB_k,l,m'^H_g_2(η)  ,

where the vector η is a compact notation to denote the set of all the downlink transmit powers η_m, k, for all the values of k and m.
All previous formulas concern the UC approach but they can be extended to the CF approach with ordinary efforts.



 §.§ Uplink data transmission

The third phase of the communication protocol amounts to uplink data transmission. 
We  denote by x_k^UL(n) the P-dimensional data vector to be transmitted by the k-th MS in the n-th sample time;
the corresponding signal received at the m-th AP is expressed as: 

    y_m(n)=∑_k=1^K√(η_k)H_k,mL_kx_k^UL(n)+w_m(n)

where 
η_k=P_t,k^UL(L_k^HL_k),
and P_t,k^UL is the uplink transmitted power by the k-th MS.



  §.§.§ CF approach
 in this case, each AP forms the following statistic, ∀ k: 

    ỹ_m,k(n)=Q_k,m^Hy_m(n)   ,

wherein now the ZF beamformer Q_k,m previously defined is used as a post-coder.
Then each AP sends to the CPU the vectors ỹ_m,k(n) via the backhaul link, and the CPU forms the following soft estimate of the data vectors transmitted by th k-th MS:

    x̂_k^UL(n)=∑_m=1^Mỹ_m,k(n),     k=1,…, K.

Plugging (<ref>) and (<ref>) into (<ref>) it can be checked that the signal contributions associated to the data-symbol  x_k^UL(n) are coherently summed.


  §.§.§ UC approach
 in this case the signal transmitted by the k-th MS is decoded only by the APs belonging to the set ℳ(k). Accordingly, the CPU performs the following soft estimate: 

    x̂_k^UL,UC=∑_m∈ℳ(k)ỹ_m,k(n),     k=1,…,K,

Substituting (<ref>) and (<ref>) into (<ref>), we get: 

    x̂_k^UL,UC   =∑_m ∈ M(k)Q_k,m^H√(η̃_k)H_k,mL_kx_k^UL(n)+
       +∑_j=1, j≠k^K√(η̃_j)∑_m ∈ M(k)Q_k,m^HH_j,mL_jx_j^UL(n)+
       +
    ∑_m ∈ M(k)Q_k,m^Hw_m(n)   .

Given the above equation, the achievable rate for the k-th MS can be easily shown to be expressed in the UC case as

    R_k =B log_2 | 𝐈+ η_k 𝐑_k^-1(∑_m ∈ M(k)𝐁_k,k,m^H)(∑_m ∈ M(k)𝐁_k,k,m)|   ,

wherein the interference covariance matrix now is written as

    [                                               𝐑_k= ∑_j ≠ kη_j(∑_m ∈ M(j)𝐁_j,k,m^H)(∑_m ∈ M(j)𝐁_j,k,m);                                                   ;                                                                                   +σ^2_w ( M(k) )   , ]

with (·) denoting cardinality. The above expression can be used also for the CF case by letting M(k)={1, 2, …, M}, for all k=1, …, K.




§ GLOBAL ENERGY EFFICIENCY MAXIMIZATION

In this section we address the problem of power control for energy efficiency maximization considering both the downlink and the uplink channel of the described cellular network. The downlink case is treated in the coming Section <ref>, while the uplink scenario is dealt with in Section <ref>.



 §.§ Downlink power control

Mathematically the problem is formulated as the optimization program:

	
    max_η∑_k=1^Kℛ_k(η)∑_m=1^M[∑_k ∈𝒦(m)δη_m,k + P_c,m]
       s.t.∑_k ∈𝒦(m)η_m,k≤ P_max,m, ∀ m=1,…,M 
           η_m,k≥ 0, ∀ m=1,…,M, k=1,…,K

	

where P_c,m>0 is the circuit power consumed at AP m, δ≥1 is the inverse of the transmit amplifier efficiency, and P_max,m is the maximum transmit power from AP m; η is a KM× 1 vector containing all the transmit power of all AP.

Problem (<ref>) is challenging due to its fractional objective, which has a non-concave numerator. This prevents the direct use of standard fractional programming methods such as Dinkelbach's algorithm to solve (<ref>) with affordable complexity. In addition, another issue is represented by the large number of optimization variables, i.e. KM. To counter both issues, we resort to the successive lower-bound maximization method, introduced in <cit.>[In <cit.> the method is termed successive upper-bound minimization since there the focus is on minimization problems.] (it has already been used in other publications, <cit.> and <cit.>), and whose details are summarized in the Appendix. In brief, this method tackles (<ref>) by alternatively optimizing the transmit powers of each AP, while keeping the transmit powers of the other APs fixed. However, as it will be shown, even with respect to the transmit powers of a single AP, the global energy efficiency problem remains challenging and for this reason, each subproblem is tackled by means of sequential optimization, as described next.

Consider thus problem (<ref>) and define the variable blocks 
η_m={η_m,k}_k=1, …, K, for the CF case, and 
η_m={η_m,k}_k ∈ K(m), for the UC case; 
the vector η_m is K-dimensional in the CF case, while in the UC approach it will contain N entries, i.e. the powers to be used to communicate with the MSs in the set K(m). The global energy efficiency maximization problem with respect to the m-th variable block refers to the maximization, with respect to η_m of the quantity

    ∑_k=1^Kℛ_k(η_m,η_-m)∑_m=1^M[∑_k ∈𝒦(m)δη_m,k + P_c,m]  ,









with η_-m representing the set of all the APs' transmit powers except the ones in η_m.
It can be easily seen that the numerator of the global energy efficiency is not concave even with respect to only the variable block η_m.
In particular, (<ref>) has a non-concave numerator, since ℛ_k(η), which depends on all transmit powers of the base stations, is the difference of two logarithmic functions with the powers η_m of AP m appearing in both logarithmic terms, and the difference of two logarithms is in general not concave.
 This means that (<ref>) can not be solved with affordable complexity by either standard convex optimization nor fractional programming tools. Therefore, plain alternating maximization is not suitable, which motivates the use of the successive lower-bound maximization method. To this end, let us observe that each summand at the numerator of the global energy efficiency can be expressed as the difference between the functions g_1 and g_2 defined in (<ref>). At this point we offer the following result.

	Both functions g_1 and g_2 in (<ref>) are concave in η_m 


As a first step, let us show that the function f: (x,y)∈{ℝ_0^+×ℝ_0^+}→√(xy)∈ℝ_0^+ is jointly concave in (x,y). To this end, the Hessian of f can be written as

    H=1/4([ -x^-3/2y^1/2 x^-1/2y^-1/2; x^-1/2y^-1/2 -x^1/2y^-3/2 ]) ,

which can be seen to be negative-semidefinite, since its determinant is zero and the elements on the diagonal are non-positive. 

Next, we observe that both g_1 and g_2 in (<ref>) are obtained by composing the function f with the log_2 |(·)| function of positive-definite argument. Then, the result follows recalling that the function log_2 |(·)| is matrix-concave and matrix-increasing over the set of positive-semidefinite matrices <cit.>. 

Based on Lemma <ref>, we argue that (<ref>) is the difference between two concave function. Thus, recalling that any concave function is upper-bounded by its first-order Taylor expansion around any given point η_m,0, a concave lower-bound of ℛ_k can be obtained as: 


















    ℛ_k   =g_1(η_m)-g_2(η_m)
       ≥ g_1(η_m)-g_2(η_m,0)-∇_η_m^Tg_2|_η_m,0(η-η_m,0)
        =ℛ_k(η-η_m,0).
  
Then, a suitable approximate problem for the implementation of the sequential optimization method is  

	
    max_η_m∑_k=1^Kℛ_k(η_m,η_m,0,η_-m)δ∑_m=1^M∑_k ∈𝒦(m)η_m,k + P_c
       s.t.∑_k ∈𝒦(m)η_m,k≤ P_max,m
           η_m,k≥ 0,   k=1,…,K

	

It can be seen that Problem (<ref>) fulfills the three properties P1, P2 and P3 of the sequential method detailed in the Appendix[The first-order Taylor expansion and its derivative coincide with the function and with its derivative when evaluated at the center of the expansion]. 
Moreover, the numerator of the objective is now concave, which enables the use of fractional programming tools to globally solve (<ref>), such as the popular Dinkelbach's algorithm <cit.>. The overall power control procedure is formulated as in Algorithm <ref> and, based on the properties of the successive lower-bound maximization method reviewed in the Appendix, we can state the following result.

Algorithm <ref> monotonically improves the global energy efficiency value after each iteration and converges to a first-order optimal point of the original global energy efficiency maximization problem in (<ref>).











  §.§.§ An alternative definition for the GEE on the downlink

The definition of the global energy efficiency reported in (<ref>) considers a circuit power consumption that does not depend on the transmit power used by each base station. However, a base station that does not transmit will consume a lower circuit power, since it will switch to idle mode. We stress that such a circumstance can arise when maximizing the energy efficiency, since the marginal increase to the system achievable rate granted by their activation is overweighted by the increased network power consumption.
To account for this circumstance, the terms P_c,m can be further detailed to depend on the actual used transmit power, namely defining for all m=1,…,M,

    P_c,m=P_c,m(1[P_T(m)>0]+0.5(1-1[P_T(m)>0]))

where 
P_T(m)=∑_k η_m,k is the power radiated by the m-th AP and 
1[P_T(m)>0] is the indicator function of the set [P_T(m)>0], being 1 when P_T(m)>0 and 0 otherwise.  According to (<ref>), we assume that the circuit power consumption, equal to P_c,m when the AP is active, is halved when it does not radiate any power. 

While this more sophisticated power consumption model accounts for the lower circuit power consumption of idle access points, it leads to a non-differentiable global energy efficiency function, due to the presence of the indicator function I in (<ref>). In this case Algorithm <ref> is still guaranteed to monotonically increase the global energy efficiency value after each iteration, but no first-order optimality can be guaranteed upon convergence, owing to the non-differentiability of (<ref>).  Nevertheless, we should remark that it is also possible to approximate the indicator function in (<ref>) by a smooth function, such as a sigmoid, thus recovering the first-order optimality property of the algorithm.



 §.§ Uplink power control

With regard to the uplink, 
the problem of global energy efficiency maximization is formulated as the optimization program:

	
    max_η∑_k=1^Kℛ_k(η_1, …, η_K)[∑_k=1^K δη_k + P_c,k] 
       s.t.  η_k(𝐋_k^H 
    	𝐋_k)≤ P_T,max, ∀ k=1,…,K 
              η_k≥ 0, ∀ k=1,…,K   ,

	

where P_T,max denotes the maximum transmit power for the MSs[For the sake of simplicity this power is assumed to be the same for all the MSs; this assumption however can be easily relaxed.]. 
It is easy to realize that the optimization problem (<ref>) has the same structure as (<ref>), thus implying that it can be solved by using the same procedure that has been developed for the downlink. It is also worth noting that while in the downlink the number of variables to be optimized is  MK in the CF case and MN in the UC case, for the uplink only K transmit powers are to be optimized, thus implying that the power optimization is much less computationally intensive on the uplink than on the downlink. 




§ SUM-RATE MAXIMIZATION


While the main focus of this work is on the maximization of the global energy efficiency function, it should be emphasized that the same approach can be adopted also for sum-rate maximization, which can be tackled as a special case of global energy efficiency maximization. Nonetheless, as detailed in the following, there are instances, when perfect CSI is available, in which the sum-rate maximization algorithm reduces to a standard convex optimization problem. 
Specifically, the global energy efficiency function reduces to the sum-rate function upon plugging δ=0 and P_c=1 in (<ref>). Thus, in principle, sum-rate maximization can be performed by applying the simplified version of Algorithm <ref> wherein standard convex optimization routines can be used to solve Problem (<ref>) when μ=0 and P_c=1 in each iteration[Note that when μ=0 and P_c=1, (<ref>) is no-longer a fractional program.]. 




















Nevertheless, sum-rate maximization has one peculiarity that should be highlighted with respect to the global energy efficiency maximization scenario. With reference to the case in which perfect CSI is available, it is possible to consider the ZF precoder 

    R_k=σ_zL_k^HL_k=N_MS/Pσ_zI ,

thus removing multi-user interference. As a result, the sum-rate function simplifies to: 

    ∑_k=1^Kℛ_k   =∑_k=1^KBlog_2|
    I+(N_MS/Pσ_z^2I)^-1
       ∑_m,m'∈ℳ(k)√(η_m,kη_m',k)C_k,m,m'|

where C_k,m,m' is: 

    C_k,m,m'=L_k^HH_k,m^HQ_k,mQ_k,m'^HH_k,m'L_k ,

which is a concave function by virtue of Lemma <ref>. In this case, since the sum-rate is already concave and can be globally maximized with polynomial complexity by using standard convex programming theory.



 



 

 





§ PERFORMANCE ANALYSIS





 §.§ Simulation setup

We start by defining the used channel model for generating the matrices 𝐇_k,m.
 According to the widely used clustered channel model for mmWave frequencies (see <cit.> and references therein), H_k,m can be expressed as 

    H_k,m= γ∑_i=1^N_cl∑_l=1^N_rayα_i,l√(L(r_i,l))a_AP(θ_i,l,k,m^AP)a_MS^H(θ_i,l,k,m^MS) + H_LOS,

where N_cl is the number of clusters, N_ray is the number of rays that we consider for each cluster, γ is a normalization factor defined as √(N_APN_MSN_clN_ray), 
H_LOS is the line-of-sight (LOS) component, α_i,l is the complex path gain distributed as 𝒞𝒩(0,1), so that its amplitude is Rayleigh-distributed, L(r_i,l) is the attenuation related to the path (i,l), a_AP and a_MS are the array responses at the m-th AP and at the k-th MS, respectively,  and they depend on the angles of arrival and departure, θ_i,l,k,m^AP and θ_i,l,k,m^MS, relative to the 
(i,l)-th path of the channel between the k-th MS and the m-th AP.
The path-loss is defined as <cit.>

    L(r)=-20log_10(4π/λ)-10nlog_10(r)-X_σ,

wherein r is the distance between the transmitter and the receiver, n is the path loss exponent, X_σ is the shadow fading term in logarithmic units with zero mean and σ^2-variance and f_0  is a fixed frequency (see also table <ref>). The 
H_LOS in (<ref>) is written as[For the ease of notation we omit the subscript k,m.] 

    H_LOS=I(d)√(N_APN_MS)e^jη√(L(d))a_AP(θ_LOS^AP)a_MS^H(θ_LOS^MS).

In the above equation, η∼𝒰(0,2π), I(d) is a 0-1 random variate indicating if a LOS link exists between the transmitter and the receiver, and d is the transmitter-receiver distance, measured in meters. Denoting by p the probability that I_LOS(d)=1, we have, for the  UMi (Urban Microcellular) scenarios <cit.>: 

    p=min(20/d,1) (1-e^-d/39)+e^-d/39.
 
As for the number of scatterers and their positions, it should be said that, while usually for every AP-MS pair, a random and independently generated set of scatterers is considered to contribute to the channel matrix (<ref>), in this paper, in order to model the possible channel correlation when the devices are closely spaced, we consider the same set of scatterers for the generation of all the channels. In particular, we assume that, in the considered area, there is a given number of random clusters, each one contributing with three rays[As specified in Section <ref>, a density of 0.4 cluster/sqm. will be considered.].  Given these clusters, in order to generate the generic channel H_k,m between the k-th MS and the m-th AP, we consider as active only those clusters falling in an ellipse built around the position of  the MS and the AP (see Fig. <ref>). In this way, on one hand we exclude far clusters from contributing to the channel, while, on the other hand, we are guaranteed that devices closely located will have correlated channels since they will be affected by similar sets of scatterers (see Fig. <ref> for a 
graphical illustration).





In the following simulation we have considered a carrier frequency of f_0=73 GHz, a bandwidth of B=200 MHz and, with regard to the above channel model,  we have simulated the UMi Open Square scenario <cit.> of size 250 × 250  sqm. In order to generate the correlated channels, a total  of 25.000 randomly deployed clusters (corresponding to a cluster density of 0.4 cluster/sqm.)  has been generated. The additive white noise at the receiver has a power spectral density of -174 dBm/Hz and the receiver noise figure has been set to F=6 dB. The simulated system has M=80 APs randomly deployed in the area to cover and equipped with an ULA with N_AP=16 antennas each. We consider both a scenario with K=6 users, that can be representative of a lightly-loaded network, and a scenario with K=16, that can be representative of a heavily-loaded scenario.  We plot results for the UC networking deployment with N=1 and N=3.
Each MS is equipped with and ULA with N_MS=8 antennas. Although the illustrated algorithms enable each MS  to transmit and receive multiple data-streams, a multiplexing order of P=1 is assumed here for the sake of simplicity. The presented results show the global energy efficiency [Mbit/Joule] as defined in (<ref>), with the circuit power consumption expressed as in (<ref>),   and the achievable rate per user [bit/s], defined as the sum-rate (i.e. the numerator of the global energy efficiency) divided by the number of MSs K.
For benchmarking purposes, we 
compare the performance achieved by the proposed power control rules (maximizing the global energy efficiency and the sum-rate) with 
that achieved by downlink uniform power control, which assumes that  
η_m,k=P_T/K
in the CF case and that 

    η_m,k=
    P_tcard(𝒦(m)) ,    k ∈𝒦(m)
    
    0,    k ∉𝒦(m)

for the UC approach, respectively.
The numerical values come from an average over 1000 independent channel scenarios as well as users and access points locations. 
The convergence criterion for the proposed sequential optimization algorithm is based on the computation of the relative tolerance, i.e. the norm of the difference between the current optimized vector and the optimized vector available at the previous iteration, divided by the norm of the current optimized vector.  The successive convex approximation routine was directly implemented by the authors with a Matlab script, whereas the routine for maximizing the sum-rate in the case of concave rate was fmincon.
All APs have been assumed to have the same maximum feasible transmit power P_max on the downlink. The transmit amplifier efficiency of each transmitter has been assumed equal to one, i.e. δ=1, while the hardware circuit power was modeled according to the model in (<ref>) for each AP, with 
P_c,m=1 W, for all m=1, …, M. For the uplink, we instead used P_c,k=0.3 W, for all k=1, …, K.
In the figures, the dashed line represent the results obtained by means of optimization procedure, instead, the solid line (it has not been reported in the legend in order to avoid a redundancy), with the same color, represent the case of uniform power allocation.  



 §.§ Numerical results



Fig. <ref> compares the global energy efficiency value versus P_max achieved by the proposed 
global energy efficiency-maximizing power control scheme (labeled as OPT), considering the UC approach in the following scenarios:


  * Perfect CSI and FD beamforming. 

  * Perfect CSI and HY beamforming, with 4 RF chains used in the BCD-SD HY beamforming algorithm. 

  * Imperfect CSI and FD beamforming, with pilot sequences of length τ_p=64 and uplink transmit power of 100 mW.

  * Imperfect CSI and HY beamforming, with pilot sequences of length τ_p=64 and uplink transmit power of 100 mW, with 4 RF chains used in the BCD-SD HY beamforming algorithm.

The plot on the left refers to the UC rule wth N=1, while on the right we have the case N=3. 
Fig. <ref> considers the same setting as Fig. <ref> but reports the downlink average rate per user achieved by the proposed rate-maximizing power control scheme (labeled as OPT). 
Figs. <ref> and <ref> report the same results as Figs. <ref> and
<ref>, respectively, with the difference that they refer to the highly-loaded scenario, i.e. with K=16 users. 
Inspecting the figures, several considerations can be made. First of all, we see that the proposed power optimization method provides better performance than the uniform power allocation scheme. This is not always true for the case with imperfect CSI (see plot on the left in Fig. <ref>), as a consequence of the fact that, in the imperfect CSI case, the optimization step was performed using the estimated channels, while the plotted curve represents the true global energy efficiency, computed using the real channel coefficients.  In other words, as far as the imperfect CSI scenarios are concerned, the metric that is optimized is different from the metric that is plotted. Next, as expected, we notice that FD beamforming and the availability of perfect CSI lead to better performance with respect to the practical situation that imperfect CSI and HY beamforming is to be accounted for. For instance, focusing on the plot on the left in Fig <ref> and considering a maximum transmit power of 0 dBW, the average rate-per-user drops from 2.65 Gbit/s of the ideal case  to 1.65 Gbit/s, 
for the case in which both HY beamforming and ICSI are taken into account. Further, when comparing the lightly-loaded scenario with the highly-loaded one, we see that the system global energy efficiency does not dramatically change in the two situations, while, conversely, in a heavily-loaded system with HY beamforming and incomplete CSI the average rate-per-user at 0dBW of maximum transmit power is approximately equal to 500 Mbit/s, with about 60%  loss with respect to the average rate-per-user attainable in a lightly loaded scenario.
A peculiar behavior that is observed in Figs. <ref> and <ref> is that in the case of perfect CSI with FD beamforming with N=1, the uniform and optimal power allocation schemes achieve the same performance. This is due to the fact that the rate is an increasing function of the total available power P_max and so, when only one user needs to be served, uniform power allocation coincides with the rate-maximizing power allocation strategy. Instead, this behavior is not observed in  Figs. <ref> and <ref>, where a visible gap is present between the performance with uniform power allocation and with the optimized power allocation also when N=1. This is expected because, unlike the rate, the global energy efficiency is not monotonically increasing with P_max.





Next, Fig. <ref> 
 shows the CDFs for the rate-per-user in case of FD and HY beamforming, respectively, considering P_max=0 dBW and the lightly loaded scenario. 
Again the two situations N=1 and N=3 are examined.  The curves confirm the findings previously commented. Additionally, it can be seen that for the practical case of incomplete CSI and HY beamforming the scenario with N=1 provides much better performance than the one with N=3 for the vast majority of the users, i.e. it appears to be convenient to have each AP serve just one user than 3 users, except than when considering the unlucky users that happen to be situated on the left tail of the rate-per-user distribution. This behavior can be explained by noticing that in the considered setting, where M>>K, the majority of the MSs has in its neighborhood several APs, and letting these APs dedicate their own resources to only MS results in overall increased performance. 
 







Finally, we turn our attention to the performance of the uplink channel. Fig. 
<ref> 
reports the average uplink rate-per-user versus the uplink transmit power. No power control is performed here. Again, the scenarios N=1 and N=3 are considered. From the figure, it is seen that also for the uplink the rate is a clearly increasing function of the transmit power only in the ideal case of perfect CSI and fully digital beamforming, while, in the other situations, a saturation effect is observed. Results also show that, when focusing on the practical scheme with incomplete CSI and HY beamforming the case N=1 provides better results than the case N=3. In particular, assuming an uplink transmit power of -10dBW a data-rate of about 1.2 Gbit/s can be achieved for the case N=1, versus a data-rate of about 600 Mbit/s for the case N=3. Again, this behavior can be explained by noticing that in a scenario with a dense AP deployment each MS is sorrounded by several APs, and it is better to let each AP to focus its own resourced on just one MS rather then letting the AP to share its resources among multiple MSs. 





§ CONCLUSION


This paper has considered a CF massive MIMO system operating at mmWave frequencies with a UC association between APs and MSs. Adopting a clustered channel model capable of taking into account the channel correlation for nearby devices, the paper has analyzed both the CF and UC approaches, by proposing a low-complexity power allocation rule aimed at global energy efficiency maximization. The presence of HY beamforming architectures at the APs has been considered, whereas, for the MSs, a simple channel-independent 0-1 beamforming structure has been considered. The obtained results have confirmed that the proposed resource allocation algorithms are effective at increasing the system energy efficiency and the system average rate-per-user,  as well as that the use of HY beamforming architectures introduces, as expected,  a considerable performance degradation. 
To the best of the authors' knowledge, this is the first paper to consider the CF and UC massive MIMO architectures at mmWave frequencies. This study can be generalized along several paths. First of all, a comparison with a massive MIMO deployment with co-located antennas should be carried out. Then, given the gap observed between the HY and the FD beamforming schemes, more sophisticated HY beamformers should be considered to improve the system performance. Finally, given the intrinsic macro-diversity guaranteed by the dense AP deployment, it would be of interest to investigate on the suitability of CF architectures to circumvent the blockage effects that cause link instability: in this sense, CF massive MIMO could be envisioned as an architecture capable  of making mmWave frequencies useful for ultra-reliable communications in indoor and small-sized outdoor environments. 






§ APPENDIX: SUCCESSIVE LOWER-BOUND MAXIMIZATION

In this appendix we briefly review the successive lower-bound optimization technique. 
The method is based on the  idea of merging the tools of alternating optimization <cit.> and sequential convex programming <cit.>. Specifically, let us consider the generic optimization problem 

    max_x∈ X f(x) ,

with f:ℝ^n→ℝ a differentiable function, and X a compact set. Similarly to the alternating optimization method, the successive lower-bound maximization partitions the variable space into J blocks, x=(x_1,…,x_J), which are cyclically optimized one at a time, while the other blocks are kept fixed. Thus, Problem (<ref>) is decomposed into M subproblems, wherein the generic subproblem is stated as

    max_x_m f(x_m,x_-m) ,

with x_-m cotaining all variable blocks except the m-th. If Problem (<ref>) is globally solved for each m=1,…,J, then we have an instance of the alternating maximization method, which, as proved in <cit.>, monotonically improves the objective of (<ref>), and yields a first-order optimal point if the solution of (<ref>) is unique for any m, and if X= X_1×…× X_J, with x_m∈ X_m for all m. However, if globally solving (<ref>) is difficult (e.g. because (<ref>) is not a convex problem), then implementing the alternating maximization method proves difficult. In this case, the successive lower-bound maximization method proposes to find a (possibly suboptimal) solution of (<ref>), by means of the sequential convex programming method. Besides leading to a computationally viable algorithm, this approach is proved to preserve the optimality properties of the alternating optimization method <cit.>, despite the fact that a possible suboptimal solution of (<ref>) is determined. 

As for sequential optimization, its basic idea is to tackle a difficult maximization problem by a sequence of easier maximization problems. To elaborate, denote by g_i(x_m) the i-th constraint of (<ref>),  for i=1,…,C, and consider a sequence of approximate problems { P_ℓ}_ℓ with objectives {f_ℓ}_ℓ and constraint functions {g_i,ℓ}_i=1^C, such that the following three properties are fulfilled, for all ℓ:

	
  (P1) f_ℓ(x_m)≤ f(x_m), g_i,ℓ(x_m)≤ g_i,ℓ(x_m), for all i and x_m;
	
  (P2) f_ℓ(x_m^(ℓ-1))=f(x_m^(ℓ-1)), g_i,ℓ(x_m^(ℓ-1))=g_i(x_m^(ℓ-1)) with x_m^(ℓ-1) the maximizer of f_ℓ-1;
	
  (P3) ∇ f_ℓ(x_m^(ℓ-1))=∇ f(x_m^(ℓ-1)), ∇ g_i,ℓ(x_m^(ℓ-1))=∇ g_i(x_m^(ℓ-1)).

In <cit.> (see also <cit.>) it is shown that, subject to constraint qualifications, the sequence {f(x_m^(ℓ))}_ℓ of the solutions of the ℓ-th Problem P_ℓ, is monotonically increasing and converges. Moreover, every convergent sequence {x_m^(ℓ)}_ℓ attains a first-order optimal point of the original Problem (<ref>). Thus, sequential optimization guarantees at the same time to monotonically improve the objective function, and to fulfill the kkt first-order optimality conditions of the original problem. Nevertheless, in order to be able to use the method, it is necessary to find lower bounds of the original objective function, which fulfill all three properties P1, P2, P3, while at the same time leading to an approximate problem that can be solved with affordable complexity.


IEEEtran





