无线通信代码搬运/复现系列(4) ：MISO 广播信道的迫零脏纸编码波束形成器设计

Le-Nam Tran, Markku Juntti, Mats Bengtsson, and Björn Ottersten, “Beamformer Designs for MISO Broadcast Channels with Zero-forcing Dirty Paper Coding,” IEEE Trans. Wireless Commun., vol. 12, no. 3, pp. 1173-1185, Mar. 2013.

我们考虑使用零强迫脏纸编码（ZF-DPC）设计多输入多输出（MISO）广播信道（MISO BCs）的波束形成器。假设有一个总功率约束（SPC），大多数先前提出的波束形成器设计基于 QR 分解（QRD），这是满足 ZF 约束的自然选择。然而，基于 QRD 的 ZF-DPC 设计的最优性仍然未知。在本文中，我们首先从分析上确立 QRD 基于设计在任何性能度量下确实是 SPC 的最优。然后，我们提出了一种针对具有每天线功率约束（PAPCs）的 ZF-DPC 的最优波束形成器设计方法，使用凸优化框架。波束形成器设计首先被表述为一个秩-1 约束优化问题。利用 ZF-DPC 方案的特殊结构，我们证明秩约束可以放宽，并仍然提供相同的解决方案。此外，我们在广播信道（BCs）和多址信道（MACs）之间的对偶框架下，提出了一种快速收敛的波束形成器设计问题算法。 更具体地说，我们展示了具有 ZF-DPC 的 BC 具有基于 ZF 的连续干扰消除（ZF-SIC）的双重 MAC。通过这种方式，ZF-DPC 的波束形成器设计被转化为 ZF-SIC 的功率分配问题，这可以更高效地解决。

多输入多输出（MIMO）广播信道（BC）代表了一种点对多点传输，其中多天线发射器，例如蜂窝网络中的基站（BS），可以同时通过相同的系统资源向多个接收器发送数据。服务多个用户的能力得益于发射器和/或接收器的多天线所提供的空间复用增益。最近的信息理论研究证明，脏纸编码（DPC）实现了 MIMO BC 的容量区域 [3] 。通过 DPC，发射器处非因果已知的多用户干扰可以在不损失信息的情况下被完美消除。尽管这种多用户编码策略是最优的，但找到相应的最优发射协方差面临着高计算复杂性。因此，寻找 DPC 的简化替代方案一直是许多研究的主要焦点。
简化下行信道波束形成器设计的最常见方法是零强制（ZF）技术。在 [4] 中，提出了一种次优传输策略，现在称为 ZF-DPC，适用于多输入单输出（MISO）广播信道，将 ZF 技术与 DPC 相结合。具体来说，设 hk 和 wk 分别为 kth 用户的信道向量和波束形成器。那么，ZF 约束要求对所有 k<j 满足 hHkwj=0 。也就是说，由用户 k+1,k+2 …引起的第 k 个用户的干扰通过 ZF 方法被消除。这些 ZF 约束结合使用 DPC 来消除非因果已知干扰（由用户 1 到 k−1) 引起的干扰）将广播信道分解为一组并行的无干扰信道。除了 ZF 约束外，系统性能通常在一些功率约束下进行优化，例如总功率约束（SPC）或每天线功率约束（PAPC）

代码翻译及中文注释

function Dr = computeHessian(omega, psi, t)
global P power_pattern h_ul V K N

Dr11 = zeros(K, K);  % 初始化Hessian矩阵的第一区块
for k = 1:K
    Dr11(k, k) = -(((h_ul{k}') * inv(V{k}' * diag(psi) * V{k} + omega(k) * h_ul{k} * h_ul{k}') * h_ul{k})^2 + 1 / (t * omega(k)^2));
end
Dr11 = real(Dr11);  % 保留实部

Dr22 = zeros(N, N);  % 初始化Hessian矩阵的第二区块
for n = 1:N
    Dr22(n, n) = 1 / (t * psi(n)^2);
    for m = 1:N
        for k = 1:K
            vn = V{k}(n, :);  % 提取V矩阵的第n行
            vn = vn';
            vm = V{k}(m, :);  % 提取V矩阵的第m行
            vm = vm';
            
            Xi = inv(V{k}' * diag(psi) * V{k});  % 计算矩阵Xi
            Gamma = inv(V{k}' * diag(psi) * V{k} + omega(k) * h_ul{k} * h_ul{k}');  % 计算矩阵Gamma
            Dr22(n, m) = Dr22(n, m) - vn' * (Gamma * vm * vm' * Gamma - Xi * vm * vm' * Xi) * vn;
        end
    end
end
Dr22 = real(Dr22);  % 保留实部

Dr12 = zeros(K, N);  % 初始化Hessian矩阵的第三区块
for n = 1:K
    for m = 1:N
        Gamma = inv(V{n}' * diag(psi) * V{n} + omega(n) * h_ul{n} * h_ul{n}');  % 计算矩阵Gamma
        vm = V{n}(m, :);
        vm = vm';
        
        Dr12(n, m) = -h_ul{n}' * Gamma * vm * vm' * Gamma * h_ul{n};  % 计算Hessian矩阵的交叉项
    end
end
Dr12 = real(Dr12);  % 保留实部

% 组合最终的Hessian矩阵
Dr = [Dr11          Dr12             -ones(K, 1);
      Dr12'         Dr22              P * power_pattern;
      ones(1, K)    zeros(1, N)       0];

代码解释

这个函数 computeHessian 用于计算二阶导数矩阵（Hessian 矩阵），它在优化问题中用于牛顿法等二阶优化算法。具体来说，这个函数处理的是一个涉及多用户多输入多输出（MIMO）系统的优化问题，使用了拉格朗日乘子法和障碍方法。

1. 输入参数与全局变量：

• omega: 一个向量，代表每个用户的优化变量（如功率分配）。
• psi: 一个向量，代表与某些矩阵 V 相关的优化变量。
• t: 一个标量，通常与障碍方法中的阻尼系数相关。
• 全局变量 P、power_pattern、h_ul、V、K、N:
  • P: 与功率分配相关的一个标量。
  • power_pattern: 一个矩阵，涉及功率模式的配置。
  • h_ul: 每个用户的上行链路信道向量。
  • V: 每个用户的预编码矩阵。
  • K: 用户数量。
  • N: RIS或天线数量。

2. Hessian矩阵的计算：

• Dr11：K × K 的矩阵，对角元素与用户的 omega 相关。

  • 计算每个用户的二阶导数，对应 omega 变量。这个矩阵的元素是每个用户信道向量的二阶导数计算结果。

• Dr22：N × N 的矩阵，涉及与 psi 变量相关的二阶导数。

  • 对于每个天线或RIS元素，计算与 psi 相关的二阶导数。
  • 包含两个部分：一个与 psi 的对角项相关的简单项，另一个是通过矩阵 Gamma 和 Xi 的复杂计算。

• Dr12：K × N 的矩阵，涉及 omega 和 psi 之间的交叉项的二阶导数。

  • 计算 omega 和 psi 之间的交叉导数，涉及信道矩阵和预编码矩阵的乘积。

3. 最终Hessian矩阵的组合：

• 最后，函数将计算得到的 Dr11、Dr22 和 Dr12 组合成一个完整的Hessian矩阵 Dr，其中包含了优化问题中涉及的所有二阶导数信息。
• Dr 的结构是一个分块矩阵，由 Dr11、Dr22、Dr12 及其他一些固定项组合而成。

这个Hessian矩阵可以用于二阶优化算法，例如牛顿法，以便在每一步优化中更快地逼近问题的最优解。

代码翻译与中文注释

function r = computeresidualerror(omega, psi, mu, t)
global P power_pattern h_ul V K N
alpha = 2;  % 一个常数参数
r3 = sum(omega) - P;  % 计算第3个残差项
r1 = zeros(K, 1);  % 初始化r1的值
r2 = zeros(N, 1);  % 初始化r2的值

for k = 1:K
    % 计算r1的第k个元素
    r1(k) = (h_ul{k}') * inv(V{k}' * diag(psi) * V{k} + omega(k) * h_ul{k} * h_ul{k}') * h_ul{k} + 1 / (t * omega(k));
end
r1 = real(r1);  % 保留r1的实部

for n = 1:N
    r2(n) = -1 / (t * psi(n));  % 计算r2的第n个元素的初始值
    for k = 1:K
        v = V{k}(n, :);  % 提取V矩阵的第n行
        v = v';
        Xi = inv(V{k}' * diag(psi) * V{k});  % 计算矩阵Xi
        % 计算r2的第n个元素
        r2(n) = r2(n) - v' * omega(k) * Xi * h_ul{k} * h_ul{k}' * Xi * v / (1 + omega(k) * h_ul{k}' * Xi * h_ul{k});
    end
end
r2 = real(r2);  % 保留r2的实部

% 组合最终的残差向量r
r = [r1 - mu; r2 + P * power_pattern * mu; r3];

代码解释

这个函数 computeresidualerror 用于计算优化过程中不同约束条件下的残差向量。该残差向量在优化算法中通常用来判断是否满足某些约束或作为停止条件。

1. 输入参数与全局变量：

• omega: 每个用户的功率分配向量。
• psi: 用于预编码矩阵 V 的相关参数。
• mu: 对偶变量（通常用于拉格朗日乘子法）。
• t: 一个缩放参数，通常用于障碍函数法或对数障碍法中。
• 全局变量 P、power_pattern、h_ul、V、K、N:
  • P: 总功率约束。
  • power_pattern: 与功率模式相关的矩阵。
  • h_ul: 每个用户的上行链路信道向量。
  • V: 预编码矩阵。
  • K: 用户数量。
  • N: 天线或RIS元素数量。

2. 残差项的计算：

• r3:

  • r3 = sum(omega) - P; 表示功率约束的残差，即所有用户的功率总和与总功率约束 P 之间的差值。

• r1:

  • r1(k) 计算与 omega(k) 相关的残差。具体来说，它使用了矩阵 V{k} 和 h_ul{k} 计算得到信道增益的相关信息，并加入一个与 t 和 omega(k) 相关的项。

• r2:

  • r2(n) 计算与 psi(n) 相关的残差，涉及 V{k}、h_ul{k} 和 omega(k) 的复杂矩阵运算。这个部分计算了每个天线或RIS元素相关的残差。

3. 残差向量的组合：

• 最终的残差向量 r 由三个部分组合而成：
  • r1 - mu：与 omega 相关的残差与对偶变量 mu 的差值。
  • r2 + P * power_pattern * mu：与 psi 相关的残差和对偶变量 mu 的加权和。
  • r3：功率约束的残差。

总结

这个函数计算了一个优化问题中的残差向量，用于在每次迭代中判断优化的进展情况。残差向量结合了功率分配、预编码矩阵和对偶变量的信息，在优化算法中用于评估当前解的质量，并决定是否需要进一步迭代。

主程序

代码解释与中文注释

这个MATLAB代码实现了基于迭代法和CVX（Convex Optimization Toolbox）两种方法的波束成形设计，并比较了这两种方法的收敛性和最终结果。代码的目标是通过最大化和速率来设计MIMO系统中的波束成形器。

1. 基本设置与初始化

clear 
clc
global P power_pattern V  K N h_ul
rng(1)

N=4;  % 天线数量
K=4;  % 用户数量
P = 10;
P = 10.^(P./10);  % 发射功率（转换为线性刻度）
power_pattern=ones(N,1);  % 初始化功率模式
power_pattern=(power_pattern)/sum(power_pattern);  % 归一化功率模式

H = 1/2*(randn(K,N)+1i*randn(K,N));  % 随机生成信道矩阵

V=cell(K,1);  % 初始化V矩阵的cell数组
h_eff=cell(K,1);  % 初始化有效信道的cell数组
V{1}=eye(N);  % 第一个用户的V矩阵为单位矩阵
h_eff{1}=H(1,:);  % 第一个用户的有效信道为其原始信道

for k=2:K
    V{k}=null(H(1:k-1,:));  % V矩阵为前k-1个用户信道的正交补
    h_eff{k}=H(k,:)*V{k};  % 计算第k个用户的有效信道
end

h_ul=cell(K,1);
for k=1:K
    h_ul{k}=(H(k,:)*V{k})';  % 计算每个用户的上行链路信道向量
end

这部分代码初始化了系统的基本参数，包括天线数量、用户数量、发射功率等。然后随机生成了一个信道矩阵，并计算了每个用户的有效信道和上行链路信道向量。

2. 迭代法波束成形设计

BETA = .5;
ALPHA = .01;
MAXITERS = 100;
Sumrate = zeros(MAXITERS,1);
omega=ones(K,1);  % 初始化功率分配
psi=ones(N,1);    % 初始化功率模式
mu=0;             % 初始化对偶变量
t=1;
mygamma = 50;
n=0;

tic
while(t<1e20)
    for iters =1:MAXITERS
        r = computeresidualerror(omega,psi,mu,t);  % 计算残差
        if (norm(r) < 1e-8)
            break;
        end
        Dr=computeHessian(omega,psi,t);  % 计算Hessian矩阵
        step = -Dr\r;  % 计算更新步长
        domega = step(1:K);  % 提取功率分配更新步长
        dpsi = step(K+(1:N));  % 提取功率模式更新步长
        dmu=step(K+N+1:end);  % 提取对偶变量更新步长
        s = 1;
        newomega = omega+s*domega; newpsi = psi+s*dpsi;
        while ((min(newomega) < 0) || (min(newpsi) < 0)),
            s = BETA*s;
            newomega = omega+s*domega;
            newpsi = psi+s*dpsi;
        end;
        
        newr = computeresidualerror(newomega,newpsi,mu+s*dmu,t);  % 计算新残差
        while (norm(newr) > (1-ALPHA*s)*norm(r)) % 回溯线搜索
            s = BETA*s;
            newomega = omega+s*domega;
            newpsi = psi+s*dpsi;
            newr = computerresidualerror(newomega,newpsi,mu+s*dmu,t);
        end
        omega = omega+s*domega; psi = psi+s*dpsi;
        mu=mu+dmu*s;
        n=n+1;
        Sumrate(n)=0;
        for k=1:K
            Sumrate(n)=Sumrate(n)+log(real(det(V{k}' * diag(psi)*V{k}+h_ul{k}*omega(k)*h_ul{k}')))-log(real(det(V{k}' * diag(psi)*V{k})));
        end
    end
    t=t*mygamma;
end
toc
Sumrate(n:end)=[];
gap=abs(Sumrate-Sumrate(end));
semilogy(1:n-1,gap)
ylabel('Residual error')
xlabel('Iteration')
saveas(gcf,'../results/convergence.png')

这部分代码使用迭代方法来设计波束成形器。通过计算残差、Hessian矩阵、以及更新步长，逐步优化功率分配 omega 和功率模式 psi，以最大化和速率。迭代的过程中会进行回溯线搜索以确保每一步更新都是有效的。

3. CVX求解比较

cvx_expert true
tic
cvx_begin  sdp quiet
    variable covarmat(sum([N:-1:N-K+1]),sum([N:-1:N-K+1])) complex hermitian
    obj=0;
    power=0;
    idx=1;
    for k=1:K
        covarmat(idx:idx+N-k,idx:idx+N-k) >= 0;
        obj=obj+log(1+real(h_eff{k} * covarmat(idx:idx+N-k,idx:idx+N-k) * h_eff{k}'));
        power=power+diag(V{k} * covarmat(idx:idx+N-k,idx:idx+N-k) * V{k}');
        idx=idx+N-k+1;
    end
    maximize( (obj));
    subject to
    real(power) <= P*power_pattern;
cvx_end
toc
Sumrate_CVX=real(obj);
abs(Sumrate(end)-Sumrate_CVX)

这部分代码使用CVX（凸优化工具箱）求解相同的波束成形优化问题，以作为对比。它直接求解优化问题的目标函数和约束条件，计算得到的和速率 Sumrate_CVX 被用来与迭代方法的结果 Sumrate(end) 进行比较，以评估迭代方法的性能和收敛性。

总结

该代码通过两种方法设计MIMO系统中的波束成形器，并比较了两种方法的性能。迭代法通过逐步优化实现波束成形设计，而CVX方法则通过求解凸优化问题直接得到结果。最终的结果比较表明了迭代方法的有效性和收敛性。