无线通信代码搬运/复现系列(5) ：多用户 MISO 下行的最佳能效传输波束成形

Oskari Tervo, Le-Nam Tran, and Markku Juntti, “Optimal Energy-Efficient Transmit Beamforming for Multi-User MISO Downlink,” IEEE Trans. Signal Process., vol. 63, no. 20, pp. 5574-5587, Oct. 2015.

本文研究了多用户多输入单输出（MISO）下行系统中用于能效最大化（EEmax）的波束成形技术。针对这一具有挑战性的非凸问题，我们首先采用分支-减少-界限（BRB）方法推导出最优解。我们还提出了两种低复杂度的近似设计。第一种设计使用著名的零强迫波束成形（ZFBF）来消除用户间干扰，从而将 EEmax 问题简化为一个凹凸分数规划。特别地，该问题通过结合丁克尔巴赫方法的闭式表达式高效求解。在第二种设计中，我们旨在使用序列凸近似（SCA）方法找到一个驻点。通过适当的变换，我们得到了一个快速收敛的迭代算法，其中在每次迭代中解决一个凸规划。我们进一步表明，每次迭代中的问题也可以近似为二阶锥规划（SOCP），从而利用计算效率高的最先进 SOCP 求解器。数值实验表明，第二种设计收敛迅速，并实现了接近最优的性能。 为了进一步提高能效，我们还考虑了联合波束形成和天线选择（JBAS）问题，并提出了两种设计。在第一种方法中，我们利用视角重构结合连续松弛来解决 JBAS 问题。在第二种方法中，引入稀疏诱导正则化来近似 JBAS 问题，然后通过 SCA 方法进行求解。数值结果表明，联合波束形成和天线选择在大量发射天线的情况下显著提高了能效。

SPmin.m
这段代码实现了一个优化算法，用于设计多用户MIMO系统的波束成形器，使得在满足每个用户信噪比（SINR）要求的前提下，最小化发射功率。该优化问题通过使用YALMIP工具箱和MOSEK求解器解决。

function [minpower, beamformer] = SPmin(channel, sinrthreshold)

[nUsers, nTx] = size(channel); % 获取用户数量和发射天线数量

% 设置优化器，选择MOSEK求解器，并关闭详细输出
ops = sdpsettings('solver', 'mosek', 'verbose', 0); 

% 定义优化变量：波束成形矩阵，大小为 nTx x nUsers，复数类型
beamformer = sdpvar(nTx, nUsers, 'full', 'complex'); 

F = []; % 初始化约束条件集合
obj = norm(vec(beamformer)); % 定义目标函数为波束成形矩阵的范数（即功率的平方根）

for iUser = 1:nUsers
    % 添加约束条件1：确保波束成形器设计的信号实部有效
    F = [F, imag(channel(iUser, :) * beamformer(:, iUser)) == 0];
    
    % 添加约束条件2：确保每个用户的SINR门限得以满足
    F = [F, real(channel(iUser, :) * beamformer(:, iUser)) >= ...
        sinrthreshold(iUser) * norm([1, channel(iUser, :) * beamformer(:, 1:nUsers ~= iUser)])];
end

% 优化问题求解
diagnotics = optimize(F, obj, ops);

% 检查求解结果是否成功
if (diagnotics.problem == 0)
    % 若成功，计算最小发射功率，并返回波束成形矩阵
    minpower = double(obj)^2;
    beamformer = value(beamformer);
else
    % 若失败，返回NaN表示求解失败
    minpower = NaN;
    beamformer = NaN;
end

代码解释

1. 输入参数:

   • channel: 信道矩阵，大小为 nUsers x nTx，表示每个用户与每个发射天线之间的信道增益。
   • sinrthreshold: 每个用户的信噪比（SINR）门限，大小为 nUsers x 1。

2. 输出参数:

   • minpower: 最小化后的发射功率。
   • beamformer: 设计得到的波束成形矩阵，大小为 nTx x nUsers。

3. 算法流程:

   • 步骤1: 通过 size 函数获取信道矩阵的尺寸，即用户数量和发射天线数量。
   • 步骤2: 使用 sdpsettings 函数设置优化求解器为 MOSEK，并关闭详细输出。
   • 步骤3: 定义优化变量 beamformer，表示波束成形矩阵，大小为 nTx x nUsers，为复数矩阵。
   • 步骤4: 初始化约束条件集合 F，并定义目标函数 obj，即波束成形矩阵的范数。
   • 步骤5: 对于每个用户，添加两个约束条件：
     • 约束1: 确保波束成形矩阵在该用户的信道上的虚部为零，确保传输信号的实部有效。
     • 约束2: 确保每个用户的SINR大于或等于给定的门限。
   • 步骤6: 使用 optimize 函数求解优化问题，若求解成功，则返回最小化后的发射功率和波束成形矩阵；若求解失败，则返回 NaN 表示失败。

总结

该代码利用凸优化工具箱YALMIP和MOSEK求解器，设计了一个优化算法来最小化MIMO系统中的发射功率，并确保每个用户的信噪比要求得以满足。通过定义合适的约束条件和目标函数，该算法能够有效地设计波束成形器，从而在满足SINR要求的前提下最小化功率消耗。

一阶近似

function output = firstorderapp(x,y,x0,y0)
% approximate \sqrt(xy) around (x0,y0)
%   Detailed explanation goes here
output = sqrt(x0*y0)+1/2*sqrt(x0/y0)*(y-y0)+1/2*sqrt(y0/x0)*(x-x0);
end

EEmaxMISOdownlinkoptimalBB
能效最大化 分支定界

function [beamformers, currentbestUBstored, currentbestobjectivestored] = EEmaxMISOdownlinkoptimalBB(channel, sinrthreshold, power, Po, PAeff, tol)

[nUsers, nTx] = size(channel);
beamformers = zeros(nTx, nUsers);

% 设置 t 的上下限
tmin = [1/(Po + power/PAeff); log(1 + sinrthreshold)];
tmax = [1/Po; log(1 + power * diag(channel * channel'))];

% 计算初始盒子
newinitalbox = computeimprovedbox(channel, power, Po, PAeff, [tmin tmax], 1e-3);

gap = 0;
ms = 3;

% 初始最优目标值和上界
currentbestobj = newinitalbox(1,1) * sum(newinitalbox(2:end, 1));
currentbestUB = newinitalbox(1,2) * sum(newinitalbox(2:end, 2));

sequenceoflowerbounds = currentbestobj;
sequenceofupperbounds = currentbestUB;

currentbestobjectivestored = currentbestobj;
currentbestUBstored = currentbestUB;

nBoxes = 1;
boxstored(:,:,1) = newinitalbox;

nColors = 20;
cc = hsv(nColors);
iColor = 1;
nIterations = 1;

% 主循环：当上界与当前最优目标值的差距大于容忍度时
while ((currentbestUB - currentbestobj) > tol)
    currentbestUB - currentbestobj
    for iBox = 1:nBoxes
        if (sequenceofupperbounds(iBox) == currentbestUB)
            % 分支操作
            lowerlimit = boxstored(:,1,iBox);
            upperlimit = boxstored(:,2,iBox);

            [~, dim] = max(upperlimit - lowerlimit);  % 找到最大维度用于分支

            % 计算新的上下限
            newupperlimit = upperlimit;
            newupperlimit(dim) = (upperlimit(dim) + lowerlimit(dim)) / 2;

            newlowerlimit = lowerlimit;
            newlowerlimit(dim) = (upperlimit(dim) + lowerlimit(dim)) / 2;

            newbox1 = [lowerlimit newupperlimit];
            newbox2 = [newlowerlimit upperlimit];

            % 改进新盒子
            improvednewbox1 = computeimprovedbox(channel, power, Po, PAeff, newbox1, 1e-3);
            improvednewbox2 = computeimprovedbox(channel, power, Po, PAeff, newbox2, 1e-3);

            % 存储盒子
            boxstored(:,:,iBox) = improvednewbox1;
            boxstored(:,:,1+nBoxes) = improvednewbox2;
            iColor = iColor + 1;

            % 计算新的上下界
            sequenceofupperbounds(iBox) = improvednewbox1(1,2) * sum(improvednewbox1(2:end, 2));
            sequenceofupperbounds(1+nBoxes) = improvednewbox2(1,2) * sum(improvednewbox2(2:end, 2));

            sequenceoflowerbounds(iBox) = improvednewbox1(1,1) * sum(improvednewbox1(2:end, 1));
            sequenceoflowerbounds(1+nBoxes) = improvednewbox2(1,1) * sum(improvednewbox2(2:end, 1));

            % 更新全局上下界
            currentbestobj = max(sequenceoflowerbounds);
            currentbestUB = max(sequenceofupperbounds);

            nBoxes = nBoxes + 1;

            % 剪枝操作
            prunedboxes = find(currentbestobj > sequenceofupperbounds);
            if (~isempty(prunedboxes))
                boxstored(:,:,prunedboxes) = [];
                sequenceoflowerbounds(prunedboxes) = [];
                sequenceofupperbounds(prunedboxes) = [];
                nBoxes = nBoxes - length(prunedboxes);
            end

            nIterations = nIterations + 1;
            currentbestobjectivestored(nIterations) = currentbestobj;
            currentbestUBstored(nIterations) = currentbestUB;
            break;
        end
    end
end

% 在最终的盒子中找到最优波束成形器
iBox = find(sequenceoflowerbounds == currentbestobj);
t = boxstored(:,1,iBox);

cvx_begin quiet
    variable beamformers(nTx, nUsers) complex
    subject to
    for iUser = 1:nUsers
        imag(channel(iUser,:) * beamformers(:, iUser)) == 0;
        real(channel(iUser,:) * beamformers(:, iUser)) >= ...
            sqrt(exp(t(iUser+1)) - 1) * norm([1 channel(iUser,:) * beamformers(:, 1:nUsers ~= iUser)]);
    end
    norm(vec(beamformers)) <= sqrt(min(power, PAeff * (1/t(1) - Po)));
cvx_end
cvx_status

end

代码解析

1. 输入参数:

   • channel: 信道矩阵，表示每个用户的信道增益。
   • sinrthreshold: 每个用户的信噪比门限。
   • power: 发射功率上限。
   • Po: 恒定功率消耗。
   • PAeff: 功率放大器效率。
   • tol: 容忍度，用于判断算法的收敛性。

2. 初始盒子的设置:

   • tmin 和 tmax 定义了 t 的初始上下界。
   • newinitalbox 通过调用 computeimprovedbox 函数计算改进后的初始盒子。

3. 主循环:

   • 循环的主要目的是在最优值和上界之间找到一个符合容忍度的解。
   • 在每次循环中，算法会分支（branching）当前的盒子，生成新的盒子。
   • 每个新盒子通过调用 computeimprovedbox 进一步改进。
   • 通过计算新的上下界，更新全局的最优目标值和上界。
   • 使用剪枝操作（pruning）删除不必要的盒子。

4. 最终求解:

   • 最后，使用CVX求解器在最终的盒子中找到最优的波束成形器。

总结

这个代码使用盒子方法和分支界定技术来求解MISO系统的能效最大化问题。通过不断分割搜索空间，并使用凸优化工具箱（CVX）求解最终的优化问题，算法能够找到一个符合指定容忍度的最优解。

该函数 computeimprovedbox 的作用是通过调整给定的盒子上下限，改进一个边界框，使其尽可能紧凑，并且满足一定的可行性条件。这个函数主要通过二分搜索的方法来调整盒子的上下限。

1. 函数头部与输入参数

function [improvedbox] = computeimprovedbox(controllerisfeasibility, channel, box, tol)

• 输入参数:

  • controllerisfeasibility: 用于检测给定参数是否可行的控制器函数。
  • channel: 信道矩阵，用于描述用户与天线之间的连接关系。
  • box: 初始盒子的上下限，第一列是下限，第二列是上限。
  • tol: 容忍度，用于二分搜索时的停止条件。

• 输出参数:

  • improvedbox: 改进后的盒子，其上下限更为紧凑。

2. 改进上限

improvedlowerlimit = box(:,1);
improvedupperlimit = box(:,2);

for dim = 1:length(improvedlowerlimit)
    t = improvedlowerlimit;
    t(dim) = improvedupperlimit(dim);
    
    [~, problem] = controllerisfeasibility{{channel, [t(1); sqrt(exp(t(2:end)) - 1)]}};
    if(problem ~= 0)  % 如果上限不可行，使用二分搜索
        tmax = improvedupperlimit(dim);
        tmin = improvedlowerlimit(dim);
        while (1)
            t(dim) = (tmax + tmin) / 2;
            [~, problem] = controllerisfeasibility{{channel, [t(1); sqrt(exp(t(2:end)) - 1)]}};
            if(problem == 0)  % 如果可行
                tmin = t(dim);
                if((tmax - tmin) < tol)
                    break;
                end
            else
                tmax = t(dim);
            end
        end
    end
    improvedupperlimit(dim) = t(dim);
end

• 首先初始化盒子的上下限。
• 对每一个维度，尝试将该维度的上限设为当前值，并检查是否可行。
• 如果不可行，通过二分搜索不断缩小上限，直到找到一个尽可能大的可行值。

3. 改进下限

for dim = 1:length(improvedlowerlimit)
    t = improvedupperlimit;
    t(dim) = improvedlowerlimit(dim);
    
    [~, problem] = controllerisfeasibility{{channel, [t(1); sqrt(exp(t(2:end)) - 1)]}};
    if(problem == 0)  % 如果下限可行，使用二分搜索
        tmax = improvedupperlimit(dim);
        tmin = improvedlowerlimit(dim);
        while (1)
            t(dim) = (tmax + tmin) / 2;
            [~, problem] = controllerisfeasibility{{channel, [t(1); sqrt(exp(t(2:end)) - 1)]}};
            if(problem == 0)  % 如果可行
                tmin = t(dim);
                if((tmax - tmin) < tol)
                    break;
                end
            else
                tmax = t(dim);
            end
        end
    end
    improvedlowerlimit(dim) = t(dim);
end

• 同样的过程，用于改进下限。
• 对每一个维度，尝试将该维度的下限设为当前值，并检查是否可行。
• 如果可行，通过二分搜索不断增大下限，直到找到一个尽可能小的可行值。

4. 返回改进后的盒子

improvedbox = [improvedlowerlimit improvedupperlimit];

• 返回改进后的盒子，该盒子比原来的更为紧凑，并且满足给定的可行性条件。

代码总结

这个函数 computeimprovedbox 的目的是通过在每个维度上调整盒子的上下限，使得盒子变得尽可能小，同时仍然满足给定的可行性条件。它使用了二分搜索的策略，在保证解的精度的同时，逐步缩小搜索空间，以找到一个更为优化的边界框。这种方法在优化问题中特别有用，能够有效地减少搜索空间，提高算法的效率。

该函数 computeEE 用于计算能效（Energy Efficiency, EE），即在多用户 MISO 系统中，信道容量与功率消耗的比值。函数接受信道矩阵、波束成形矩阵以及功率消耗相关参数作为输入，返回系统的能效。（目标函数）

function energyefficiency = computeEE(channel,beamformer,Po,PAeff)
T = channel*beamformer;
signal = diag(T); % diagonal entries are h_k*w_k
signal = abs(signal).^2; % |h_k*w_k|^2 for each k
T(logical(eye(size(T)))) = 1; % replace the diagonals by 1
interference = sum(T.*conj(T),2); % norm of each row is the interference for each user
sumrate = sum(log(1 + signal./interference)); % rats in nat/s/Hz
energyefficiency = real(sumrate/(Po + (1/PAeff)*norm(beamformer(:))^2));

这个函数 Algorithm3_SCA 实现了一个基于顺序凸近似（SCA）的迭代算法，用于在多用户MISO系统中设计能效最大化的波束成形器。目标是最大化系统的能效（Energy Efficiency, EE），同时满足每个用户的信噪比（SINR）要求。

function [EEbeamforming, beamformeropt, achievedpower, userrate] = Algorithm3_SCA(channel, sinrthreshold, power, Po, PAeff, beamformerinit)
    [nUsers, nTx] = size(channel);
    mybeta0 = zeros(nUsers,1);
    myvartheta0 = zeros(nUsers,1);
    myalphanext = zeros(nUsers,1);

    % 初始化参数
    for iUser = 1:nUsers
        mybeta0(iUser) = norm([1, channel(iUser,:)*beamformerinit(:, 1:nUsers~=iUser)])^2;
        myvartheta0(iUser) = (real(channel(iUser,:)*beamformerinit(:,iUser)) / sqrt(mybeta0(iUser)))^2 + 1;
        myalphanext(iUser) = log(myvartheta0(iUser));
    end

    z0 = (Po + (1/PAeff) * norm(beamformerinit(:))^2)^2;
    t0 = sum(myalphanext)^2 / z0;
    nIterations = 20;

    seqobj = zeros(nIterations+1, 1);
    EEbeamforming = zeros(nIterations+1, 1);
    EEbeamforming(1) = computeEE(channel, beamformerinit, Po, PAeff);
    seqobj(1) = sqrt(t0);

    % 定义优化变量
    beamformer = sdpvar(nTx, nUsers, 'full', 'complex');
    myvartheta = sdpvar(nUsers, 1);
    mybeta = sdpvar(nUsers, 1); 
    myalpha = sdpvar(nUsers, 1);
    sdpvar t z1 z
    obj = t;

    % 设置优化选项
    opts = sdpsettings('solver', 'mosek', 'verbose', 0);

    % 迭代优化过程
    for iIteration = 1:nIterations
        F = [];
        for iUser = 1:nUsers
            % 添加约束条件 (8)
            F = [F, imag(channel(iUser,:) * beamformer(:,iUser)) == 0];

            % 添加第二个约束条件，确保信噪比满足要求
            F = [F, cone([1, channel(iUser,:)*beamformer(:,1:nUsers~=iUser)], ...
                real(channel(iUser,:)*beamformer(:,iUser)) / sqrt(sinrthreshold(iUser)))];

            % 添加约束条件 (26b)
            F = [F, myvartheta(iUser) >= exp(myalpha(iUser))];

            % 添加约束条件 (27b)
            F = [F, cone([1, channel(iUser,:)*beamformer(:,1:nUsers~=iUser), (mybeta(iUser)-1)/2], (mybeta(iUser)+1)/2)];

            % 添加约束条件 (30b)，基于一阶近似
            F = [F, real(channel(iUser,:)*beamformer(:,iUser)) >= firstorderapp(myvartheta(iUser)-1, mybeta(iUser), myvartheta0(iUser)-1, mybeta0(iUser))];
        end

        % 添加约束条件 (22c)
        F = [F, cone([z1; (z-1)/2], (z+1)/2)];
        F = [F, cone([sqrt(1/PAeff)*beamformer(:); (z1-Po-1)/2], (z1-Po+1)/2)];

        % 添加基于一阶近似的目标函数约束
        F = [F, sum(myalpha) >= firstorderapp(t, z, t0, z0)];

        % 添加功率约束 (5b)
        F = [F, cone(beamformer(:), sqrt(power))];

        % 求解优化问题
        diagontics = optimize(F, -obj, opts);

        if(diagontics.problem == 0 || diagontics.problem == 4)
            % 如果求解成功，更新参数
            beamformer0 = value(beamformer);
            t0 = value(t);
            z0 = value(z);
            myvartheta0 = value(myvartheta);
            mybeta0 = value(mybeta);
            seqobj(iIteration+1) = sqrt(value(obj));
            EEbeamforming(iIteration+1) = computeEE(channel, beamformer0, Po, PAeff);
            achievedpower = norm(beamformer0(:))^2;
            userrate = zeros(nUsers, 1);
            for iUser = 1:nUsers
                userrate(iUser) = log(1 + abs(channel(iUser,:) * beamformer0(:,iUser))^2 ...
                    / (norm([1, channel(iUser,:)*beamformer0(:,1:nUsers~=iUser)])^2));
            end
        else
            % 如果求解失败，返回NaN
            EEbeamforming = NaN; 
            achievedpower = NaN; 
            seqobj = NaN; 
            beamformer0 = NaN;
        end
    end

    beamformeropt = beamformer0;
end

代码总结

1. 初始化: 函数首先根据初始波束成形矩阵计算相关的参数，如 mybeta0 和 myvartheta0，并初始化目标函数的初值。

2. 定义优化问题: 在每次迭代中，使用YALMIP和MOSEK设置优化问题的约束和目标函数，并根据SCA方法逐步改进波束成形矩阵。

3. 迭代过程: 每次迭代更新波束成形矩阵、SINR相关参数以及目标函数，直到达到设定的最大迭代次数。

4. 结果输出: 返回优化后的能效值、最终的波束成形矩阵、达成的总功率以及用户速率。

该函数能够有效地优化多用户MISO系统的波束成形器，使得在满足SINR约束的同时，能效达到最大化。

这个函数 Algorithm1_BRB 实现了基于分支定界（Branch and Bound, BRB）的算法，用于在多用户MISO系统中寻找最优的波束成形器设计，以最大化能效。该算法通过逐步分割（Branching）和约束（Bounding）搜索空间，并通过剪枝（Pruning）来减少搜索空间的大小。

1. 函数头部与输入参数

function [optbeamformers, currentbestUBstored, currentbestobjectivestored] = Algorithm1_BRB(controllerisfeasibility, channel, sinrthreshold, power, Po, PAeff, tol)

• 输入参数:

  • controllerisfeasibility: 用于检测给定参数是否可行的控制器函数。
  • channel: 信道矩阵，描述每个用户和每个发射天线之间的信道增益。
  • sinrthreshold: 每个用户的信噪比门限。
  • power: 发射功率上限。
  • Po: 恒定功率消耗。
  • PAeff: 功率放大器效率。
  • tol: 容忍度，用于判断算法的收敛性。

• 输出参数:

  • optbeamformers: 计算得到的最优波束成形矩阵。
  • currentbestUBstored: 迭代过程中存储的上界序列。
  • currentbestobjectivestored: 迭代过程中存储的目标函数最优值序列。

2. 初始盒子设置

[nUsers, nTx] = size(channel);

tmin = [1/(Po + power/PAeff); log(1 + sinrthreshold)];
tmax = [1/Po; log(1 + power * diag(channel * channel'))];

newinitalbox = computeimprovedbox(controllerisfeasibility, channel, [tmin tmax], 1e-3);

currentbestobj = newinitalbox(1,1) * sum(newinitalbox(2:end,1));
currentbestUB = newinitalbox(1,2) * sum(newinitalbox(2:end,2));

• 这里初始化了上下界 tmin 和 tmax，并调用 computeimprovedbox 函数生成初始的搜索盒子。
• currentbestobj 和 currentbestUB 分别初始化为目标函数的当前最优值和上界。

3. 主循环：分支和界定

for iIter = 1:maxIteration
    [~,iBoxes] = max(sequenceofupperbounds);
    iBox = iBoxes(randi(length(iBoxes)));  % 随机选择一个最大上界的盒子

    % 分支操作
    lowerlimit = boxstored(:,1,iBox);
    upperlimit = boxstored(:,2,iBox);
    [~,dim] = max(upperlimit - lowerlimit);  % 选择最大维度进行分支

    newupperlimit = upperlimit;
    newupperlimit(dim) = (upperlimit(dim) + lowerlimit(dim)) / 2;

    newlowerlimit = lowerlimit;
    newlowerlimit(dim) = (upperlimit(dim) + lowerlimit(dim)) / 2;

    newbox1 = [lowerlimit newupperlimit];
    newbox2 = [newlowerlimit upperlimit];

    % 改进新盒子
    t = newbox1(:,1);
    [~, problem] = controllerisfeasibility{{channel, [t(1); sqrt(exp(t(2:end)) - 1)]}};
    if(problem == 0)
        improvednewbox1 = computeimprovedbox(controllerisfeasibility, channel, newbox1, 1e-3);
    end

    t = newbox2(:,1);
    [~, problem] = controllerisfeasibility{{channel, [t(1); sqrt(exp(t(2:end)) - 1)]}};
    if(problem == 0)
        improvednewbox2 = computeimprovedbox(controllerisfeasibility, channel, newbox2, 1e-3);
    end

    % 存储盒子
    boxstored(:,:,iBox) = improvednewbox1;
    boxstored(:,:,1+nBoxes) = improvednewbox2;

    % 更新上下界
    sequenceofupperbounds(iBox) = improvednewbox1(1,2) * sum(improvednewbox1(2:end,2));
    sequenceofupperbounds(1+nBoxes) = improvednewbox2(1,2) * sum(improvednewbox2(2:end,2));

    sequenceoflowerbounds(iBox) = improvednewbox1(1,1) * sum(improvednewbox1(2:end,1));
    sequenceoflowerbounds(1+nBoxes) = improvednewbox2(1,1) * sum(improvednewbox2(2:end,1));

    currentbestobj = max(sequenceoflowerbounds);
    currentbestUB = max(sequenceofupperbounds);

    nBoxes = nBoxes + 1;

    % 剪枝操作
    prunedboxes = find(currentbestobj > sequenceofupperbounds);
    boxstored(:,:,prunedboxes) = [];
    sequenceoflowerbounds(prunedboxes) = [];
    sequenceofupperbounds(prunedboxes) = [];
    nBoxes = nBoxes - length(prunedboxes);

    currentbestobjectivestored(iIter+1) = currentbestobj;
    currentbestUBstored(iIter+1) = currentbestUB;

    % 检查收敛条件
    if(abs(currentbestUB - currentbestobj) < tol)
        break
    end
end

• 这里的主要工作是逐步分支和更新搜索盒子，然后通过剪枝来减少不必要的计算。
• 在每次迭代中，算法会选择一个具有最大上界的盒子进行分支，并进一步更新其上界和下界。

4. 求解最优波束成形器

iBox = find(sequenceoflowerbounds == currentbestobj);
t = boxstored(:,1,iBox);
opts = sdpsettings('solver', 'mosek', 'verbose', 0);
beamformers = sdpvar(nTx, nUsers, 'full', 'complex');
F = [];

for iUser = 1:nUsers
    F = [F, imag(channel(iUser,:) * beamformers(:,iUser)) == 0];
    F = [F, cone([1, channel(iUser,:)*beamformers(:,1:nUsers~=iUser)], ...
        real(channel(iUser,:)*beamformers(:,iUser)) / sqrt(exp(t(iUser+1))-1))];
end

F = [F, cone(beamformers(:), sqrt(min(power, PAeff * (1/t(1) - Po))))];
diagnotics = optimize(F, [], opts);

if (diagnotics.problem == 0) || (diagnotics.problem == 4)
    optbeamformers = value(beamformers);
else
    optbeamformers = NaN;
end

• 这里通过CVX求解最终的波束成形优化问题，以获得最优的波束成形矩阵。
• 约束条件确保每个用户的信号实部有效且满足SINR要求，同时确保波束成形矩阵的功率在给定的上限内。

总结

这个 Algorithm1_BRB 函数通过分支定界算法逐步搜索最优波束成形器设计，最终在满足功率和SINR约束的前提下，最大化系统能效。它结合了迭代优化和剪枝策略，有效地减少了搜索空间，提高了计算效率。

主程序

clear 
clc
rng('default')
nUsers = 3; % number of users
nTx = 4; % number of tx antennas

%% power consumption
power = 46; % dB
power = 10.^((power-30)/10); % to linear scale
PAeff = 0.38; %  power amplifier efficiency
Pdyn = 42; % dynamic power consumption (dBm)
Pdyn = 10^((Pdyn-30)/10); % to linear scale
Psta = 33; % dBm
Psta = 10^((Psta-30)/10); % to linear scale
Po = nTx*Pdyn+Psta; % total circuit power


sinrthreshold = ones(nUsers,1); % SINR requirements (i.e., \bar{\gamma}_k in (4b)

%% path loss
d = 0.1+0.7*rand(nUsers,1) ; % users's distance: .1km to .8 km
PL = 128.1+37.6*log10(d); % pass loss model (in dB)
PL = 10.^(-PL/10);
B = 10e6; % bandwidth
No = -203.975; % dBW/Hz
noise_p = No+10*log10(B)+30; % (~104 dBm)

noise_p = 10^((noise_p-30)/10); % to linear scale

%% main loop here

% generate small scale fading channel
channel = sqrt(PL).*(sqrt(1/2)*(randn(nUsers,nTx)+1i*randn(nUsers,nTx)));
channel = channel/sqrt(noise_p); % normalize the channel with noise power

% the noise power in SINR expressions now becomes 1

%% SPmin power
[minpowerSPmin,beamformerSPmin]= SPmin(channel,sinrthreshold);
EEofSPmin = computeEE(channel,beamformerSPmin,Po,PAeff);

%% Algorithm 3 (SCA)

[EE_SCA_seq,beamformerEEmax_SCA]= Algorithm3_SCA(channel,sinrthreshold,power,Po,PAeff,beamformerSPmin);

%% Algorithm 1 (BRB)
% In Algorithm 1, we need to REPEATEDLY check if (10) is feasible or not.
% For this purpose, we first create an optimizer for (10)

channelpara = sdpvar(nUsers, nTx,'full','complex'); % treat channel in (10) as optimization variable 
ops = sdpsettings('solver','mosek','verbose',0);
beamformer = sdpvar(nTx,nUsers,'full','complex');
u = sdpvar(nUsers,nUsers,'full','complex');
t = sdpvar(nUsers+1,1); % also treat t as opt variable
myalpha = sdpvar(nUsers,1);

F = [u == channelpara*beamformer];

for iUser=1:nUsers
    F = [F,imag(u(iUser,iUser)) == 0]; % this is (10c)
%     F = [F,real(u(iUser,iUser)) >= t(iUser+1)*myalpha(iUser)]; % (10b) 
%     F = [F,cone([1 u(iUser,1:nUsers~=iUser)],myalpha(iUser))]; % (10b)
    F = [F,cone([1 u(iUser,1:nUsers~=iUser)],real(u(iUser,iUser))/t(iUser+1))]; % (10b)
end
F = [F,cone(beamformer(:),sqrt(min(power,PAeff*(1/t(1)-Po))))]; % (10d)

controllerisfeasibility = optimizer(F,[],ops,{channelpara,t},beamformer); % indicate channelpara and t are parameters of the optimizer
%}

tol = 0.005; % error tolerance

[optimalbeamformerEEmax,globalUB,globalLB] = Algorithm1_BRB(controllerisfeasibility,channel,sinrthreshold,power,Po,PAeff,tol);
computeEE(channel,optimalbeamformerEEmax,Po,PAeff)


semilogx(EE_SCA_seq*B/1e6*log2(exp(1)),'b')
hold on
semilogx(globalLB*B/1e6*log2(exp(1)),'r')
semilogx(globalUB*B/1e6*log2(exp(1)),'b')
ylabel('Energy Efficiency(Mb/J)')
xlabel('Iteration count')
legend('Alg. 3 (SCA)','Alg. 1 (BRB - Lower bound)','Alg. 1 (BRB - Upper bound)')
saveas(gcf,'../results/convergence.png')