## 3 Proposed Optimization Algorithm
### BS Transmit Power Allocation
#### 3.1.2 Two-Layer Approach for Solving Problem (P1)
This bilinear term, which involves the multiplication of a binary variable and a continuous variable, necessitates the introduction of additional auxiliary variables and constraints for accurate linearization. Linearization: To linearize the bilinear term \(\gamma_{m,j}z_{j}\), we introduce a set of auxiliary variables \(b_{m,j}\). Subsequently, the problem (40) can be reformulated as an equivalent expression:

\[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]2 (43) s.t. (19)-(21) \[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]3 \[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]4 (46) \[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]5

To address the bilinear constraint (46) in problem (43), we can utilize the binary nature of \(\gamma_{m,j}\) to transform it. Specifically, the constraint \(b_{m,j}=\gamma_{m,j}z_{j}\) can be equivalently represented as the following set of linear constraints:

\[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]6 \[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]7 (49) \[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]8

Here, \(K>0\) is a large enough constant. Constraint (49) implies that, if \(\gamma_{m,j}=0\), then \(b_{m,j}=0\); constraints (48) and (50) enforce that, if \(\gamma_{m,j}=1\), then \(b_{m,j}=z_{j}\). Therefore, constraints (48)-(50) serve as linearization constraints for \(b_{m,j}=\gamma_{m,j}z_{j}\) and the equivalence holds here as long as \(K\geq z_{j\in\mathcal{N}}\left\{z_{j}\right\}\). By linearizing (46) into (48)-(50), we can finally reformulate problem (P2) into a MILP format that can be solved efficiently as follows:\[v_{2}=\sum_{n\in\mathcal{N}}\mu_{n}\rho_{n}\left(\sum_{n\in \mathcal{N}}\frac{\sum\limits_{m\in\mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{ \sum\limits_{m\in\mathcal{M}}\gamma_{m,n}d_{m,n}}-\frac{\sum\limits_{m\in \mathcal{M}}\gamma_{m,n}u_{m,n}d_{m,n}}{\sum\limits_{m\in\mathcal{M}}\gamma_{ m,n}d_{m,n}}\right),\]9 (51) s.t.