\section{Proof of Theorem}
\begin{proof}
Proof goes here.
\end{proof}
\section{Conclusion}
\label{Sec:Conc}
We propose a learning-based \gls{mpc}\xspace policy for user association in multi-band heterogeneous networks with mobile users.
Using the proposed novel loss function---linear combination of \gls{nmse}\xspace and \gls{ndcg}\xspace loss---the forecaster-aided approach is a sample-efficient and interpretable alternative to model-free \gls{rl}\xspace.
We show that the service rates for the edge users can be greatly improved and our mobility-aware policy significantly reduces the handover\xspace rate.
Moreover, using our implementation of the Frank-Wolfe method, the \gls{mpc}\xspace problem can be solved efficiently motivating the usefulness of our approach in a real deployment.
There is also considerable scope for further work.
For example, designing a forecaster robust to out-of-distribution user mobility patterns, testing the performance on more sophisticated BS deployments with multiple tiers and more frequency bands, and a more rigorous study that models the handover\xspace mechanism beyond the physical layer.
Extending the proposed techniques to MIMO systems and higher frequency bands (for example, millimeter wave and sub-terahertz) to study the impact of beamforming and high pathloss would also be suitable for the \gls{mpc}\xspace framework.

\section{Acknowledgements}
\label{Sec:Acks}
The authors would like to thank Ruichen Jiang at UT Austin for the insightful discussions on optimization algorithms which lead to the development of our computationally efficient optimization algorithm and for his very helpful feedback on the Frank-Wolfe method which significantly reduced the simulation time.
\section{Introduction}
\label{Sec:Intro}
Cellular networks continue to trend towards multi-tier and multi-band deployments to meet intense consumer demands for faster and more ubiquitous data connectivity.
The use of low-power small cell technology enables flexible deployments and helps boost network capacity through increased spatial reuse.
Meanwhile, more spectrum is being released, for instance almost 500 MHz recently in the 3.7 - 4.2 GHz ``C-band" and more to be auctioned in the ``beachfront'' bands from 1 to 2.6 GHz.  Assigning users to bands and cells is an increasingly complicated problem, with large implications on the overall network coverage and throughput, particularly in the context of mobility.  The conventional strategy of simply maximizing \gls{snr}\xspace or \gls{sinr}\xspace for user association tends to lead to an imbalanced load distribution and can severely degrade both a given user's as well as the overall service rate \cite{Seven_Andrews13}.  Our goal in this paper is to develop an optimal band assignment and \gls{bs}\xspace selection strategy for multi-band heterogeneous networks that accounts for the bandwidth and propagation characteristics of each band, along with the traffic load,  user mobility, and handover overhead.

\subsection{Motivation, Background, and Scope of Paper}
A multi-band deployment increases the overall bandwidth as well as number of connection options for subscribers, however, it comes at the cost of increased complexity, especially in the context of mobility \cite{andrews14overview}.
Traditional \emph{user association}---meaning a user's \gls{bs}\xspace selection and band assignment---policies that maximize received signal power or \gls{sinr}\xspace lead to preferential association towards the lowest frequency bands---which offer better penetration and lower pathloss---and to tower-mounted macrocells with higher transmit power.
Such association rules are highly sub-optimal in terms of throughput, because higher carrier frequencies typically offer wider bandwidths while small cells are much more lightly loaded.
Furthermore, heterogeneous deployments with small coverage areas are more prone to frequent handovers\xspace and consequently poor service rates due to transmission delays, ping-pong effects, and connection drops \cite{Arshad_Velocity17, feriani22multiobjective}.

As a result, there has been considerable work investigating user association and cell load balancing.
These include stochastic geometry studies which analytically characterize the throughput gain from load balancing \cite{singh13offloading, kassir22analysis, Modeling_Lin13} and optimization-based user association strategies \cite{shen14distributed, ye13user, kim12distributed, ye16user, liu19joint, teng21joint}.
In \cite{ye13user, shen14distributed}, authors propose log-utility-based association strategies to achieve network-wide proportional fairness for static users, and \cite{kim12distributed} presents $\alpha$-optimal load balancing strategy for inhomogeneous user traffic.
The user association problem has also been studied in conjunction with power allocation \cite{liu19joint}, massive MIMO \cite{ye16user}, and \gls{bs}\xspace activation strategies \cite{teng21joint}, under the network utility maximization framework.

Accounting for user mobility within the optimization framework is generally intractable, with limited studies in the literature, mostly using simplistic mobility models.
However, a variety of studies underscore the benefits of mobility-aware user association and handover\xspace management \cite{addali19dynamic, lin13towards, hsueh17equivalent, Arshad_Velocity17}.
The authors in \cite{lin13towards, hsueh17equivalent} characterize handover\xspace rate as a function of user velocity and \gls{bs}\xspace density.
By leveraging similar stochastic geometry-based analysis, \cite{Arshad_Velocity17} proposes a \gls{bs}\xspace skipping handover\xspace policy for high-speed users in a two-tier cellular network.
For a comprehensive survey of mobility models---both trace-based and probabilistic---and mobility-aware performance characterizations, please refer to \cite{tabassum19fundamentals}.

Meanwhile, in actual cellular deployments, proactive load balancing is usually accomplished by biasing user association towards small cells \cite{rakotomanana16optimum, ghosh12heterogeneous} by some tunable dB value.
For example, a small cell could be assigned a 5-10 dB \gls{snr}\xspace bias relative to a macrocell \cite{QualcommHeNets}.
However, such rule-based strategies are not only suboptimal, but they require carefully tuning the bias values and struggle to cope with time-varying traffic patterns from mobile users and for large number of connections \cite{lin22embracing}.
As the industry trends towards intelligent and open networks \cite{polese22understanding, lin22embracing}, network measurements across \glspl{bs}\xspace and users can be leveraged to make data-driven user association and load balancing decisions.
As such, mobility-aware user association has been considered recently as a \gls{rl}\xspace problem \cite{zhang18loadbalancing, Handover_Wang18, chinchali18cellular, zhao19deep, xu19load, gupta21load, wu21load, feriani22multiobjective, lacava22programmable}.
Among these, \cite{zhao19deep, gupta21load} aim to optimize the network utility within the \gls{rl}\xspace framework, \cite{zhang18loadbalancing, Handover_Wang18} maximize the sum-rate, and \cite{xu19load, feriani22multiobjective} optimize cell individual offset parameters.
However, such model-free \gls{rl}\xspace techniques are known to be sample-inefficient \cite{mai2022sample}, and there has been a push towards scalable and sample-efficient formulations \cite{wu21load, alcaraz22model,agarwal2022taskdriven}.
The authors in \cite{wu21load} use transfer learning and a model-based \gls{rl}\xspace solution was explored in \cite{alcaraz22model} for load balancing and network optimization.

In this work, we focus on the most widely deployed frequency bands, namely the ``sub-6 GHz" or ``frequency range 1 (FR1)" range, which still carry the vast majority of cellular traffic even in 5G deployments today, despite the introduction of millimeter wave bands.
At these lower frequencies, initial access and control channel communication can be reliably performed without the requirement for highly directional beamforming, which greatly simplifies the user association and handover\xspace decisions.
A key challenge for load balancing with highly directional transmissions is that it leads to directional and time-varying interference, load-based association decisions problematic and non-convex.
Furthermore, since such bands will be lightly loaded for the forseeable future, it may not be necessary to perform load balancing on them.
Thus, we restrict our attention to the lower bands in this work.

\subsection{Contributions}

The overarching technical contribution of this paper is the development of a novel framework for mobility-aware \gls{bs}\xspace and frequency band association in a large-scale cellular network.  The key elements and unique aspects of our approach are summarized as follows.

\textbf{Optimization framework for mobility and handover\xspace aware association.}
We formulate a novel optimization problem to maximize network utility based on the users' service rates.
In \secref{subsec:service_rate} we define the user service rate, which depends on the load on the serving \gls{bs}\xspace, the \gls{sinr}\xspace between the user and \gls{bs}\xspace, and the signaling overhead in case of a triggered handover\xspace.
A key insight is that in a mobile network, the user association task can be decomposed into forecasting user trajectories and thus their future rates, and then a convex utility maximization problem based on those predicted rate values.
As such, we propose a learning-based \gls{mpc}\xspace approach in \secref{sec:mpc_formulation}, where we only need data-driven learning for forecasting, while we can employ interpretable utility maximization for band assignment and \gls{bs}\xspace selection.

\textbf{Sample-efficient forecaster design.}
In \secref{sec:forecaster_design}, we design a forecaster to model and predict the users' mobility patterns, which are then used by the \gls{mpc}\xspace controller to solve for associations that maximize the network utility.
The forecaster is trained in a supervised fashion, which is generally more sample-efficient and stable compared to model-free \gls{rl}\xspace.
To further improve the sample-efficiency, in we identify important characteristics of a ``good" forecast for the given utility maximization problem and develop a novel loss function to train a task-aware forecaster, instead of using the task-agnostic mean-squared error loss.
We compare the forecaster-based approach with \gls{rl}\xspace in terms of generalization ability and sample-efficiency and show that the proposed approach can learn the system dynamics induced by user mobility with $100\times$ less training samples.

\textbf{Computationally efficient optimization algorithm.}
In \secref{sec:frank_wolfe_algo}, we propose an optimization algorithm to quickly solve the constrained utility maximization problem for user association.
We exploit the structure of the constraint set to derive an analytical solution for the direction-finding step in the projection-free Frank-Wolfe method.
Leveraging the analytical solution and backtracking line-search to dynamically adjust the step-size in every iteration, the execution time for solving the utility maximization can be reduced by two orders of magnitude compared to a generic convex solver.

\textbf{Realistic validation of rate gain and handover\xspace reduction.}
In \secref{sec:results}, we evaluate the performance of the proposed forecaster-aided \gls{mpc}\xspace solution through extensive simulations using an industry-grade network simulator developed by Meta (formerly Facebook) that accurately models a commercial cellular network.  We demonstrate that our \gls{mpc}\xspace approach can provide a more uniform user experience and improve the 5th percentile service rate by about $3.5\times$ compared to conventional strategies.
Furthermore, we show that by using our forecaster-based mobility-aware user association strategy, the median number of handovers\xspace in a time slot can be simultaneously reduced by about $7\times$ compared to a handover\xspace agnostic strategy.

\textit{Notation}: ${\mathbf{A}}_{:,i}$ denotes the $i$-th column of matrix ${\mathbf{A}}$, ${\mathbf{A}}_{i,:}$ denotes the $i$-th row of matrix ${\mathbf{A}}$, ${\mathbf{A}}^T$ denotes the transpose of a matrix, $[{\mathbf{A}}]_{i,j}$ denotes the $(i,j)$-th element of ${\mathbf{A}}$, and ${\mathbf{A}}^\dagger$ denotes the Hermitian transpose of a matrix. The $i$-th element of vector ${\mathbf{a}}$ is denoted by $a_i$. $\indicator{\mathcal{A}}$ denotes the indicator function over set $\mathcal{A}$, and $|\mathcal{A}|$ denotes the cardinality of set $\mathcal{A}$.
\section{Introduction}
\label{Sec:Intro}
Cellular networks continue to trend towards multi-tier and multi-band deployments to meet intense consumer demands for faster and more ubiquitous data connectivity.
The use of low-power small cell technology enables flexible deployments and helps boost network capacity through increased spatial reuse.
Meanwhile, more spectrum is being released, for instance almost 500 MHz recently in the 3.7 - 4.2 GHz ``C-band" and more to be auctioned in the ``beachfront'' bands from 1 to 2.6 GHz.  Assigning users to bands and cells is an increasingly complicated problem, with large implications on the overall network coverage and throughput, particularly in the context of mobility.  The conventional strategy of simply maximizing \gls{snr}\xspace or \gls{sinr}\xspace for user association tends to lead to an imbalanced load distribution and can severely degrade both a given user's as well as the overall service rate \cite{Seven_Andrews13}.  Our goal in this paper is to develop an optimal band assignment and \gls{bs}\xspace selection strategy for multi-band heterogeneous networks that accounts for the bandwidth and propagation characteristics of each band, along with the traffic load,  user mobility, and handover overhead.

\subsection{Motivation, Background, and Scope of Paper}
A multi-band deployment increases the overall bandwidth as well as number of connection options for subscribers, however, it comes at the cost of increased complexity, especially in the context of mobility \cite{andrews14overview}.
Traditional \emph{user association}---meaning a user's \gls{bs}\xspace selection and band assignment---policies that maximize received signal power or \gls{sinr}\xspace lead to preferential association towards the lowest frequency bands---which offer better penetration and lower pathloss---and to tower-mounted macrocells with higher transmit power.
Such association rules are highly sub-optimal in terms of throughput, because higher carrier frequencies typically offer wider bandwidths while small cells are much more lightly loaded.
Furthermore, heterogeneous deployments with small coverage areas are more prone to frequent handovers\xspace and consequently poor service rates due to transmission delays, ping-pong effects, and connection drops \cite{Arshad_Velocity17, feriani22multiobjective}.

As a result, there has been considerable work investigating user association and cell load balancing.
These include stochastic geometry studies which analytically characterize the throughput gain from load balancing \cite{singh13offloading, kassir22analysis, Modeling_Lin13} and optimization-based user association strategies \cite{shen14distributed, ye13user, kim12distributed, ye16user, liu19joint, teng21joint}.
In \cite{ye13user, shen14distributed}, authors propose log-utility-based association strategies to achieve network-wide proportional fairness for static users, and \cite{kim12distributed} presents $\alpha$-optimal load balancing strategy for inhomogeneous user traffic.
The user association problem has also been studied in conjunction with power allocation \cite{liu19joint}, massive MIMO \cite{ye16user}, and \gls{bs}\xspace activation strategies \cite{teng21joint}, under the network utility maximization framework.
Meanwhile, in actual cellular deployments, proactive load balancing is usually accomplished by biasing user association towards small cells \cite{rakotomanana16optimum, ghosh12heterogeneous} by some tunable dB value.
For example, a small cell could have a 5-10 dB \gls{snr}\xspace bias relative to a macrocell \cite{QualcommHeNets}.

Accounting for user mobility within the optimization framework is generally intractable, with limited studies in the literature, mostly using simplistic mobility models.
However, a variety of studies underscore the benefits of mobility-aware user association and handover\xspace management \cite{addali19dynamic, lin13towards, hsueh17equivalent, Arshad_Velocity17}.
The authors in \cite{lin13towards, hsueh17equivalent} characterize handover\xspace rate as a function of user velocity and \gls{bs}\xspace density.
By leveraging similar stochastic geometry-based analysis \cite{Arshad_Velocity17} proposes a \gls{bs}\xspace skipping handover\xspace policy for high-speed users in a two-tier cellular network.
For a comprehensive survey of mobility models---both trace-based and probabilistic---and mobility-aware performance characterizations, please refer to \cite{tabassum19fundamentals}.

Mobility-aware user association has been considered recently as a \gls{rl}\xspace problem \cite{zhang18loadbalancing, Handover_Wang18, zhao19deep, xu19load, gupta21load, wu21load, feriani22multiobjective}.
Among these, \cite{zhao19deep, gupta21load} aim to optimize the network utility within the \gls{rl}\xspace framework, \cite{zhang18loadbalancing, Handover_Wang18} maximize the sum-rate, and \cite{xu19load, feriani22multiobjective} optimize cell individual offset parameters.
However, such model-free \gls{rl}\xspace techniques are known to be sample-inefficient \cite{mai2022sample}, and there has been a push towards scalable and sample-efficient formulations \cite{wu21load, alcaraz22model,agarwal2022taskdriven}.
The authors in \cite{wu21load} use transfer learning and a model-based \gls{rl}\xspace solutions was explored in \cite{alcaraz22model} for load balancing and network optimization.

In this work, we focus on the most widely deployed frequency bands, namely the ``sub-6 GHz" or ``frequency range 1 (FR1)" range, which still carry vast majority of cellular traffic even in 5G deployments today, despite the introduction of millimeter wave bands.  At these lower frequencies, initial access and control channel communication can be reliably performed without the requirement for highly directional beamforming, which greatly simplifies the user association and handover\xspace.   A key challenge for load balancing with highly directional transmissions is that the association decision for a given UE has a subsequent affect on the interference, making \gls{sinr}\xspace and load-based association decisions problematic and nonconvex.  Furthermore, since such bands will be lightly loaded for the forseeable future, it may not be necessary to perform load balancing on them.  Thus, we restrict our attention to the lower bands in this work.

\subsection{Contributions}

The overarching technical contribution of this paper is the development of a novel framework for mobility-aware \gls{bs}\xspace and frequency band association in a large-scale cellular network.  The key elements and unique aspects of our approach are summarized as follows.

\textbf{Optimization framework for mobility and handover\xspace aware association.}
We formulate a novel optimization problem to maximize network utility based on the users' service rates.
In \secref{subsec:service_rate} we define the user service rate, which depends on the load on the serving \gls{bs}\xspace, the \gls{sinr}\xspace between the user and \gls{bs}\xspace, and the signaling overhead in case of a triggered handover\xspace.
A key insight is that in a mobile network, the user association task can be decomposed into forecasting user trajectories and thus their future rates, and then a convex utility maximization problem based on those predicted rate values.
As such, we propose a learning-based \gls{mpc}\xspace approach in \secref{sec:mpc_formulation}, where we only need data-driven learning for forecasting, while we can employ interpretable utility maximization for band assignment and \gls{bs}\xspace selection.

\textbf{Sample-efficient forecaster design.}
In \secref{sec:forecaster_design}, we design a forecaster to model and predict the users' mobility patterns, which are then used by the \gls{mpc}\xspace controller to solve for associations that maximize the network utility.
The forecaster is trained in a supervised fashion, which is generally more sample-efficient and stable compared to model-free \gls{rl}\xspace.
To further improve the sample-efficiency, in we identify important characteristics of a ``good" forecast for the given utility maximization problem and develop a novel loss function to train a task-aware forecaster, instead of using the task-agnostic mean-squared error loss.
We compare the forecaster-based approach with \gls{rl}\xspace in terms of generalization ability and sample-efficiency and show that the proposed approach can learn the system dynamics induced by user mobility with $100\times$ less training samples.

\textbf{Computationally efficient optimization algorithm.}
In \secref{sec:frank_wolfe_algo}, we propose an optimization algorithm to quickly solve the constrained utility maximization problem for user association.
We exploit the structure of the constraint set to derive an analytical solution for the direction-finding step in the projection-free Frank-Wolfe method.
Leveraging the analytically solution and backtracking line-search to dynamically adjust the step-size in every iteration, the execution time for solving the utility maximization can be reduced by two orders of magnitude compared to a generic convex solver.

\textbf{Realistic validation of rate gain and handover\xspace reduction.}
In \secref{sec:results}, we evaluate the performance of the proposed forecaster-aided \gls{mpc}\xspace solution through extensive simulations using an industry-grade network simulator developed by Meta (formely Facebook) that accurately models a commercial cellular network.  We demonstrate that our \gls{mpc}\xspace approach can provide a more uniform user experience and improve the 5th percentile service rate by about $3.5\times$ compared to conventional strategies.
Furthermore, we show that by using our forecaster-based mobility-aware user association strategy, the median number of handovers\xspace in a time slot can be simultaneously reduced by about $7\times$ compared to a handover\xspace agnostic strategy.

\textit{Notation}: ${\mathbf{A}}_{:,i}$ denotes the $i$-th column of matrix ${\mathbf{A}}$, ${\mathbf{A}}_{i,:}$ denotes the $i$-th row of matrix ${\mathbf{A}}$, ${\mathbf{A}}^T$ denotes the transpose of a matrix, $[{\mathbf{A}}]_{i,j}$ denotes the $(i,j)$-th element of ${\mathbf{A}}$, and ${\mathbf{A}}^\dagger$ denotes the Hermitian transpose of a matrix. The $i$-th element of vector ${\mathbf{a}}$ is denoted by $a_i$. $\indicator{\mathcal{A}}$ denotes the indicator function over set $\mathcal{A}$, and $|\mathcal{A}|$ denotes the cardinality of set $\mathcal{A}$.
\section{Learning-based Model Predictive Control}
\label{sec:mpc_formulation}
\begin{figure}[t!]
    \centering
    \includegraphics[width=6in,                 keepaspectratio]{Figures/forecaster_mpc.pdf}
    \caption{Forecaster-based MPC design. Current and past capacity measurements (computed from the \gls{sinr}\xspace's reported during earlier slots) are input to the forecaster which outputs a prediction of future capacities for a future horizon $H$. These are then used to generate new user associations which maximize the cumulative log-utility.}
    \label{fig:forecaster_mpc}
\end{figure}
The genie-aided optimization in \eqref{opt:util_max_genie} is non-causal because it requires future rate and \gls{hit}\xspace measurements, both of which are unknown at time $n$.
In this section, we propose a learning-based \gls{mpc}\xspace solution where we forecast system behavior and optimize the forecast to produce the best control decision at the current time \cite{mpc_james}.
For the considered mobile network,
the two main sources of uncertainty are (i) user mobility which affects both $c_{ij}[n]$ and \gls{hit}\xspace and (ii) the implemented handover\xspace procedure which primarily affects \glspl{hit}.

The temporal variations in \gls{sinr}\xspace measurements---consequently in $c_{ij}[n]$---are largely a function of users' mobility patterns.
These are predictable to some extent, and for a reasonable time interval $T_s$ the next step in the trajectory is significantly correlated with the past trajectory \cite{liang03predictive}.
Thus, to exploit these temporal correlations, we use a recurrent deep neural network architecture---specifically a \gls{lstm}\xspace---to forecast the future measurements using past and current samples.
We denote by $\hat{c}_{ij}[n]$ the prediction of $c_{ij}[n]$.

On the other hand, in addition to current channel conditions, the \gls{hit}\xspace also depends on the specific handover\xspace procedure implemented by the operator.
Moreover, a handover\xspace is not merely a physical layer task and involves control signal exchange across layers \cite{tayyab19survey}.
As a result, the \gls{hit}\xspace cannot be reliably predicted from only \gls{sinr}\xspace measurements and would require some degree of feedback from higher layers.
Thus, we instead model \gls{hit}\xspace as a random variable.

Distributing the logarithm over the product and using $\sum_{i}x_{ij}[n] = \ell_j[n]$, the objective \eqref{opt:util_max_genie_obj} can be written as
\begin{align}
    \sum_{n=1}^{H}\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}} x_{ij}[n] \log(c_{ij}[n]) + x_{ij}[n]\log\left(1 - \frac{\tho{i}[n]}{T_s}\right) - \sum_{n=1}^H\sum_{j\in\mathcal{B}}\ell_j[n]\log(\ell_j[n]).
\end{align}
The first term here is the logarithmic utility of the achievable rates, and as discussed earlier, we use a forecaster to predict these for $n>1$.
The third term is the entropy of the load variable $\ell_j[n]$ and is maximized when all \glspl{bs}\xspace are equally loaded.
The second term corresponds to the penalty for each handover\xspace. Recall that $\tho{i}[n]$ is a function of the previous association of user $i$ as well as the new association, and is hard to reliably predict from physical layer \gls{sinr}\xspace measurements.
Thus, we model $\tho{i}[n]$ as a random variable and use the convex surrogate penalty
\begin{align}
    \label{eq:ho_surrogate_penalty}
    \frac{(x_{ij}[n] - x_{ij}[n-1])^2}{2}{\mathbb{E}}\left[\log\left(1 - \frac{\tho{i}[n]}{T_s}\right)\right].
\end{align}
Since, $x_{ij}[n]$ and $x_{ij}[n-1]$ are indicator variables, the quadratic term is $1$ if user $i$ goes through a handover\xspace.
With these modifications, the user association \gls{mpc}\xspace can be expressed as
\begin{subequations} \label{opt:util_max_mpc}
\begin{align}
\underset{\{x_{ij}[n]\}}{\mathrm{maximize}} ~&\sum_{n=1}^{H}\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}} x_{ij}[n] \log(\bar{c}_{ij}[n])  + \frac{\eta}{2}\left(x_{ij}[n] - x_{ij}[n-1]\right)^2 \nonumber\\
&-\sum_{n=1}^H\sum_{j\in\mathcal{B}}\ell_j[n]\log(\ell_j[n])\label{opt:util_max_mpc_obj}\\
\mathrm{subject~to}~\nonumber~\\
& \bar{c}_{ij}[n] = c_{ij}[n]\mathbbm{1}(n=1) + \hat{c}_{ij}[n]\mathbbm{1}(n>1) \label{opt:util_max_mpc_rate_def}\\
& \ell_j[n] = \sum_{u\in\mathcal{U}}x_{uj}[n] \quad \forall j\in\mathcal{B} \label{opt:util_max_mpc_load_def}\\
& \sum_{j\in\mathcal{B}}x_{ij}[n] = 1 \quad \forall i\in \mathcal{U} \label{opt:util_max_mpc_ue_cons}\\
& x_{ij}[n]\in\{0,1\} \quad \forall i\in \mathcal{U},\forall j\in\mathcal{B}. \label{opt:util_max_mpc_indicator_cons}
\end{align}
\end{subequations}
Here, $\eta$ is the expected value defined in \eqref{eq:ho_surrogate_penalty} and is a negative constant.
The association indicators $x_{ij}[0]$ refer to the user association from the previous time slot and are not optimization variables.
The above \gls{mpc}\xspace uses the current measurements for $n=1$ and the forecasts $\hat{c}_{ij}[n]$ for $n = 2, \dots, H$.
The forecasting horizon is a hyperparameter which depends on the quality of the forecaster and the available compute power.
Typically, the reliability of the forecaster degrades as the prediction horizon increases, and as such larger $H$ can lead to undesirable control decisions.
The computational complexity of the \gls{mpc}\xspace increases dramatically with $H$ due to the quadratic handover\xspace penalty; a major decision variable when choosing the right $H$.
The above \gls{mpc}\xspace has a concave objective---sum of an affine function, a negative quadratic, and entropy---with linear constraints and can be solved using a generic convex solver such as \cite{cvxpy1}.
We later present a variation of the Frank-Wolfe algorithm which can much more efficiently solve \eqref{opt:util_max_mpc}.

\section{Forecaster Design}
\label{sec:forecaster_design}
\begin{figure}[t!]
    \centering
    \includegraphics[width=6in,                 keepaspectratio]{Figures/forecaster_arch.pdf}
    \caption{Neural network architecture for the forecaster.}
    \label{fig:forecaster_arch}
\end{figure}
Next, we describe the design of the forecaster which uses the time-series of current and past rates to predict the future measurements.
Our goal is to use this forecast to solve for association indicators $x_{ij}[n]$ that maximize the network utility according to \eqref{opt:util_max_mpc}.
Hence, we design the forecaster to directly predict the achievable rates $c_{ij}[n]$, which are related to the \gls{sinr}\xspace measurements through \eqref{eq:achievable_rate_def}.
However, the task-aware loss function and the supervised learning procedure described in this section can be used to forecast any metric of interest.
For example, one could train an \gls{sinr}\xspace forecaster and then use \eqref{eq:achievable_rate_def} to get a prediction of the user rates if this better suits their particular use case.

We denote by ${\mathbf{c}}_i[n]$ the $|\mathcal{B}|\times 1$ vector of achievable rates of user $i$ from all \glspl{bs}\xspace in $\mathcal{B}$.
To exploit the spatial and temporal correlations in ${\mathbf{c}}_i[n]$, we use an \gls{lstm}\xspace-based architecture shown in \figref{fig:forecaster_arch}.
The input to the forecaster is the history of $N$ rate measurement vectors for a user $i$.
The forecaster then generates a rate prediction for $H$ future slots---also known as the forecasting horizon---denoted by $\hat{{\mathbf{c}}}_i[n+1], \dots, \hat{{\mathbf{c}}}_i[n+H]$.
This architecture has the following benefits.
First, the forecaster predicts the trajectory for each user individually and does not depend on the number of users in the network.
In fact, with this architecture, the learning-based \gls{mpc}\xspace formulation in \eqref{opt:util_max_mpc} can even support a time-varying $\mathcal{U}$.
Second, during both the training and inference phases, measurements across users can be batched together and processed in parallel by the forecaster.

Since ${\mathbf{c}}_i[n]$ are continuous variables, training the forecaster lies within the supervised regression paradigm, where the typical approach is to minimize the mean-squared error (MSE) between the prediction and true value.
However, our goal is to use the forecast from the trained model to solve for the load balancing user association indicators $x_{ij}[n]$.
Hence, we instead design a task-aware training loss function with the following two objectives in mind.
First, for a user $i$, the forecaster must preserve the \textit{ranking} induced by the capacity measurements.
Second, the forecaster must accurately predict the achievable rate between user $i$ and \gls{bs}\xspace $j$, relative to other \glspl{bs}\xspace.
In order words, it is satisfactory for the forecaster to be off, by say $1$ Mbps, as long as the forecast for every relevant \gls{bs}\xspace is off by approximately the same amount.
As such, we use the following loss function, which is a sum of the \gls{nmse}\xspace and \gls{ndcg}\xspace, to train the forecaster.
\begin{align}
    \label{eq:joint_loss}
    \mathcal{L}({\mathbf{c}}_i, \hat{{\mathbf{c}}}_i) &= \frac{\doublebars{{\mathbf{c}}_i - \hat{{\mathbf{c}}}_i}^2}{\doublebars{{\mathbf{c}}_i}^2} + 1 - \mathrm{NDCG}({\mathbf{c}}_i, \hat{{\mathbf{c}}}_i).
\end{align}
The first term in \eqref{eq:joint_loss} is the \gls{nmse}\xspace which penalizes any shift from the target values in the aforementioned relative sense, rather than in the absolute sense.
The second term is the \gls{ndcg}\xspace ranking loss \cite{wang13theoretical} which is defined as
\begin{align}
    \label{eq:ndcg_def}
    \mathrm{NDCG}({\mathbf{c}}_i, \hat{{\mathbf{c}}}_i) & \triangleq \frac{\mathrm{DCG}}{\mathrm{IDCG}}~~\textrm{where,}\\
    \mathrm{DCG} &= \sum_m \frac{g({\mathbf{c}}_{\hat{m}}) - 1}{\log_2(m)}\\
    \mathrm{IDCG} &= \sum_m \frac{g({\mathbf{c}}_{[m]}) - 1}{\log_2(m)}.
\end{align}
Here, $g(x) = 2^{\varsigma(x)}$, $\varsigma(\cdot)$ represents the sigmoid function, $(\cdot)_{[m]}$ denotes the $m$-th largest element of a vector and $\hat{m}$ is the index of the $m$-th largest element of the prediction $\hat{{\mathbf{c}}}_i$.
The \gls{ndcg}\xspace metric is used in the ranking literature to quantify the performance of a recommendation system.
In such systems, it is more desirable to recommend most relevant results first where the relevance of each recommendation is measured by a \textit{relevance score}, which is given by ${\mathbf{c}}_i$ for the load balancing task.
Unlike other ranking measures, \gls{ndcg}\xspace allows for a graded relevance score, and thus, is more suited for our use case.
\Gls{dcg} is the weighted sum of the degree of relevance of the ranked items.
The weight is a decreasing function of the rank (position) of the object, and is therefore called discount.
In practice, the logarithmic decade with the rank $m$ has shown promising results.
\gls{ndcg}\xspace normalizes \gls{dcg}\xspace by \gls{idcg}\xspace, which is simply the \gls{dcg}\xspace measure of the best ranking result i.e. the true ranking induced by target ${\mathbf{c}}_i$.

For the load balancing task we can interpret the achievable rate $c_{ij}[n]$ as the relevance score of \gls{bs}\xspace $j$ for user $i$.
The higher the offered rate, the better is the ``rank" of the \gls{bs}\xspace.
Similar to the ranking problem, we want to preserve the ranking among the \glspl{bs}\xspace induced by ${\mathbf{c}}_i$.
Conventionally, \gls{dcg}\xspace is defined without the sigmoid function.
However, in our experiments we observe that reducing the dynamic range of ${\mathbf{c}}_i$ using the sigmoid function helps train a better forecaster.
Both the \gls{nmse}\xspace and \gls{ndcg}\xspace are upper-bounded by $1$ and contribute equally to the loss.

Computing \gls{ndcg}\xspace requires a sorting operation.
Thus, the loss function \eqref{eq:joint_loss} is discontinuous---hence non-differentiable---and not amenable to optimization.
We instead use the smooth surrogate loss function proposed in \cite{pobrotyn2021neuralndcg} which has the following properties.
First, it introduces a temperature parameter to approximate the sorting operation.
Thus, it makes it infinitely differentiable with non-sparse gradients.
Second, it approaches the true \gls{ndcg}\xspace loss as the temperature parameter approaches zero.
In our experiments, we observe that using the smooth surrogate loss greatly improves the quality of the learned forecaster as compared to the indicator function.

\section{Solving the Control Problem using the Frank-Wolfe Method}
\label{sec:frank_wolfe_algo}
Even though the user association \gls{mpc}\xspace in \eqref{opt:util_max_mpc} is a convex program, solving it using a generic convex solver can be slow, especially for large cellular networks.
Thus, we now present a computationally efficient algorithm using the Frank-Wolfe method \cite{jaggi13revisiting}, also known as the conditional gradient method, that is a low-memory, projection-free, first-order optimization method used for non-linear constrained optimizations.
While competing methods such as the projected gradient descent need a projection step into the constraint set, the Frank-Wolfe method instead solves a linear approximation---the first-order Taylor expansion---in the same set, thus always staying in the feasible set.
In our case, we can analytically solve this linear optimization, making the Frank-Wolfe method a good fit for the presented user association formulation.

Consider the optimization problem in \eqref{opt:fw_generic} where $f(\cdot)$ is a differentiable function and $\mathcal{C}$ represents a convex compact set.
\begin{subequations}
\label{opt:fw_generic}
\begin{align}
    \underset{{\mathbf{z}}}{\mathrm{maximize}} &~f({\mathbf{z}}) \\
    \mathrm{subject~to}~&{\mathbf{z}} \in \mathcal{C}.
\end{align}
\end{subequations}
Such constrained optimization problems can be solved iteratively using the Frank-Wolfe method.
In each iteration there are two key steps: the direction-finding step, and the update step.
In each iteration $k$th, we consider the linear approximation of the problem given by the first-order Taylor approximation of $f(\cdot)$ at the current solution ${\mathbf{z}}_k$, and move in the direction the maximizes this linear approximation.
Specifically, if ${\mathbf{s}}_k$ represents the solution of the following linear program
\begin{subequations}
\label{opt:fw_linear_oracle_general}
\begin{align}
    \underset{{\mathbf{s}}}{\mathrm{maximize}}~&{\mathbf{s}}^T\nabla f({\mathbf{z}}_k) \\
   \mathrm{subject~to}~&{\mathbf{s}}\in\mathcal{C},
\end{align}
\end{subequations}
then the update direction for the $k$th iteration is ${\mathbf{d}}_k = {\mathbf{s}}_k - {\mathbf{z}}_k$.
Given the update direction,
\begin{align}
    {\mathbf{z}}_{k+1} & = {\mathbf{z}}_k + \alpha {\mathbf{d}}_k
\end{align}
where $\alpha$ represents the update step-size.
The algorithm terminates if the Frank-Wolfe gap $h_k = {\mathbf{d}}^{T}_{k}\nabla f({\mathbf{z}}_k)$ falls below some tolerance.

For the proposed \gls{mpc}\xspace, let $f(\{x_{ij}[n]\})$ denote the objective in \eqref{opt:util_max_mpc_obj}.
Specifically,
\begin{align}
    \label{eq:fw_obj}
    f(\{x_{ij}[n]\}) &=  \sum_{n=1}^{H}\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}} x_{ij}[n] \log\left(\frac{\bar{c}_{ij}[n]}{\sum_{u\in\mathcal{U}}x_{uj}[n]}\right)  + \frac{\eta}{2}\left(x_{ij}[n] - x_{ij}[n-1]\right)^2,
\end{align}
where $\bar{c}_{ij}[n]$ is as defined in \eqref{opt:util_max_mpc_rate_def} and the partial derivative of $f(\cdot)$ with respect to $x_{ij}[n]$ is given by
\begin{align}
    \label{eq:fw_partial_derivative}
    \frac{\partial f}{\partial x_{ij}[n]} =& \log(\bar{c}_{ij}[n]) - 1 - \log\left(\sum_{u\in\mathcal{U}}x_{uj}[n]\right) \nonumber\\
    & + \eta(x_{ij}[n] - x_{ij}[n-1]) - \eta(x_{ij}[n+1] - x_{ij}[n])\mathbbm{1}(n<H).
\end{align}
Note that we use braces within the parenthesis to emphasize that $f(\cdot)$ is a multivariate function.
Then, at the $k$th iteration, the direction-finding linear program in \eqref{opt:fw_linear_oracle_general} can be written as,
\begin{subequations}
\label{opt:fw_linear_oracle}
\begin{align}
    \{s^{(k)}_{ij}[n]\} =& \underset{\{s_{ij}[n]\}}{\arg\max} \sum_{n=1}^H\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}} s_{ij}[n]\frac{\partial f}{\partial x^{(k)}_{ij}[n]} \\
   \mathrm{subject~to}\nonumber\\
   &\sum_{j\in\mathcal{B}}s_{ij}[n] = 1 \quad\forall i\in\mathcal{U}, n=1, \dots, H\\
    & s_{ij}[n] \in [0, 1] \quad\forall i\in \mathcal{U}, j\in\mathcal{B}, n=1, \dots, H,
\end{align}
\end{subequations}
where $\{s^{(k)}_{ij}[n]\}$ represents the solution that maximizes the first-order Taylor approximation of $f(\cdot)$ at $\{x^{(k)}_{ij}[n]\}$.
The solution to the above constrained linear optimization is straight forward.
\begin{align}
    \label{eq:fw_linear_oracle_sol}
    j^*_{i}[n] &= \underset{j}{\arg\max} \frac{\partial f}{\partial x^{(k)}_{ij}[n]} \\
    s^{(k)}_{ij}[n] &= \mathbbm{1}(j = j^*_{i}[n]).
\end{align}
In other words, user $i$ chooses the \gls{bs}\xspace $j^*$ which maximizes the partial derivative at the $k$-th iteration.
Given $s^{(k)}_{ij}[n]$, the update direction is defined as $d^{(k)}_{ij}[n] = s^{(k)}_{ij}[n] - x^{(k)}_{ij}[n]$ and update step is given by
\begin{align}
    \label{eq:fw_update}
    x^{(k+1)}_{ij}[n] &= x^{(k)}_{ij}[n] + \alpha^{(k)}d^{(k)}_{ij}[n],
\end{align}
where $\alpha^{(k)}$ is the step-size chosen according to the backtracking line-search rule \cite{pedregosa2020linearly}.
The pseudo-code for our Frank-Wolfe implementation is given in \algoref{algo:fw_btls}.
The function $Q(\cdot)$ in line 10 of \algoref{algo:fw_btls} is a quadratic upper bound on the curvature of $f(\cdot)$ and is defined as
\begin{align*}
    Q(\alpha, L) &= f(\{x_{ij}[n]\}) + \alpha h^{(k)} - \frac{\alpha^2L}{2}\sum^H_{n=1}\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}}d^{(k)}_{ij}[n],
\end{align*}
where $h^{(k)}$ (line 6) denotes the Frank-Wolfe gap at iteration $k$ and $L$ is the Lipschitz constant as defined in \cite{pedregosa2020linearly}.

\begin{algorithm}
\caption{Frank-Wolfe algorithm for user association \gls{mpc}\xspace ($\beta=0.9$, $\delta=2.0$)}\label{algo:fw_btls}
\begin{algorithmic}[1]
\Require $\bar{c}_{ij}[n],~x_{ij}[0]$
\Ensure $x_{ij}[n], n=1,2, \dots, H$
\State Randomly initialize $x^{(0)}_{ij}[n]$ such that $\sum_{i}x^{(0)}_{ij}[n] = 1$
\For{$k=1, 2, \dots$}
\State $j^*_i[n] = \arg\max_j~\partial f(x^{(k)})/\partial x_{ij}[n]$
\State $s^{(k)}_{ij}[n] = \mathbbm{1}(j = j^*_i[n])$
\State $d^{(k)}_{ij}[n] = s^{(k)}_{ij}[n] - x^{(k)}_{ij}[n]$
\State $h^{(k)} = \sum_{n=1}^H\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}} d^{(k)}_{ij}[n]\cdot\partial f(x^{(k)})/\partial x_{ij}[n]$ \Comment{Known as the Frank-Wolfe gap.}
\State // BTLS begins
\State $h^{(k)} = \beta h^{(k-1)}$
\State $\alpha^{(k)} = \min\left\{1, h^{(k)}/\left(h^{(k)}\sum_{n=1}^H\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}}\left(d^{(k)}_{ij}[n]\right)^2\right)\right\}$
\For {$f(\{x_{ij}^{(k)}[n] + \alpha^{(k)} d_{ij}^{(k)}[n]\}) < Q(\alpha^{(k)}, h^{(k)})$}
\State $h^{(k)} = \delta h^{(k)}$
\EndFor
\State //BTLS ends
\State $x^{(k+1)}_{ij}[n] = x^{(k)}_{ij}[n] + \alpha^k\cdot d^{(k)}_{ij}[n]$
\EndFor
\State\Return $x^{(k)}_{ij}[n]$
\end{algorithmic}
\end{algorithm}

\section{Formulating Load Balancing as Utility Maximization}
\label{Sec:problem_formulation}

The system model discussed in the previous section accounts for multiple sources of heterogeneity.
Each \gls{bs}\xspace could have a different transmit power $P_j$, the bandwidth $W_j$ can vary across \glspl{bs}\xspace thus affecting the receiver noise and the achievable rate as the pre-logarithmic term, and each \gls{bs}\xspace can operate on a different carrier frequency with vastly different propagation characteristics.
We now formulate a network utility maximization problem with the intention of fairly utilizing the available heterogeneous resources.
The metric of interest in our formulation is the user's \textit{service rate} which is a function of \gls{sinr}\xspace, \textit{load} on the serving \gls{bs}\xspace, and average throughput degradation experienced by the user due to a triggered handover\xspace.
To ground intuition, we first present a genie-aided optimization which requires non-causal measurements.
In the next section, we propose a learning-based user association solution which leverages a forecaster, modeled as a deep neural network, to estimate these non-causal measurements.
In the following, we describe in detail the service rate and subsequently piece together the genie-aided utility maximization.

\subsection{User Service Rate}
\label{subsec:service_rate}
In a cellular network, each \gls{bs}\xspace typically serves multiple users and it must share its available time and frequency resources among these users.
As a result, user's service rate depends on the load on its serving \gls{bs}\xspace and on the deployed resource scheduling algorithm.
Furthermore, while seamless handovers\xspace are necessary to support user mobility---a key feature of cellular networks---they involve relatively complicated signaling procedures and come with a costly overhead \cite{andrews14overview}.
Frequent handovers\xspace can lead to increased connection drop rates, higher transmission delays, and induce ping-pong effects \cite{feriani22multiobjective}.
These can be especially detrimental for heterogeneous deployments which typically have a smaller average cell size.
Hence, we define a user's service rate as a function of all the aforementioned factors.

The achievable data rate between user $i$ and \gls{bs}\xspace $j$ at time $n$ can be expressed as
\begin{align}
    \label{eq:achievable_rate_def}
    c_{ij}[n] &= W_j\log_2(1 + \mathrm{SINR}_{ij}[n]).
\end{align}
Let $\tho{i}[n]$ denote the \gls{hit}\xspace associated with the network's handover\xspace procedure and $T_s$ denote the total slot duration.
Although not explicit from the notation, the \gls{hit}\xspace is a function of the source and target \glspl{bs}\xspace and the quality of channel between the user and the two \glspl{bs}\xspace.
Furthermore, \gls{hit}\xspace also depends on the proprietary handover\xspace procedure implemented by the operator.
We define the service rate for user $i$ from \gls{bs}\xspace $j$ during slot $n$ as
\begin{align}
    \label{eq:service_rate_def}
    r_{ij}[n] &= \frac{c_{ij}[n]}{\ell_j[n]}\left(1 - \frac{\tho{i}[n]}{T_s}\right),
\end{align}
where $\ell_j[n] = \sum_{i\in\mathcal{U}}x_{ij}[n]$ represents the total number of users being served by \gls{bs}\xspace $j$---the load on \gls{bs}\xspace $j$---and $1/\ell_{j}[n]$ represents the fraction of resources allocated to each user by the scheduler at \gls{bs}\xspace $j$.
A round robin resource scheduler would achieve this exactly and a proportional fair scheduler would achieve this on average \cite{vishwanath02opportunistic}.
The multiplicative term $\left(1 - \frac{\tho{i}[n]}{T_s}\right)$ represents the fraction of time remaining for data transmission after handover\xspace related signaling.
$\tho{i}[n]$ is zero if user $i$ does not change association during slot $n$.
Note that one can easily include an \gls{sinr}\xspace gap from the Shannon capacity to account for discrete modulation with finite constellation size.

\subsection{Genie-Aided User Association}
\label{subsec:genie_aided_mpc}
We now formulate the genie-aided network utility maximization problem.
There is a considerable body of literature exploring different applications of various network utility functions such as network-wide resource allocation \cite{yigal07fairness, gupta22system}, network-wide max-min fairness \cite{rasek20joint}, logarithmic utility for balanced load distribution \cite{ye13user}, and $\alpha$-optimal user association \cite{kim12distributed}.
We use the logarithmic utility function which maximizes the product of service rates to achieve a healthy balance between network sum-rate and fairness and naturally achieves load balancing \cite{srikant2014comnets}.
The optimization takes as input the achievable rates $c_{ij}[n]$ and the \glspl{hit} $\tho{i}$, and yields association indicators $x_{ij}[n]$ which maximize the sum of the logarithm of the service rates $r_{ij}[n]$.

A genie-aided network utility maximization can thus be expressed as
\begin{subequations} \label{opt:util_max_genie}
\begin{align}
\underset{\{x_{ij}[n]\}}{\mathrm{maximize}} ~&\sum_{n=1}^{H}\sum_{j\in\mathcal{B}}\sum_{i\in\mathcal{U}} x_{ij}[n] \log\left(\frac{c_{ij}[n]}{\ell_j[n]}\cdot\left(1 - \frac{\tho{i}[n]}{T_s}\right)\right) \label{opt:util_max_genie_obj}\\
\mathrm{subject~to}~\nonumber\\
& \sum_{j\in\mathcal{B}}x_{ij}[n] = 1 \quad \forall i\in \mathcal{U} \label{opt:util_max_genie_ue_cons}\\
& x_{ij}[n]\in\{0,1\} \quad \forall i\in \mathcal{U},\forall j\in\mathcal{B} \label{opt:util_max_genie_indicator_cons}\\
& \ell_j[n] = \sum_{u\in\mathcal{U}}x_{uj}[n] \quad \forall j\in\mathcal{B}. \label{opt:util_max_genie_load_def}
\end{align}
\end{subequations}
The constraint \eqref{opt:util_max_genie_ue_cons} and the binary integer constraint \eqref{opt:util_max_genie_indicator_cons} together ensure that each user $i$ is associated to exactly one \gls{bs}\xspace in time slot $n$.
Constraint \eqref{opt:util_max_genie_load_def} defines the total load on \gls{bs}\xspace $j$ in slot $n$ as the total number of users being served by $j$.
The achievable rates $c_{ij}[n]$ are related to the \gls{sinr}\xspace by Shannon's formula as shown in \eqref{eq:achievable_rate_def}.
Note that the handover\xspace penalty $\left(1 - \frac{\tho{i}[n]}{T_s}\right)$, is a function of $x_{ij}[1], \dots, x_{ij}[n]$.
In order to account for this dependency, our objective is to maximize the logarithmic utility for a time horizon $H$.
Here, $n=1$ represents the current time slot and $n>1$ represents the future.
As a result, the above optimization is genie-aided because it is solved over a future time horizon assuming non-causal knowledge of $c_{ij}[n]$ and $\tho{i}[n]$.
Nonetheless, it serves as an upper bound on the achievable logarithmic utility.

Given the unknown system dynamics induced by user mobility and handover\xspace procedures, and the \textit{cumulative utility maximization} objective, \eqref{opt:util_max_genie} is a good candidate for an \gls{rl}\xspace-based solution.
Although, a powerful technique, model-free \gls{rl}\xspace is known to be sample-inefficient \cite{wu21load, alcaraz22model}, and does not scale well as the  number of users and \glspl{bs}\xspace increase.
As a result, it often converges to highly sub-optimal policies for large-scale networks with \glspl{bs}\xspace and users densities close to realistic deployments.
We instead propose a learning-based \gls{mpc}\xspace approach where a forecaster is trained to model the system dynamics in a supervised fashion, which is relatively more sample-efficient and stable, and has been quite successful in other fields \cite{agarwal2022taskdriven}.
In \secref{sec:results}, we will compare the proposed \gls{mpc}\xspace approach to \gls{rl}\xspace in terms of their sample-efficiency for training and their generalization ability to different user drops.

\section{Simulation Results}
\label{sec:results}

\begin{table}
\caption{Wireless simulation parameters}
\label{tab:sim_param}
\centering
\begin{tabular}{|c|c|}
\hline
Transmit Power           & $46$ dBm                                                                     \\ \hline
Frequency Bands          & [$3750$ MHz, $750$ MHz]                                                        \\ \hline
Bandwidth                & \begin{tabular}[c]{@{}c@{}}$3700$ MHz: $40$ MHz\\ $700$ MHz: $10$ MHz\end{tabular} \\ \hline
Pathloss and shadowing & Urban macrocell model from 3GPP \cite{3GPP2017} \\ \hline
Number of \gls{bs}\xspace                & $19\times 3\times 2 = 114$                                                                    \\ \hline
Number of users                & $1200$                                                              \\ \hline
User distribution               & Uniform or GMM with $10$ cluster centers                                                                \\ \hline
Inter-site separation              & $500$ m                                                                 \\ \hline
Noise                    & $-174$ + $10$ dBm/Hz                                                           \\ \hline
BS Height                & $25$ m                                                                       \\ \hline
User Height                & $1.5$ m                                                                      \\ \hline
\end{tabular}
\end{table}
We now evaluate the forecaster-aided user association strategy in \eqref{opt:util_max_mpc} through Monte-Carlo simulations.
The proposed approach is compared against a max-$C_{\mathrm{Shannon}}$\xspace association scheme, a genie-aided \gls{mpc}\xspace, and a \gls{hit}\xspace agnostic scheme, in terms of user service rates and number of handovers\xspace per time slots.
These baselines are described in detail in \secref{subsec:service_rate_ho}.
In \secref{subsec:mpc_vs_rl}, we compare the proposed forecaster-aided \gls{mpc}\xspace approach against model-free \gls{rl}\xspace from \cite{gupta21load} in terms of training sample-efficiency and generalization capabilities.
Finally, \secref{subsec:fw_vs_cvxpy} compares the execution time of our Frank-Wolfe implementation with a generic convex solver from CVXPY, and the max-$C_{\mathrm{Shannon}}$\xspace scheme.

\subsection{Evaluation Environment}
The wireless environment for the simulations was developed by Meta using their industry-grade system simulator.
We consider a $19$ site hexagonal deployment with an inter-site distance of $500$ m and $1200$ active users.
Each site supports three sectors and two frequency bands.
Therefore, the total number of \glspl{bs}\xspace is $|\mathcal{B}| = 19\times3\times2 = 114$.
For each Monte-Carlo iteration, $1200$ users are dropped either uniformly or according to a Gaussian mixture model (GMM) with $10$ cluster centers and diagonal covariance matrix, to simulate user clustering and inhomogeneous traffic conditions.
For GMM, during each Monte-Carlo iteration, the cluster centers are chosen uniformly, and the standard deviation in each dimension is $150$ m.
Users move with random speed $v$ drawn from a normal distribution with mean $15$ and variance $9$ \cite{3GPP_RRC_2020}, and their trajectory is simulated by a Gauss-Markov process \cite{tabassum19fundamentals}.
The simulator generates spatially consistent channel coefficients using statistical ray-tracing, which is crucial for a reliable mobility simulation.
It uses the urban macrocell pathloss model adopted by 3GPP \cite[page 26]{3GPP2017}, which incorporates the effects of blockage, the multi-slope nature of pathloss, and the elevation difference between the BSs and the UEs.
The \gls{hit}\xspace $\tho{i}$ is modeled as a discrete random variable.
If a handover\xspace is initiated, it is successful with a probability of $0.8$ and suffers a \gls{hit}\xspace of $20$ ms.
Otherwise, with a probability of $0.2$, it results in a link failure and suffers a \gls{hit}\xspace of $90.76$ ms.
The total slot duration is $T_s = 100$ ms.
The link failure probability and interruption times were obtained through a network simulator-3 (ns-3) simulation\footnote{https://github.com/CollinBrady1993/ns-3-dev-git/tree/Handover} and the parameters for wireless simulation are listed in Table \ref{tab:sim_param}.

\subsection{Service Rates and Handovers}
\label{subsec:service_rate_ho}
\begin{figure}[t!]
    \centering
    \subfloat[Uniform user drop]{%
        \includegraphics[width=3in,keepaspectratio]
        {Figures/rate_cdf_uniform_drop_max_cap_tx_pwr_corrected.pdf}

        \label{fig:rate_cdf_uniform}}
    \quad
    \subfloat[Clustered user drop with $10$ centers]{%
        \includegraphics[width=3in,keepaspectratio]
        {Figures/rate_cdf_clustered_drop_10_max_cap_tx_pwr_corrected.pdf}%

        \label{fig:rate_cdf_clustered}}
    \caption{Service rate CDF for the (a) uniform user drop and (b) clustered user drop. The load balancing gain over the max-$C_{\mathrm{Shannon}}$\xspace scheme is evident at the lower percentiles and more profound for an inhomogeneous traffic scenario.}
    \label{fig:rate_cdf}
\end{figure}

We will now compare the performance of the proposed forecaster-aided \gls{mpc}\xspace solution \eqref{opt:util_max_mpc} with the max-$C_{\mathrm{Shannon}}$\xspace user association, a \gls{hit}\xspace agnostic load balancing scheme, and to a genie-aided \gls{mpc}\xspace in terms of the user  service rate and the number of handovers\xspace triggered in a time slot.
Under the max-$C_{\mathrm{Shannon}}$\xspace scheme, the parent \gls{bs}\xspace assigns the user to the best \gls{bs}\xspace in terms of the achievable rate $c_{ij}[n]$.
For the \gls{hit}\xspace agnostic scheme, the central controller solves the utility maximization in \eqref{opt:util_max_mpc} for $H=1$ and $\eta=0$.
In other words, it only uses the current measurements (mobility-unaware) and does not penalize handovers\xspace.
This is similar to the user association strategies proposed in \cite{ye13user, shen14distributed}.
The genie-aided \gls{mpc}\xspace solves for $x_{ij}[n]$ using \eqref{opt:util_max_mpc} but assumes non-causal knowledge of future rates.

\figref{fig:rate_cdf} shows the service rate \gls{cdf}\xspace for the four user association schemes under different user drop scenarios where
\figref{fig:rate_cdf_uniform} presents uniform user drop and \figref{fig:rate_cdf_clustered} presents clustered user drop which emulates inhomogeneous traffic conditions.
The genie \gls{mpc}\xspace and the forecaster \gls{mpc}\xspace are evaluated for a forecasting horizon of $H=3$.
Compared to the max-$C_{\mathrm{Shannon}}$\xspace scheme, the log-utility-based association schemes provide a more uniform user experience by judiciously distributing users across \glspl{bs}\xspace.
Under the log-utility, the users with a strong wireless channel---for example, those close to the cellular tower---are offloaded to a \gls{bs}\xspace or frequency band with a relatively poor rate, leaving the lower frequency band for the users at the cell edge so that they can benefit from the lower pathloss and better penetration.
In \figref{fig:rate_cdf}, this results in significantly higher service rates at the 5th percentile compared to max-$C_{\mathrm{Shannon}}$\xspace, with a slight degradation beyond the 95th percentile.
We will refer to the improvement in service rate compared to max-$C_{\mathrm{Shannon}}$\xspace as \textit{rate gain}.
In \figref{fig:rate_cdf_uniform} and \figref{fig:rate_cdf_clustered}, we observe that compared to the uniform drop scenario, the 5th percentile rate gain is higher for the clustered user drop where the user traffic is inhomogeneous.
Due to this additional spatial heterogeneity, the user service rates suffer greater degradation under the load-agnostic max-$C_{\mathrm{Shannon}}$\xspace association strategy where \glspl{bs}\xspace close to the user clusters end up serving most of the users in the network while others remain under-utilized.
As a result, load-aware association is more beneficial for the clustered user drop where the 5th percentile rate gain is about $3.5\times$ as opposed to $1.5\times$ for the uniform user drop.
\begin{figure}[t!]
    \centering
    \subfloat[Uniform user drop]{%
        \includegraphics[width=3in,keepaspectratio]
        {Figures/percentile_rate_uniform_drop_max_cap_tx_pwr_corrected.pdf}

        \label{fig:rate_percentile_uniform}}
    \quad
    \subfloat[Clustered user drop with $10$ centers]{%
        \includegraphics[width=3in,keepaspectratio]

        {Figures/percentile_rate_gain_clustered_drop_10.pdf}%

        \label{fig:rate_percentile_clustered}}

    \caption{Rate gain from load balancing across different percentiles for (a) uniform user drop and (b) clustered user drop. The forecaster \gls{mpc}\xspace achieves service rates similar to genie \gls{mpc}\xspace. The \gls{hit}\xspace agnostic approach suffers throughput degradation from handovers\xspace.}
    \label{fig:rate_percentile}
\end{figure}

The coverage improvement, as measured by the 5th percentile service rate, can be more clearly observed in \figref{fig:rate_percentile_uniform} and \figref{fig:rate_percentile_clustered} which present the rate gain from the three log-utility-based association schemes at different percentiles for both uniform and clustered user drops.
However, the gain diminishes at the higher percentiles because, compared to the max-$C_{\mathrm{Shannon}}$\xspace scheme, under the log-utility `good' users suffer a slight degradation in \gls{sinr}\xspace to help improve the overall network performance.
Furthermore, \figref{fig:rate_percentile} shows that the forecaster \gls{mpc}\xspace achieves service rates similar to the genie \gls{mpc}\xspace, especially for the uniform user drop scenario.
Among the three load-aware association strategies, the mobility-aware \gls{mpc}\xspace approaches benefit from knowledge of future user trajectories and thus, can avoid handover\xspace related degradation to users' service rates.
On the other hand, the \gls{hit}\xspace agnostic scheme is plagued with excessive handovers\xspace and as a result the edge user performance suffers.
In \figref{fig:rate_percentile_clustered}, the \gls{hit}\xspace agnostic scheme achieves a 5th percentile rate gain of about $2.6\times$ over max-$C_{\mathrm{Shannon}}$\xspace.
This means, compared to the proposed forecaster-aided \gls{mpc}\xspace, \gls{hit}\xspace agnostic scheme suffers a $39\%$ loss in edge user service rate due to additional handovers\xspace.

\begin{figure}[t!]
    \centering
    \subfloat[Uniform user drop]{%
        \includegraphics[width=3in,keepaspectratio]
        {Figures/ho_cdf_uniform_drop_max_cap_tx_pwr_corrected.pdf}

        \label{fig:ho_cdf_uniform}}
    \quad
    \subfloat[Clustered user drop with $10$ centers]{%
        \includegraphics[width=3in,keepaspectratio]
        {Figures/ho_cdf_clustered_drop_10_max_cap_tx_pwr_corrected.pdf}

        \label{fig:ho_cdf_clustered}}
    \caption{\gls{cdf}\xspace of handovers\xspace triggered in a time slot for (a) uniform user drop and (b) clustered user drop. Forecaster \gls{mpc}\xspace significantly reduces handovers\xspace compared to \gls{hit}\xspace agnostic by leveraging the available forecast of user rates.}
    \label{fig:ho_cdf}
\end{figure}
The handover\xspace performance of the four association schemes is compared in \figref{fig:ho_cdf} where
\figref{fig:ho_cdf_uniform} and \figref{fig:ho_cdf_clustered} present the \gls{cdf}\xspace of the number of handovers\xspace initiated in a time slot.
The \gls{hit}\xspace agnostic scheme initiates the most number of handovers\xspace, as evident from the \gls{cdf}\xspace, and the proposed forecaster-aided \gls{mpc}\xspace scheme reduces the median number of handovers\xspace by $7\times$ for the clustered user drop and by $5\times$ for the uniform user drop.
Compared to the genie-aided scheme, the median of handover\xspace \gls{cdf}\xspace for the proposed scheme is about $1.5\times$ greater.
Note, although not considered in this work, reducing handovers\xspace has benefits beyond improving user service rate.
For example, unnecessary handovers\xspace and the associated signaling overhead can lead to increased session latency.
Similarly, reducing unwanted handovers\xspace will reduce energy consumption which is especially beneficial for mobile devices.
These results illustrate that using the proposed forecaster design with the novel loss function, the \gls{mpc}\xspace in \eqref{opt:util_max_mpc} can achieve performance close to the genie-aided scheme in terms of both the user  service rate and the number of handovers\xspace.

\subsection{\gls{mpc}\xspace vs \gls{rl}\xspace: Sample-efficiency \& Generalization}
\label{subsec:mpc_vs_rl}
In this section, we will compare the proposed \gls{mpc}\xspace approach to model-free \gls{rl}\xspace in terms of the training sample-efficiency and their ability to generalize to different user drops.
We use the \gls{mdp}\xspace model and the neural network architecture proposed in \cite{gupta21load}.
However, we consider a simulation environment with $10\times$ more users and \glspl{bs}\xspace.
As such, unlike our previous work \cite{gupta21load} we had little success with the deep recurrent Q-network (DRQN) algorithm and instead use the proximal policy optimization (PPO) algorithm \cite{schulman17proximal} for training.
The \gls{mdp}\xspace corresponding to mobility-aware user association  can be defined as the 5-tuple $\langle\mathcal{S},\mathcal{X},\mathcal{T},\mathcal{R},\gamma\rangle$.
$\mathcal{S}$ is the state space and $\mathcal{X}$ represents the action space.
At each time step $n$, the central controller receives a representation of the network's current \textit{state}, $S[n] \in \mathcal{S}$, and on that basis executes an \textit{action}, ${\mathbf{X}}[n] \in \mathcal{X}$, causing the environment state to transition to a new state $S[n+1]$ with probability ${\mathbb{P}}(S[n+1]|S[n], {\mathbf{X}}[n]) = \mathcal{T}(S[n],{\mathbf{X}}[n],S[n+1])$.
In the next time step the central controller receives the \textit{reward} $R[n+1] = \mathcal{R}(S[n], {\mathbf{X}}[n], S[n+1])$ and the new state $S[n+1]$.
Over an episode of $T$ time slots, the controller collects a cumulative reward $\sum_{n=1}^T\gamma^nR[n+1]$ where $\gamma \in (0, 1]$
is known as the discount factor introduced mainly to have a
finite cumulative reward, especially as $T$ tends to infinity.
The objective of the controller is to maximize the expected cumulative reward.
The specific formulations of the state, action, and reward are given below and summarized in \figref{fig:mlb_mdp}.

\begin{figure}[t!]
    \centering
    \includegraphics[width=6in,                 keepaspectratio]{Figures/MLB-MDP.pdf}
    \caption{\gls{mdp}\xspace for mobility-aware user association.}
    \label{fig:mlb_mdp}
\end{figure}

\textbf{State:} The central controller chooses the new association based on the $|\mathcal{U}|\times|\mathcal{B}|$ \gls{sinr}\xspace measurements $\bm{\Phi}[n]$ received for each user-\gls{bs}\xspace pair and the $|\mathcal{B}|\times 1$ BS load $\bm{\ell}[n]$.
Thus, the state at time $n$ is defined as a tuple
\begin{align}
    \label{eq:state_def}
    S[n] = (\bm{\Phi}[n], \bm{\ell}[n]).
\end{align}

\textbf{Action:} Observing the \gls{sinr}\xspace for user-\gls{bs}\xspace pairs and the load on each \gls{bs}\xspace, the central controller returns the new association indicators $\{x_{ij}[n]|i\in\mathcal{U}, j\in\mathcal{B}\}$ such that they satisfy \eqref{opt:util_max_genie_ue_cons} and \eqref{opt:util_max_genie_indicator_cons}.
The action at time $n$ is thus denoted by a $|\mathcal{U}|\times|\mathcal{B}|$ binary matrix ${\mathbf{X}}[n]$ such that the sum of each row is $1$.

\textbf{Reward:} Given the association matrix ${\mathbf{X}}[n]$, the serving BS for user $i$ is given by $j^* = \arg\max_jx_{ij}[n]$ and it receives utility $\log(r_{ij^*}[n])$ where $r_{ij^*}[n]$ is given by \eqref{eq:service_rate_def}.
The \gls{hit}\xspace $T_{{\rm HO},i}[n]$ is non-zero if user $i$ changes serving \glspl{bs}\xspace, $x_{ij^*}[n] \neq x_{ij^*}[n-1]$.
A handover\xspace can either be successful or can result in an link failure.
If an link failure event occurs, the user re-establishes its connection with the network via the radio resource control (RRC) protocol, which results in a large interruption time of the order of hundreds of milliseconds.
The reward received by the controller at time $n+1$ is the sum of utility received by all users.
\begin{align}
\label{eq:reward_def}
    R[n+1] &= \sum_{j\in\mathcal{B}}\sum_{i\in \mathcal{U}}x_{ij}[n] \log\left(r_{ij}[n]\right)
\end{align}
Note that $\sum_{n=1}^TR[n+1]$ is equal to the objective in \eqref{opt:util_max_genie_obj}.

\textbf{State transition probabilities:} The state $S[n]$ has two parts.
The evolution of $\bm{\Phi}[n]$ is governed by the motion of users and is independent of actions taken by the controller.
However, the load on each BS depends on the past action ${\mathbf{X}}[n-1]$ and the results of the handover\xspace procedure.
Neither user mobility nor the handover\xspace procedure are explicitly modeled and the \gls{rl}\xspace agent must learn the underlying system dynamics and determine the optimal policy in a model-free manner.
\begin{figure}
    \centering
    \includegraphics[width=4in,keepaspectratio]{Figures/forecaster_genie_rate_gain_vs_num_samples.pdf}%

    \caption{Comparing training sample efficiency of user rate forecaster and \gls{rl}\xspace agent using the PPO algorithm. Even with $100\times$ more training samples, the PPO agent achieves 5th percentile service rate similar to max-$C_{\mathrm{Shannon}}$\xspace.}
    \label{fig:sample_efficiency}
\end{figure}

For this \gls{mdp}\xspace model, an \gls{rl}\xspace agent is trained for centralized control using the PPO algorithm. \figref{fig:sample_efficiency} presents the growth of 5th percentile rate gain of model-free \gls{rl}\xspace and the proposed forecaster \gls{mpc}\xspace with the number of training samples.
We train the PPO agent and the forecaster the same dataset.
A training sample for the PPO agent refers to the $4$-tuple $(S[n], {\mathbf{X}}[n], R[n+1], S[n+1])$, as is typical for \gls{rl}\xspace.
For the forecaster, a training sample refers to a labeled sample for supervised training as defined in \secref{sec:forecaster_design}, where the input feature is the time-series of current and past capacity measurements and the target label is the measurement for future time steps.
Note, that one \gls{rl}\xspace training sample includes measurements from all users in the network.
However, the forecaster predicts the trajectory for each user independently and thus, treats measurements from each user as a separate training sample.
We observe in \figref{fig:sample_efficiency} that the forecaster needs about $10^4$ training samples to achieve performance close to the genie.
However, even with $100\times$ more training samples, the PPO agent can only achieve a edge user service rate similar to max-$C_{\mathrm{Shannon}}$\xspace, i.e. a 5th percentile rate gain of $1$.
In terms of training time, training the PPO agent with $10^6$ samples took about 3 weeks, where the forecaster can be trained in a few hours.
Both these models were trained on a machine with an Intel(R) Xeon(R) Gold 6226R CPU @ 2.90GHz and NVIDIA RTX A5000 GPU.
Note that in \figref{fig:sample_efficiency}, even an untrained forecaster can improve the edge user service rate by $2.7\times$.
This is because it has access to the true capacity measurements for the current time step.
Thus, even with random predictions from an untrained forecaster and reliable measurements for the current slot, the \gls{mpc}\xspace scheme leads to a well-balanced user distribution across bands and \glspl{bs}\xspace.
This is another advantage of our model-based approach over model-free \gls{rl}\xspace.
The gap between untrained forecaster-aided \gls{mpc}\xspace and the genie-aided \gls{mpc}\xspace represents the forecasting gain at $H=3$.

\begin{figure}
    \centering
    \subfloat[Uniform drop]{%
        \includegraphics[width=3.15in,
                        keepaspectratio]
        {Figures/rate_cdf_uniform_drop_ppo_max_cap_tx_pwr_corrected.pdf}%

        \label{fig:ppo_rate_cdf_uniform}}
    \quad
    \subfloat[Clustered drop]{%
        \includegraphics[width=3.15in,
                        keepaspectratio]
        {Figures/rate_cdf_clustered_drop_10_ppo_max_cap_tx_pwr_corrected.pdf}%

        \label{fig:ppo_rate_cdf_clustered}}

    \caption{Comparing the generalization capability of \gls{rl}\xspace and the rate forecaster. The PPO agent and the forecaster are trained on (a) uniform user drop scenario. These trained models are then used for inference on (b) a clustered user drop. PPO agent exhibits poor generalization while the forecaster reliably learns the user mobility patterns and generalizes across user drops.}
    \label{fig:ppo_rate_cdf}
\end{figure}

Next, we train the forecaster using $10^4$ training samples and the PPO agent using $10^6$ training samples from the uniform user drop scenario.
\figref{fig:ppo_rate_cdf_uniform} shows the \gls{cdf}\xspace of the service rates achieved by the two user association strategies on the uniform user drop.
We observe that the forecaster-aided \gls{mpc}\xspace achieves performance close to the genie and the PPO agent achieves performance similar to the max-$C_{\mathrm{Shannon}}$\xspace strategy.
These trained models are then tested on clustered user drops and the corresponding service rate \glspl{cdf} are presented in \figref{fig:ppo_rate_cdf_clustered}.
While the forecaster \gls{mpc}\xspace can effectively translate its learned mobility model to the clustered drop scenarios, the PPO agent suffers significant degradation and exhibits poor generalization.
This is probably because the PPO agent learns some specific attribute of the uniform drop dataset which does not translate to the clustered drop scenario.
Such uninterpretability and poor generalization of \gls{rl}\xspace methods have been well documented and methods to combat such behavior are an active area of research.

\subsection{Computation Efficiency of Frank-Wolfe}
\label{subsec:fw_vs_cvxpy}
In the next section we evaluate the computational efficiency our optimization algorithm which is based on the Frank-Wolfe method.
We compare our Frank-Wolfe-based implementation to a splitting conic solver (SCS) implementation from CVXPY \cite{cvxpy1}, a standard software tool for convex optimization.
In \tabref{tab:obj_value} we present the user service rates, the cell sum-rates, and the log-utility achieved by the two optimization algorithms for different forecasting horizons.
The cell sum-rate for a \gls{bs}\xspace $j$ is defined as $\sum_i x_{ij}[n]r_{ij}[n]$.
The two algorithms achieve similar performance on all three metrics and this validates that our algorithm converges to the correct global optima.

In \tabref{tab:compute_time}, we present the per time slot execution time for solving the proposed \gls{mpc}\xspace in \eqref{opt:util_max_mpc} for $H=1$ and $H=2$.
The execution times in \tabref{tab:compute_time} are normalized by the execution time of the max-$C_{\mathrm{Shannon}}$\xspace scheme.
Note, the tables present execution time for solving the \gls{mpc}\xspace which does not include the interference time from the forecaster.
We observe that on our fairly large network with $114$ \glspl{bs}\xspace and $1200$ users, the SCS solver is quite computationally expensive.
However, our Frank-Wolfe implementation reduces the execution time by two orders of magnitude for a time horizon of $H=1$ and $H=2$.
This is because the direction-finding step in the Frank-Wolfe method \eqref{opt:fw_linear_oracle} poses a closed-form solution which was discussed in \secref{sec:frank_wolfe_algo}.
As such, the solution to the linear program reduces to a simple maximization over an array.
Using the backtracking line-search procedure proposed in \cite{pedregosa2020linearly} also aides in faster convergence by reducing the number of iterative steps.
These effects together make the proposed \gls{mpc}\xspace solution practically viable.

\begin{table}
\caption{Average service rate and objective achieved by splitting conic solver (SCS) from CVXPY and our optimization algorithm based on the Frank-Wolfe method.}
\label{tab:obj_value}
\centering
\begin{tabular}{|c|ccc|c|c|}
\hline
\multirow{2}{*}{} & \multicolumn{3}{c|}{User service rate (Mbps)}                                            & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Mean cell\\ sum-rate\end{tabular}}        & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Mean\\ objective value\end{tabular}} \\ \cline{2-4}
                  & \multicolumn{1}{c|}{5th percentile} & \multicolumn{1}{c|}{50th percentile} & 90th percentile &            &                                                                       \\ \hline
CVXPY + SCS       & \multicolumn{1}{c|}{0.69}        & \multicolumn{1}{c|}{2.06}        & 4.54       & 25.98 & 17410.72                                                                        \\ \hline
Our method      & \multicolumn{1}{c|}{0.70}        & \multicolumn{1}{c|}{2.04}        & 4.66       & 26.09 & 17411.60                                                                       \\ \hline
\end{tabular}
\end{table}

\begin{table}
\caption{Average per time slot execution time relative to the max-$C_{\mathrm{Shannon}}$\xspace scheme}
\label{tab:compute_time}
\centering
\begin{tabular}{|c|c|c|}
\hline
                      & \begin{tabular}[c]{@{}c@{}}Forecaster-aided  \\ \gls{mpc}\xspace for $H=1$\end{tabular} & \begin{tabular}[c]{@{}c@{}}Forecaster-aided  \\ \gls{mpc}\xspace for $H=2$\end{tabular}\\ \hline
CVXPY + SCS           & $137.5$          & $187.5$          \\ \hline
Our Method & $1.25$           & $3.26$           \\ \hline
\end{tabular}
\end{table}
\section{System Model}
\label{Sec:system_model}
We consider a downlink cellular network which consists of $K$ different carrier frequency bands.
Let $\mathcal{B}_k, \text{ where } k\in\{1, 2, \dots, K\}$, denote the set of all \glspl{bs}\xspace operating on the $k$-th frequency band and $\mathcal{U}$ denote the set of all users in the network.
We characterize a \gls{bs}\xspace by its physical location \textit{and} its carrier frequency band.
Thus, in this formulation, radio units operating on different carrier frequencies but mounted on the same physical tower are treated as different co-located \glspl{bs}\xspace, as shown in \figref{fig:network_struct}.
Note that $\mathcal{B}_k$'s are disjoint sets and the set of all \glspl{bs}\xspace is denoted as $\mathcal{B} = \bigcup_{k=1}^K\mathcal{B}_k$.
As described later in \secref{sec:forecaster_design}, the set $\mathcal{U}$ can vary over time.

\begin{figure}[t!]
    \centering
    \includegraphics[width=4in,                 keepaspectratio]{Figures/network_struct.pdf}
    \caption{Network structure and flow of \gls{sinr}\xspace measurements. The \glspl{bs}\xspace receive \gls{sinr}\xspace measurements from the user in their cell. The central controller receives \gls{sinr}\xspace measurements from the \glspl{bs}\xspace and communicates the association indicators back to the \glspl{bs}\xspace and users.}
    \label{fig:network_struct}
\end{figure}

We will assume slotted time indexed by $n$.
The users are mobile and each is initially associated with a parent \gls{bs}\xspace.
It is assumed that user (re-)association decisions are made on a larger timescale compared to the scheduling of time and frequency resources.
Consequently, many scheduling slots are subsumed under a user association slot \cite{ye13user, shen14distributed}.
With this distinction, henceforth, a time slot will refer to a user association slot unless stated otherwise.
The \gls{sinr}\xspace used for association is averaged over the association interval, and thus, the small-scale fading is also averaged out.
The \gls{sinr}\xspace of the signal delivered to user $i$ by \gls{bs}\xspace $j\in\mathcal{B}_k$ during time slot $n$ can be expressed as
\begin{align}
    \mathrm{SINR}_{ij}[n] &= \frac{P_j g_{ij}[n]}{\sum_{j'\in\mathcal{B}_k, j'\neq j}P_{j'}g_{ij'}[n] + W_j\sigma^2},
\end{align}
where $P_j$ denotes the transmit power at \gls{bs}\xspace $j$, $g_{ij}[n]$ represents the channel gain between user $i$ and \gls{bs}\xspace $j$ and includes the pathloss, shadowing, and the antenna gains, and $\sigma^2$ denotes the power spectral density of the thermal noise.
Since the \gls{bs}\xspace index $j$ also identifies its frequency band and different frequency bands can operate on different bandwidths, we denote by $W_j$ the bandwidth at \gls{bs}\xspace $j$.

At the beginning of a time slot $n$, each user measures the \gls{sinr}\xspace from its neighboring \glspl{bs}\xspace and reports these measurements to the parent \gls{bs}\xspace.
The \glspl{bs}\xspace then forward these \gls{sinr}\xspace measurements to a central network controller, which processes these measurements and computes the association indicators $x_{ij}[n]$, which are defined as
\begin{align}
    x_{ij}[n] &= \begin{dcases}
                    1, & \textrm{user $i$ associates to \gls{bs}\xspace $j$ in slot $n$} \\
                    0, & \textrm{otherwise}.
                 \end{dcases}
\end{align}
These new association decisions are communicated to the \glspl{bs}\xspace by the controller and then forwarded to the users by their current parent \gls{bs}\xspace.
The remaining time in slot $n$ is then used for handover\xspace and data transmission.
A handover\xspace is triggered for user $i$ if $x_{ij}[n] \neq x_{ij}[n-1]$.
Note that in our formulation, switching frequency bands on the same site also counts as a handover\xspace.
\figref{fig:network_struct} summarizes the signal flow in a time slot.
Such closed-loop control is a major thrust for open radio access networks (O-RAN) where the RAN intelligent controller (RIC) can centrally control hundreds of \glspl{bs}\xspace while efficiently handling feedback and real-time centralized optimization \cite{lin22embracing, lacava22programmable}.
Furthermore, the central availability of data and analytics in the network enables using intelligent data-driven approaches to user association and generally network resource management by taking a user-centric approach \cite{polese22understanding}.

\section{System Model}
\label{Sec:SysMod}

A downlink cellular network consisting of $K$ different frequency bands is considered.
Let $\mathcal{B}_k,~ k\in \{1, 2, \cdots, K\}$, denote the set of all base stations (BS) operating in the $k$-th frequency band and $\mathcal{U}$ be the set of all UEs in the network.
In this formulation, a BS is characterized by its physical deployment location \textit{and} the carrier frequency band used for transmissions.
In other words, a physical site that supports two distinct frequency bands counts as two BSs, Fig. \ref{fig:network_struct}.
It is easy to see that $\mathcal{B}_k$ are disjoint sets and the set of all BSs can be denoted as $\mathcal{B} = \bigcup_{k = 1}^K\mathcal{B}_k$.

The UEs are mobile and measure the received signal strength from their neighboring BSs, which can span multiple frequency bands, by measuring {\rm RSRP}s.
\begin{figure}[t!]
    \centering
    \includegraphics[width=4in,                 keepaspectratio]{Figures/network_struct.pdf}
    \caption{Network structure and flow of RSRP and SINR measurements. The BSs receive RSRP measurement from the UEs in their cell. The central controller receives SINR measurements and load information from the BSs and communicates the association indicators back to the BSs and UEs.}
    \label{fig:network_struct}
\end{figure}
A slotted time model is considered where each UE reports the {\rm RSRP} from top-$m$ BSs to its serving BS, at the beginning of time slot $t$. For each UE $i$, the serving BS constructs the set $\{\mathrm{RSRP}_{ij}[t]| j\in\mathcal{B}\}$ by filling in zeros for the BSs for whom the measurements were not reported.
Hence, $m$ should be large enough to include significant interferers. Given the RSRP measurements for a UE $i$, its SINR with BS $j\in \mathcal{B}_k$ can be written as
\begin{align}
\label{eq:sinr_definition}
    \mathrm{SINR}_{ij}[t] &= \frac{\mathrm{SNR}{ij}[t]}{\sigma^2W_j + \sum_{j'\in \mathcal{B}_k, j'\neq j}\mathrm{SNR}{ij'}[t]},
\end{align}
where $W_j$ represents the bandwidth of BS $j$ and $\sigma^2$ is the power spectral density of additive noise.
Note that SINR is different from RSRP because it accounts for interference and depends on the bandwidth, whereas RSRP is a function of transmit power, pathloss, small scale fading, and large scale fading \cite[Section 8.1]{3GPP_36_873}.
Let $\bm{\Phi}[t]$ denote the $|\mathcal{U}|\times|\mathcal{B}|$ matrix of SINRs such that $[\bm{\Phi}[t]]_{i,j} = \mathrm{SINR}_{ij}[t]$.
$\ell_j[t]$ denote the load on BS $j$ (defined shortly), and $\bm{\ell}[t]$ be the vector of length $|\mathcal{B}|$ representing the load on each BS in $\mathcal{B}$.
The BSs communicate their load and SINRs from their UEs to a central network controller in the RIC platform.
The controller then decides the UE-BS association for the next time slot and communicates it to all the BSs, which then communicate it to the UEs.
A HO event occurs if a UE changes its serving BS from one time slot to the next.
Note that in this formulation HO occurs whenever the UE switches its band or site.

Since each BS generally serves more than one UE, its time and frequency resources need to be shared and the data rate for a UE depends on the load of the serving BS. The channel is assumed to be static during a time slot and all BSs are assumed to have full-buffers. The resources at each BS are distributed equally among the associated UEs, which is optimal in terms of log-utility achieved at each BS with static channels \cite{Ye13}. The data rate received by a UE $i$ associated to BS $j\in\mathcal{B}_k$ can be expressed as
\begin{align}
    \label{eq:loaded_capacity_definition}
    c_{ij}[t] &= \frac{W_j}{\ell_j[t]}\log_2(1 + \mathrm{SINR}_{ij}[t]).
\end{align}
In case an HO is initiated, the UE's data transmission is interrupted for the related control signal exchange between the source and the target BSs \cite[Section 6]{3GPP_38_133}. Let $\tho{i}[t]$ be the HO interruption time experienced by UE $i$ during time slot $t$ and $T_s$ be the slot duration. The throughput for UE $i$ associated to BS $j\in\mathcal{B}_k$, during time slot $t$, can be expressed as
\begin{align}
    \label{eq:HO_capacity_definition}
    r_{ij}[t] &= \left(1 - \frac{\tho{i}[t]}{T_s}\right)c_{ij}[t]\\
                &= \left(1 - \frac{\tho{i}[t]}{T_s}\right)\frac{W_j}{\ell_j[t]}\log_2(1 + \mathrm{SINR}_{ij}[t]).
\end{align}
The following sections mathematically formulate the long-term utility maximization problem and the MDP for MLB.