无线通信MISO系统的多用户的发射功率最小化问题（附Matlab代码，使用SOCP、SDP以及KKT）_matlab求解socp问题

引言

前面几篇博客介绍并复现RIS辅助无线通信系统中的相关优化问题，为了更深度地理解无线通信中的优化问题，打算利用三篇博客介绍并仿真发射功率最小化问题、和速率最大化问题与能效问题等三大经典问题，本文首先介绍多用户MISO系统基站发射功率最小化问题。博客和代码都为自己原创，希望大家多点赞收藏，支持创作。

系统模型

考虑一个下行多用户多输入单输出（MU-MISO）系统，其中基站配置有$N$根发射天线，以及$K$个单天线用户。该系统的发射功率最小化问题可以表示为：
$$\begin{gathered} \underset{{\mathbf{w}}_1,\dots ,{\mathbf{w}}_K}{\min }\sum _{k=1}^K{\Vert {\mathbf{w}}_k\Vert }^2 \\ \text{s.t.}\frac{{\left|{\mathbf{h}}_k^H{\mathbf{w}}_k\right|}^2}{{\sum }_{j\ne k}{\left|{\mathbf{h}}_k^H{\mathbf{w}}_j\right|}^2+{\sigma }_k^2}⩾{\gamma }_k,k=1,\dots K. \end{gathered}$$
其中$p_k$表示分配给用户$k$的发射功率。

求解方法

在上述优化问题的基础上，文本介绍三种无线通信里最常用的求解算法：

1. SOCP(二阶锥规划)
2. SDP(半正定规划)
3. KKT(Karush-Kuhn-Tucker最优化条件)

SOCP

首先介绍SOCP，对于SINR约束条件，不等式两边同时乘以分子，并除以SINR阈值，可转化为一个二阶锥约束：
$${\left|{\mathbf{h}}_k^H{\mathbf{w}}_k\right|}^2+\frac{1}{{\gamma }_k}{\left|{\mathbf{h}}_k^H{\mathbf{w}}_k\right|}^2⩾\sum _{j=1}^K{\left|{\mathbf{h}}_k^H{\mathbf{w}}_j\right|}^2+{\sigma }^2$$
令$\mathbf{W}=\left[{\mathbf{w}}_1,\dots ,{\mathbf{w}}_K\right]$，将右边简化为范数的形式：
$$\sqrt{1+\frac{1}{{\gamma }_k}}\left|{\mathbf{h}}_k^H{\mathbf{w}}_k\right|⩾\Vert {\mathbf{h}}_k^H\mathbf{W}\sigma \Vert$$
在此使用一个常用的小技巧[1]，由于${\mathbf{w}}_k$和进行一定相位旋转后$e^{j{\theta }_k}{\mathbf{w}}_k$（${\theta }_k\in \mathbb{R}$）后的绝对值是相等的，因此考虑将内积${\mathbf{h}}_k^H{\mathbf{w}}_k$做一定的相位旋转，使得内积为正实数，即：$\left|{\mathbf{h}}_k^H{\mathbf{w}}_k\right|={\mathbf{h}}_k^H{\mathbf{w}}_k\ge 0$，因此上述约束可以进一步转化为：
$$\begin{gathered} \sqrt{1+\frac{1}{{\gamma }_k}}\Re ({\mathbf{h}}_k^H{\mathbf{w}}_k)⩾\Vert {\mathbf{h}}_k^H\mathbf{W}\sigma \Vert \\ \Im ({\mathbf{h}}_k^H{\mathbf{w}}_k)=0 \end{gathered}$$

为简化表达，将目标函数可以转化为：
$$\sum _{k=1}^K{\Vert {\mathbf{w}}_k\Vert }^2=\text{trace}\left({\mathbf{W}}^H\mathbf{W}\right)=\text{vec}{\left(\mathbf{W}\right)}^H\text{vec}\left(\mathbf{W}\right)={\Vert \text{vec}\left(\mathbf{W}\right)\Vert }^2$$
因此原问题可以等价转化为以下优化问题：
$$\begin{gathered} \underset{\mathbf{W}}{\min }\Vert \text{vec}\left(\mathbf{W}\right)\Vert \\ \text{s.t.}\sqrt{1+\frac{1}{{\gamma }_k}}\Re ({\mathbf{h}}_k^H{\mathbf{w}}_k)⩾\Vert {\mathbf{h}}_k^H\mathbf{W}\sigma \Vert , \\ \Im ({\mathbf{h}}_k^H{\mathbf{w}}_k)=0,k=1,\dots K \end{gathered}$$
目标函数是凸函数，第一个约束条件是二阶锥规划，是凸约束，第二个约束条件是线性约束，因此可以使用CVX进行求解。

SDP

同样，原问题可以简单地转化为半正定规划问题，令${\mathbf{B}}_k={\mathbf{w}}_k{\mathbf{w}}_k^H$为半正定的Hermitian矩阵，${\mathbf{H}}_k={\mathbf{h}}_k{\mathbf{h}}_k^H$，由矩阵论的基本结论可知
$${\Vert {\mathbf{w}}_k\Vert }^2=\text{trace}\left({\mathbf{B}}_k\right)$$
因此原问题可以等价地转化为：
$$\begin{gathered} \underset{{\mathbf{B}}_1,\dots ,{\mathbf{B}}_k}{\min }\sum _{k=1}^K\text{trace}\left({\mathbf{B}}_k\right) \\ \text{ s}\text{.t}\text{. trace}\left({\mathbf{H}}_k{\mathbf{B}}_k\right)-{\gamma }_k\sum _{j\ne k}\text{trace}\left({\mathbf{H}}_k{\mathbf{B}}_j\right)⩾{\gamma }_k{\sigma }^2, \\ {\mathbf{B}}_k⪰0,\forall k \end{gathered}$$
由于矩阵的迹为线性函数，因此可以使用CVX求解以上问题，但是需要注意的是当${\mathbf{B}}_k$秩为1时，进行简单的特征值分解即可，若秩大于1，需要利用高斯随机化过程等方法恢复。

KKT

和以上两类方法对比，KKT算法稍微复杂一些，但是可以获得解的结构，更为直观。首先，对于原问题，其拉格朗日函数为：
$$\mathcal{L}\left({\mathbf{w}}_1,\dots ,{\mathbf{w}}_K,{\lambda }_1,\dots ,{\lambda }_K\right)=\sum _{k=1}^K{\Vert {\mathbf{w}}_k\Vert }^2+\sum _{k=1}^K{\lambda }_k\left(\sum _{i\ne k}\frac{1}{{\sigma }^2}{\left|{\mathbf{h}}_k^H{\mathbf{w}}_i\right|}^2+1-\frac{1}{{\gamma }_k{\sigma }^2}{\left|{\mathbf{h}}_k^H{\mathbf{w}}_k\right|}^2\right)$$
其中${\lambda }_k\ge 0$是对应于第$k$个SINR约束的拉格朗日乘子。在此，在最优解处，满足KKT中的稳定性条件，$\frac{\partial \mathcal{L}}{\partial {\mathbf{w}}_k}=0$，即：（需要用到常用的矩阵求导公式：$\frac{\partial {\Vert {\mathbf{w}}_k\Vert }^2}{\partial {\mathbf{w}}_k}=2{\mathbf{w}}_k$，$\frac{\partial {\mathbf{w}}_k^H{\mathbf{h}}_k{\mathbf{h}}_k^H{\mathbf{w}}_k}{\partial {\mathbf{w}}_k}=2{\mathbf{h}}_k{\mathbf{h}}_k^H{\mathbf{w}}_k$）
w k + ∑ i ≠ k λ i σ 2 h i h i H w k − λ k γ k σ 2 h k h k H w k = 0 ⇔ ( I N + ∑ i = 1 K λ i σ 2 h i h i H ) w k = λ k σ 2 ( 1 + 1 γ k ) h k h k H w k ⇔ w k = ( I N + ∑ i = 1 K λ i σ 2 h i h i H ) − 1 h k λ k σ 2 ( 1 + 1 γ k ) h k H w k ⏟ 标量 

$$\begin{aligned} \mathbf{w}_k & +\sum_{i \neq k} \frac{\lambda_i}{\sigma^2} \mathbf{h}_i \mathbf{h}_i^H \mathbf{w}_k-\frac{\lambda_k}{\gamma_k \sigma^2} \mathbf{h}_k \mathbf{h}_k^H \mathbf{w}_k=\mathbf{0} \\ & \Leftrightarrow \quad\left(\mathbf{I}_N+\sum_{i=1}^K \frac{\lambda_i}{\sigma^2} \mathbf{h}_i \mathbf{h}_i^H\right) \mathbf{w}_k=\frac{\lambda_k}{\sigma^2}\left(1+\frac{1}{\gamma_k}\right) \mathbf{h}_k \mathbf{h}_k^H \mathbf{w}_k \\ & \Leftrightarrow \quad \mathbf{w}_k=\left(\mathbf{I}_N+\sum_{i=1}^K \frac{\lambda_i}{\sigma^2} \mathbf{h}_i \mathbf{h}_i^H\right)^{-1} \mathbf{h}_k \underbrace{\frac{\lambda_k}{\sigma^2}\left(1+\frac{1}{\gamma_k}\right) \mathbf{h}_k^H \mathbf{w}_k}_{\text {标量 }} \end{aligned}$$ wk​​+i=k∑​σ2λi​​hi​hiH​wk​−γk​σ2λk​​hk​hkH​wk​=0⇔(IN​+i=1∑K​σ2λi​​hi​hiH​)wk​=σ2λk​​(1+γk​1​)hk​hkH​wk​⇔wk​=(IN​+i=1∑K​σ2λi​​hi​hiH​)−1hk​标量  σ2λk​​(1+γk​1​)hkH​wk​​​​
由于$\frac{{\lambda }_k}{{\sigma }^2}\left(1+\frac{1}{{\gamma }_k}\right){\mathbf{h}}_k^H{\mathbf{w}}_k$是一个标量，因此最优${\mathbf{w}}_k$的方向一定是平行于向量${\left({\mathbf{I}}_N+{\sum }_{i=1}^K\frac{{\lambda }_i}{{\sigma }^2}{\mathbf{h}}_i{\mathbf{h}}_i^H\right)}^{-1}{\mathbf{h}}_k$，即足以有的波束向量${\mathbf{w}}_1^{\star },\dots ,{\mathbf{w}}_K^{\star }$具有以下的形式
$${\mathbf{w}}_k^{\star }=\underset{\text{=波束功率}}{\underbrace{\sqrt{p_k}}}\underset{={\tilde{\mathbf{w}}}_k^{\star }=\text{ 波束方向 }}{\underbrace{\frac{{\left({\mathbf{I}}_N+{\sum }_{i=1}^K\frac{{\lambda }_i}{{\sigma }^2}{\mathbf{h}}_i{\mathbf{h}}_i^H\right)}^{-1}{\mathbf{h}}_k}{\Vert {\left({\mathbf{I}}_N+{\sum }_{i=1}^K\frac{{\lambda }_i}{{\sigma }^2}{\mathbf{h}}_i{\mathbf{h}}_i^H\right)}^{-1}{\mathbf{h}}_k\Vert }}}\text{ for }k=1,\dots ,K$$
$p_k$表示分配给用户$k$的功率，${\tilde{\mathbf{w}}}_k^{\star }$表示归一化的波束方向。可以看出，各用户的波束方向可以直接根据计算得到，而需要计算每个用户所分配的功率。最优的功率分配在原问题SINR约束都取等号时取得，即$\frac{1}{{\gamma }_k}{p}_k{\left|{\mathbf{h}}_k^H{\tilde{\mathbf{w}}}_k^{\star }\right|}^2-{\sum }_{i\ne k}{p}_i{\left|{\mathbf{h}}_k^H{\tilde{\mathbf{w}}}_i^{\star }\right|}^2={\sigma }^2$（$\forall k=1,\dots ,K$），因此在这里只需要联合求解$K$个方程构成的方程组：
$$\left[\begin{array}{c}{p}_1\\ ⋮\\ p_K\end{array}\right]={\mathbf{M}}^{-1}\left[\begin{array}{c}{\sigma }^2\\ ⋮\\ {\sigma }^2\end{array}\right]\text{ 其中 }[\mathbf{M}{]}_{ij}=\{\begin{array}{ll}\frac{1}{{\gamma }_i}{\left|{\mathbf{h}}_i^H{\tilde{\mathbf{w}}}_i^{\star }\right|}^2,& i=j\\ -{\left|{\mathbf{h}}_i^H{\tilde{\mathbf{w}}}_j^{\star }\right|}^2,& i\ne j\end{array}$$
最后需要解决的是拉格朗日乘子${\lambda }_k$的值即可。根据稳态性条件处的推导，可知：
$${\lambda }_k=\frac{{\sigma }^2}{\left(1+\frac{1}{{\gamma }_k}\right){\mathbf{h}}_k^H{\left({\mathbf{I}}_N+{\sum }_{i=1}^K\frac{{\lambda }_i}{{\sigma }^2}{\mathbf{h}}_i{\mathbf{h}}_i^H\right)}^{-1}{\mathbf{h}}_k}$$
通过该方程，以及固定点算法不断迭代即可收敛得到最终的拉格朗日乘子，固定点算法是初始化一组${\lambda }_0$后去根据上式更新$\lambda$，不断迭代，直至收敛。

仿真结果：

以上方法和迫零并等功率分配算法进行对比，得到如下仿真结果:

可以看出几种方法得到的性能都近似，强于迫零算法。

参考文献

[1] Björnson E, Bengtsson M, Ottersten B. Optimal multiuser transmit beamforming: A difficult problem with a simple solution structure [lecture notes][J]. IEEE Signal Processing Magazine, 2014, 31(4): 142-148.

完整程序源码及图像下载地址

功率最小化问题Matlab代码.zip：下载地址