无线通信系统中的功率分配

作者：Robert Gowers，Roger Hill，Sami Al-Izzi，Timothy Pollington和Keith Briggs

来自Boyd和Vandenberghe的《凸优化》，第196页的练习4.20

凸优化可以用于最大化无线通信系统的最小信号干扰加噪声比（SINR）。考虑一个具有n个发射机的系统，每个发射机的功率为:p_jgeq 0。为了传输到n个接收机。让G_{ij}geq 0表示从发射机j到接收机i的路径增益。这些路径增益形成矩阵Gin mathbb{R}^{n times n}。

每个接收机被分配给一个发射机，使接收机i处的信号功率为S_i=G_{ii}p_i，干扰功率为I_i=sum_{kneq i} G_{ik}p_k。给定每个接收机的噪声功率sigma_i，接收机i处的SINR为gamma_i=frac{S_i}{I_i+sigma_i}。

目标是在一定的功率约束下最大化系统的最小SINR。这些约束条件如下：

i - 每个发射机的功率:p_j leq P_j^{text{max}}

ii - 如果发射机被划分为m个不重叠的组：K_1, …, K_m，它们与总功率P_l^{text{gp}}共享一个共同的电源，则sum_{kin K_l}p_k leq P_l^{text{gp}}。

iii - 每个接收机可以接收的最大功率为P_i^{text{rc}}，sum_{k=1}^{n}G_{ik}p_k leq P_i^{text{rc}}。

目标函数可以重写为：最小化：math:max_{i=1,…,n}frac{I_i + sigma_i}{S_i}

然而，由于这是一个拟凸的目标函数，我们不能直接使用CVXPY来解决它。相反，我们必须使用二分法。首先我们进行目标的重写，：math:alpha = gamma^{-1} geq 0，作为一个约束条件：

\(I_i+\sigma_i \leq S_i\alpha\)

然后我们选择初始下界和上界：math:L_0 和：math:U_0 为：math:alpha，应该选择这样的初始下界和上界使得：math:L < alpha^* < U，其中：math:alpha^* 是：math:alpha 的最优值。从一个初始值：math:alpha_0 = frac{1}{2}(L_0+U_0) 开始，使用任意目标函数来检查：math:`alpha_0`的可行性。从可行性中确定新的上界和下界：

如果：math:alpha_0 是可行的，那么：math:L_1 = L_0，：math:U_1 = alpha_0 和：math:alpha_1 = frac{1}{2}(L_1+U_1)。

如果：math:alpha_0 是不可行的，那么：math:L_1 = alpha_1，：math:U_1 = U_0 和：math:alpha_1 = frac{1}{2}(L_1+U_1)。

这个二分法过程一直重复，直到：math:U_N - L_N < epsilon，其中：math:epsilon 是所需的容忍度。

#!/usr/bin/env python3
# 作者：R. Gowers, S. Al-Izzi, T. Pollington, R. Hill & K. Briggs

import cvxpy as cp
import numpy as np

def maxmin_sinr(G, P_max, P_received, sigma, Group, Group_max, epsilon = 0.001):
    # 根据路径增益矩阵的大小找到 n 和 m
    n, m = np.shape(G)

    # 检查输入的大小
    if m != np.size(P_max):
        print('错误：P_max 的维度与增益矩阵的维度不匹配\n')
        return '错误：P_max 的维度与增益矩阵的维度不匹配\n', np.nan, np.nan, np.nan

    if n != np.size(P_received):
        print('错误：P_received 的维度与增益矩阵的维度不匹配\n')
        return '错误：P_received 的维度与增益矩阵的维度不匹配', np.nan, np.nan, np.nan

    if n != np.size(sigma):
        print('错误：σ 的维度与增益矩阵的维度不匹配\n')
        return '错误：σ 的维度与增益矩阵的维度不匹配', np.nan, np.nan, np.nan

    # I = np.zeros((n,m))
    # S = np.zeros((n,m))

    delta = np.identity(n)
    S = G*delta  # 信号功率矩阵
    I = G-S  # 干扰功率矩阵

    # 分组矩阵：按照发射机的数量分组
    num_groups = int(np.size(Group,0))

    if num_groups != np.size(Group_max):
        print('错误：Group 矩阵中的组数与 Group_max 的维度不匹配\n')
        return ('错误：Group 矩阵中的组数与 Group_max 的维度不匹配',
                np.nan, np.nan, np.nan, np.nan)

    # 将组的最大功率归一化为范围 [0,1]
    Group_norm = Group/np.sum(Group,axis=1).reshape((num_groups,1))

    # 创建标量优化变量 p：n 个发射机的功率
    p = cp.Variable(shape=n)
    best = np.zeros(n)

    # 设置子级集的上下界
    u = 1e4
    l = 0

    # alpha 定义广义线性二次问题的子级集，在这种情况下，α 是最小 SINR 的倒数
    alpha = cp.Parameter(shape=1)

    # 设置双分度可行性检验的约束条件
    constraints = [I*p + sigma <= alpha*S*p, p <= P_max, p >= 0, G*p <= P_received, Group_norm*p <= Group_max]

    # 定义目标函数，在这个案例中只希望测试解的可行性
    obj = cp.Minimize(alpha)

    # 现在检查解是否位于 u 和 l 之间
    alpha.value = [u]
    prob = cp.Problem(obj, constraints)
    prob.solve()

    if prob.status != 'optimal':
        # 在这种情况下，等级集 u 在解的下方
        print('上下界之间无最优解\n')
        return '错误：上下界之间无最优解', np.nan, np.nan, np.nan

    alpha.value = [l]
    prob = cp.Problem(obj, constraints)
    prob.solve()

    if prob.status == 'optimal':
        # 在这种情况下，等级集 l 在解的下方
        print('上下界之间无最优解\n')
        return '错误：上下界之间无最优解', np.nan, np.nan, np.nan

    # 双分度算法开始
    maxLoop = int(1e7)
    for i in range(1,maxLoop):
        # 首先检查 u 是否在可行域内，l 是否不在可行域内，如果不是这样，循环在此结束
        # 将 α 设置为区间的中点
        alpha.value = np.atleast_1d((u + l)/2.0)

        # 根据指定的容差测试区间的大小
        if u-l <= epsilon:
            break

        # 形成并求解问题
        prob = cp.Problem(obj, constraints)
        prob.solve()

        # 如果问题是可行的，则 u -> α，如果不是，则 l -> α，当达到容差时，新的 α 可能超出界限，最佳值取最后一个可行值作为最优值
        if prob.status == 'optimal':
            u = alpha.value
            best = p.value
        else:
            l = alpha.value

        # 最终条件检查区间是否收敛到顺序 ε，即最优子级集的范围 <=ε
        if u - l > epsilon and i == (maxLoop-1):
            print("解未收敛到顺序 epsilon")

    return l, u, float(alpha.value), best

示例

作为一个简单的示例，我们将考虑一种情况，其中 \(n=5\)，并且 \(G_{ij} = 0.6`（当 :math:`i=j\)）和 :math:`0.1`（其他情况）。

对于所有发射机 \(P_j^{\text{max}} = 1\)，并且发射机分为两组，每组 \(P_l^{\text{gp}} = 1.8\)。第一组包含发射机 1 和 2，而第二组包含发射机 3，4 和 5。

对于所有接收机 \(P_i^{\text{rc}} = 4\) 和 \(\sigma_i = 0.1\)。

np.set_printoptions(precision=3)

# 在这种情况下，我们将使用一个信号权重为 0.6 和干扰权重为 0.1 的增益矩阵
G = np.array([[0.6,0.1,0.1,0.1,0.1],
              [0.1,0.6,0.1,0.1,0.1],
              [0.1,0.1,0.6,0.1,0.1],
              [0.1,0.1,0.1,0.6,0.1],
              [0.1,0.1,0.1,0.1,0.6]])

# 在这种情况下，m=n，但是这个问题可以推广到我们想要 n 个接收机和 m 个发射机的情况
n, m = np.shape(G)

# 设置每个发射机和接收机的最大功率饱和级别
P_max = np.array([1.]*n)

# 归一化的接收功率，总可能性是所有发射机的所有功率，所以是 1/n
P_received = np.array([4.,4.,4.,4.,4.])/n

# 设置噪声水平
sigma = np.array([0.1,0.1,0.1,0.1,0.1])

# 分组矩阵：按照发射机的数量分组
Group = np.array([[1.,1.,0,0,0],[0,0,1.,1.,1.]])

# 最大归一化功率组，组数乘以 1
Group_max = np.array([1.8,1.8])

# 现在运行优化问题
l, u, alpha, best = maxmin_sinr(G, P_max, P_received, sigma, Group, Group_max)

print('最小 SINR={:.4g}'.format(1/alpha))
print('功率={}'.format(best))

Minimum SINR=1.148
Power=[0.8 0.8 0.8 0.8 0.8]