无线通信代码搬运/复现系列(3) ： MIMO 广播信道的加权和速率最大化：使用脏纸编码和迫零方法

最新推荐文章于 2024-08-20 10:48:49 发布
最新推荐文章于 2024-08-20 10:48:49 发布
阅读量421 收藏 15
点赞数 10
分类专栏： 无线通信 文章标签： 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/White__River/article/details/141318763
版权

L. -N. Tran, M. Juntti, M. Bengtsson and B. Ottersten, “Weighted Sum Rate Maximization for MIMO Broadcast Channels Using Dirty Paper Coding and Zero-forcing Methods,” in IEEE Transactions on Communications, vol. 61, no. 6, pp. 2362-2373, June 2013,

我们考虑为最大化连续零强迫脏纸编码（SZF-DPC）的加权和速率（WSR）设计预编码器。对于这个问题，现有的预编码器设计通常假设有总功率约束（SPC），并依赖于奇异值分解（SVD）。基于 SVD 的设计被认为是最优的，但需要高复杂度。我们首先提出了一种低复杂度的最优预编码器设计，用于在 SPC 下的 SZF-DPC，采用 QR 分解。然后，我们提出了一种高效的数值算法，以找到满足每天线功率约束（PAPCs）的最优预编码器。为此，PAPCs 的预编码器设计被表述为一个具有协方差矩阵秩约束的优化问题。解决这个问题的一个著名方法是放宽秩约束并解决放宽后的问题。有趣的是，对于 SZF-DPC，我们能够证明秩放宽是紧的。因此，PAPCs 的最优预编码器设计是通过解决放宽后的问题来计算的，我们提出了一种定制的内点法，具有超线性收敛速率。还提出了两种次优预编码器设计，并与最优设计进行了比较。 我们还表明，所提出的数值方法适用于寻找块对角化方案的最优前置编码器。

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

代码解释

这段代码实现了一个障碍方法（Barrier Method）优化算法，用于求解具有总功率约束和逐个天线功率约束（Per-Antenna Power Constraint, PAPC）的MIMO（多输入多输出）系统中的预编码矩阵优化问题。其目的是最大化信道中的和速率（Sum Rate）。

1. 基本变量初始化
  • channel: 输入的信道矩阵。
  • Po: 各个天线的功率约束。
  • nTx: 发射天线的数量。
  • nUsers: 用户的数量。
  • effectivechannel: 每个用户的有效信道矩阵，经过处理得到的信道表示。
  • B: 存储用于正交化信道的矩阵。
  • A: 用于表示每个天线与用户相关联的矩阵，计算基于特定用户的复合信道。
2. 阶段I和阶段II

  • 阶段I：

    • 初始化变量，并通过计算牛顿步（Newton Step）来更新原始变量 Omega 和对偶变量 mymu。
    • 通过回溯线搜索确定适当的步长，确保更新的变量满足约束。
    • 继续迭代直到达到容忍的误差阈值。

  • 阶段II：

    • 在阶段I的基础上，通过进一步增加障碍函数的权重 t，继续进行优化。
    • 每次迭代都会更新 Omega 和 mymu，直到达到设定的最大值或达到收敛标准。
3. 辅助函数
  • computeresidualnorm: 计算残差的范数，用于判断算法的收敛性。
  • computeinfeasiblestartNewtonstep: 计算用于更新原始和对偶变量的牛顿步。
  • computeobjval: 计算目标函数的值，即和速率的负值。
  • computemodifiedobjval: 计算带有障碍项的目标函数的值，同时计算对应的梯度。

总结

这段代码实现了一个基于障碍方法的复杂优化算法，主要用于MIMO系统的预编码矩阵设计。在求解过程中，算法需要解决具有逐个天线功率约束的优化问题，通过逐步更新原始和对偶变量，最大化系统的和速率。

function [SumRate,Precoders,gap] = szfdpcpapcbarriermethod(channel,Po)
global A effectivechannel
nTx = size(channel{1},2);
nUsers = length(channel);

effectivechannel = cell(nUsers,1); % 用户的有效信道
B = cell(nUsers,1);
B{1} = eye(nTx);
A = cell(nUsers,nTx); % (33)中的矩阵A

h_bar=[];
Nr=zeros(1,nUsers);
for iUser=1:nUsers
    h_bar = [h_bar;channel{iUser}];
    if( iUser < nUsers)
        B{iUser+1} = null(h_bar); % B{iUser} 是从用户1到(iUser-1)的复合信道的正交基
    end
    effectivechannel{iUser} = channel{iUser}*B{iUser}; 
    for iTx=1:nTx
        A{iUser,iTx} = B{iUser}'*diag(1:nTx==iTx)*B{iUser}; % 计算 (33) 中的矩阵 A
    end    
end

initialOmega=cell(nUsers,1); % (31)中 \Omega 的初始值
for iUser=1:nUsers
    initialOmega{iUser} = eye(size(B{iUser},2)); % 设置为单位矩阵（也可以是其他值）
end

MAXITERS = 100;
ALPHA = 0.01; % 回溯线搜索参数
BETA = 0.5; % 回溯线搜索参数
NTTOL = 1e-5;
myphi = 50;
TOL = 1e-5;

%% 初始化
Omega=cell(nUsers,1);
Omega_next=cell(nUsers,1);
DeltaOmega=cell(nUsers,1);

t=1; % 算法2中的初始 t0 值
for iUser = 1:nUsers
    Omega{iUser} = initialOmega{iUser};
    DeltaOmega{iUser} = zeros(size(Omega{iUser}));
end
mymu = zeros(nTx,1);
mymunext = zeros(nTx,1);
iCenteringstep = 1;
iTotaliteration = 0;
sumratebest = [];

%% 阶段I
for iIteration = 1:MAXITERS
    iTotaliteration = iTotaliteration+1;
    modifiedobjval = computemodifiedobjval(t,Omega);
    sumrate(iTotaliteration) = computeobjval(Omega);
    sumratebest = [sumratebest min(sumrate)];
    sumratebarrier(iTotaliteration) = modifiedobjval;
    [DeltaOmega,Deltamu] = computeinfeasiblestartNewtonstep(t,Omega,mymu,Po); 
    residualnorm = computeresidualnorm(t,Omega,mymu,effectivechannel,A,Po);
    
    if (residualnorm < NTTOL),
        break; % 如果 residualnorm 很小，则退出 for 循环
    end;
    s=1;
    while(1) % 回溯线搜索
        feasibility_flag = 0;
        for iUser = 1:nUsers
            feasibility_flag = feasibility_flag+(min(real(eig(Omega{iUser}+s*DeltaOmega{iUser}))) <= 0);
            Omega_next{iUser} = Omega{iUser}+s*DeltaOmega{iUser};
        end
        mymunext = mymu+s*Deltamu;
        r_norm_next = computeresidualnorm(t,Omega_next,mymunext,effectivechannel,A,Po);
        
        minimum_flag = (r_norm_next > (1-ALPHA*s)*residualnorm);
        if(feasibility_flag==0) && (minimum_flag==0)
            break
        else
            s=s*BETA;
        end
        
    end
    % 更新原始变量和对偶变量
    mymu = mymunext;
    for iUser=1:nUsers
        Omega{iUser}=Omega_next{iUser};
    end
end
% 第一步居中步骤结束
iCenteringstep = iCenteringstep + 1;
t=t*myphi;

%% 阶段II 
while(t<1e8)
    % 内循环开始
    for iIteration=1:MAXITERS
        iTotaliteration=iTotaliteration+1;
        [modifiedobjval,Grad_barrier] = computemodifiedobjval(t,Omega);
        [DeltaOmega,Deltamu] = computeinfeasiblestartNewtonstep(t,Omega,mymu,Po);
        sumrate(iTotaliteration) = computeobjval(Omega);
        sumratebest=[sumratebest min(sumrate)];
        sumratebarrier(iTotaliteration) = modifiedobjval;
        
        % 计算牛顿减量 
        Newtondecrement = 0;
        for iUser=1:nUsers
            Newtondecrement = Newtondecrement + ...
                real(trace(Grad_barrier{iUser}' * DeltaOmega{iUser})); 
        end
        
        % 回溯线搜索
        if (abs(Newtondecrement) < NTTOL),
            % SumRate_barrier_centering(iCenteringstep) = modifiedobjval;
            break;
        end;
        
        s=1;
        while(1)
            feasibility_flag=0;
            for iUser=1:nUsers
                feasibility_flag=feasibility_flag+(min(real(eig(Omega{iUser}+s*DeltaOmega{iUser}))) <= 0);
            end
            
            for iUser=1:nUsers
                Omega_next{iUser}=Omega{iUser}+s*DeltaOmega{iUser};
            end
            mymunext=mymu+s*Deltamu;

            val_next = computemodifiedobjval(t,Omega_next);
            
            minimum_flag =(val_next > modifiedobjval + s*ALPHA*Newtondecrement);
            if(feasibility_flag==0) && (minimum_flag==0)
                break
            else
                s=s*BETA;
            end
            
        end
        % 更新原始变量和对偶变量
        mymu = mymunext;
        for iUser=1:nUsers
            Omega{iUser} = Omega_next{iUser};
        end
    end
    % 内循环结束
    iCenteringstep = iCenteringstep+1;
    t = t*myphi;
end

gap = sumratebest - sumratebest(end);
SumRate = -sumrate(end);
Precoders = Omega;

function residualnorm = computeresidualnorm(t,Omega,Deltamu,effectivechannel, A, Po)
nUsers=length(Omega);
Grad_orginal = cell(nUsers,1);
Grad_barrier = cell(nUsers,1);
rdual = cell(nUsers,1);
for iUser = 1:nUsers
    [Nr]=size(Omega{iUser},1);
    Grad_orginal{iUser} = -effectivechannel{iUser}'*effectivechannel{iUser}...
        *inv(eye(Nr)+Omega{iUser}*effectivechannel{iUser}'*effectivechannel{iUser});
    Grad_barrier{iUser} = t*Grad_orginal{iUser}-inv(Omega{iUser});
    rdual{iUser} = Grad_barrier{iUser};
end
nTx=size(Omega{1},2);
v=zeros(nTx,1);
for i = 1:nTx
    v(i) = Po(i);
    for iUser = 1:nUsers
        rdual{iUser} = rdual{iUser}+Deltamu(i)*A{iUser,i};
        v(i) = v(i)-real(trace(A{iUser,i}*Omega{iUser}));
    end
end
r = zeros(nUsers,1);
for iUser = 1:nUsers
    r(iUser) = r(iUser)+norm(rdual{iUser},'fro');
end
residualnorm = sum(r) + norm(v);



function [DeltaOmega,Deltamu] = computeinfeasiblestartNewtonstep(t,Omega,mymu,Po)
global effectivechannel A
nUsers=length(Omega);
Grad_orginal=cell(nUsers,1);
nTx=size(Omega{1},1);
myXi=cell(nUsers,nTx+1);
DeltaOmega=cell(nUsers,1);

for iUser=1:nUsers
    [Nr]=size(Omega{iUser},1);
    Grad_orginal{iUser}= -effectivechannel{iUser}'*effectivechannel{iUser}*inv(eye(Nr)+Omega{iUser}*effectivechannel{iUser}'*effectivechannel{iUser});
    
    P = t*Omega{iUser}*(-Grad_orginal{iUser});
    
    C = t*Omega{iUser}*(-Grad_orginal{iUser})*Omega{iUser}+Omega{iUser};
    for i=1:nTx
        C = C - mymu(i)*Omega{iUser}*A{iUser,i}*Omega{iUser};
    end
    myXi{iUser,1}=dlyap(P,-P'/t,C); % 求解 (42) 的第一个方程
    
    for i=1:nTx
        P=t*Omega{iUser}*(-Grad_orginal{iUser});
        C=-Omega{iUser}*A{iUser,i}*Omega{iUser};
        myXi{iUser,i+1}=dlyap(P,-P'/t,C); % 求解 (42) 的其余方程
    end
end

myPsi=zeros(nTx,nTx);
mypsi=zeros(nTx,1);

for i=1:nTx
    mypsi(i)=Po(i);
    for iUser=1=nUsers
        mypsi(i) = mypsi(i)-trace(A{iUser,i} * myXi{iUser,1})-trace(A{iUser,i}*Omega{iUser});
        for j=1:nTx
            myPsi(i,j)=myPsi(i,j)+trace(A{iUser,i} * myXi{iUser,j+1});
        end
    end
end

% 对偶变量的牛顿步
Deltamu=real(myPsi\mypsi);

% 原始变量的牛顿步

for iUser=1=nUsers
    DeltaOmega{iUser}=myXi{iUser,1};
    for i=1:nTx
        DeltaOmega{iUser}=DeltaOmega{iUser}+Deltamu(i)*myXi{iUser,i+1};
        DeltaOmega{iUser} = (DeltaOmega{iUser}+DeltaOmega{iUser}')/2; % 消除数值误差
    end
end

%%
function [val]=computeobjval(Omega)
global effectivechannel
nUsers = length(Omega);
val = 0;
for iUser = 1:nUsers
    Nr = size(effectivechannel{iUser},1);
    val = val+log(det(eye(Nr) + effectivechannel{iUser} ...
        * Omega{iUser} *  effectivechannel{iUser}'));
end
val = -real(val);

%% 
function [ modifiedobjval ,Grad_barrier] = computemodifiedobjval(t,Omega)
global effectivechannel
nUsers = length(Omega);
Grad_orginal = cell(nUsers,1);
Grad_barrier = cell(nUsers,1);

modifiedobjval = 0;
for iUser = 1:nUsers
    [Nr,Nr_bar] = size(effectivechannel{iUser});
    modifiedobjval = modifiedobjval+t*log(det(eye(Nr) + effectivechannel{iUser} * Omega{iUser} *...
        effectivechannel{iUser}'))+log(det(Omega{iUser}));
    Grad_orginal{iUser} = -effectivechannel{iUser}'*effectivechannel{iUser}*...
        inv(eye(Nr_bar)+Omega{iUser}*effectivechannel{iUser}'*effectivechannel{iUser});
    Grad_barrier{iUser} = t*Grad_orginal{iUser}-inv(Omega{iUser});
end
modifiedobjval = -real(modifiedobjval);


clear all
clc
rng(1)
nTx = 6; % number of users
nRxarray = [2 2 2];   %  number of receive antennas for each user, i.e.,
%                         %  nRxarray(i) is the number of rx. antennas at the ith user
                        
% nRxarray = [2 2 ];   %  number of receive antennas for each user, i.e.,
nUsers = length(nRxarray); % number of users


P = 10; % sum power, dB scale
P = 10.^(P./10);
Po = P/nTx*ones(nTx,1); % maximum power per antenna 

channel = cell(nUsers,1); % channel matrices of users

tic
for iUser=1:nUsers
    channel{iUser} = sqrt(1/2)*(randn(nRxarray(iUser),nTx) + 1i*randn(nRxarray(iUser),nTx));    
end
[sumrate,precoder,gap] = szfdpcpapcbarriermethod(channel,Po);
toc
semilogy(gap)
xlabel('Iteration')
ylabel('Residual error')
saveas(gcf,'../results/convergence.png')
%% YALMIP code for solving (31)
effectivechannel = cell(nUsers,1);
effectivechannel{1} = channel{1};
B = cell(nUsers,1);
B{1} = eye(nTx);
h_bar = [];
A = cell(nUsers,nTx); % matrix A in (33)
for iUser=1:nUsers
    h_bar = [h_bar;channel{iUser}];
    if( iUser < nUsers)
        B{iUser+1} = null(h_bar);
    end
    effectivechannel{iUser} = channel{iUser}*B{iUser};
    for iTx=1:nTx
        A{iUser,iTx}=B{iUser}'*diag(1:nTx==iTx)*B{iUser}; % compute matrix A in (33)
    end    
end


X=cell(nUsers,1);
F=[];
for iUser = 1:nUsers
    X{iUser} = sdpvar(nTx-sum(nRxarray(1:iUser-1)),nTx-sum(nRxarray(1:iUser-1)),'hermitian','complex');
    F=[F,X{iUser}>=0];
end
obj=0;
for iUser=1:nUsers
    obj=obj+logdet(eye(nRxarray(iUser)) + effectivechannel{iUser}*X{iUser}*effectivechannel{iUser}');
end

for iTx=1:nTx
    y = 0;
    for iUser=1:nUsers
        y = y + real(trace(A{iUser,iTx}*X{iUser}));
    end
    F=[F,y <= Po(iTx)];
end

ops = sdpsettings('solver','sdpt3','verbose',0);

tic

sol = optimize(F,-obj,ops);

sumrateyalmip = real(double(obj));
precodersoftware = cell(nUsers,1);
for iUser = 1:nUsers
    precodersoftware{iUser} = B{iUser}*double(X{iUser})*B{iUser}';
end
toc

norm(sumrate-sumrateyalmip)

在这里插入图片描述

sjtuSecondSilly
关注
  • 10
    点赞
  • 踩
  • 15
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
mimo-ofdm信道matlab代码,MIMO-OFDM无线通信及MATLAB实现代码
weixin_29868343的博客
03-17 1393
文件名大小更新时间MIMO-OFDM无线通信及MATLAB实现代码02019-05-13MIMO-OFDM无线通信及MATLAB实现代码\add_CFO.m5512019-01-19MIMO-OFDM无线通信及MATLAB实现代码\add_CP.m2952019-01-19MIMO-OFDM无线通信及MATLAB实现代码\add_pilot.m7892019-01-19MIMO-OFDM无线通信及...
通信与网络中的新一代无线标准：IEEE 802.11ac和LTE-Advanced
10-20
新一代无线通信标准IEEE 802.11ac和LTE-Advanced的出现，无线网络的传输速度得到了显著提升。这两个标准分别在无线局域网（WLAN）和蜂窝通信领域占据领先地位，它们通过技术创新实现了更高的数据吞吐能力，以满足...
无线发射源代码
03-08
由单片机控制无线模块发射数据，经测试开阔地带可达到450米，数据传输可靠
无线共享代码
weixin_30915951的博客
01-26 111
@echo off :: BatchGotAdmin :------------------------------------- REM --> Check for permissions >nul 2>&1 "%SYSTEMROOT%\system32\cacls.exe" "%SYSTEMROOT%\system32\config\system" ...
一段WiFi代码
zhilin_tang的专栏
04-22 5540
#include #include #include #include #include #include #include #include "hardware_legacy/wifi.h" #include "libwpa_client/wpa_ctrl.h" #define LOG_TAG "WifiHW" #include "cutils/log.h"
基于MATLAB的MIMO预编码设计：优化迫零算法（附完整代码与分析）
forest_LL的博客
02-24 5521
图中“RFChain” 全称为RadioFrequencyChain，代表射频链路。此MIMO预编码包含了基带预编码W（改变幅度和相位）和射频预编码F（改变相位）。用户k端收到的信号为：上式中代表基站到用户的信道矩阵，F代表射频预编码，W代表基带预编码，s代表发射信号向量，代表加性高斯白噪声。发射流的最大数量为K，则：用户k端受到的信干噪比（SINR）为：上式子中P代表基站的发射功率。其他参数的意义上面已解释。
传统多用户MISO系统中的加权和速率WSR最大化问题，（优化AP的beamformingW）
qq_45542321的博客
01-16 2143
传统多用户MISO系统中的加权和速率WSR最大化问题，（优化AP的beamformingW）
MIMO-OFDM无线通信技术及MATLAB实现（1）无线信道：传播和衰落
tybAly的博客
05-03 4458
目录 大尺度衰落 一般路径损耗模型 Okumura/Hata模型 IEEE 802.16d模型 小尺度衰落 小尺度衰落的参数 时间色散衰落与频率色散衰落 衰落信道的统计特性和产生 大尺度衰落 一般路径损耗模型 function PL = PL_free(fc,dist,Gt,Gr) %自由空间路径损耗模型 lamda=3e8/fc; tmp=lamda./(4*pi*dist); if nargin&gt......
多用户MIMO干扰网络下的加权和速率优化算法_MIMO_加权和速率最大化_干扰网络
09-23
多用户MIMO干扰网络下的加权和速率优化算法
MIMO-OFDM无线通信技术及MATLAB实现 无线信道：传播和衰落.zip
07-19
MIMO技术通过在发射端使用多个天线发送信号，同时在接收端使用同样数量或更多天线接收，从而显著提高数据传输速率和系统容量。OFDM技术则是将高速数据流分割成多个低速子载波，并在每个子载波上进行调制，以克服多径...
本科毕设的全部代码（在Matlab平台实现）-主要研究的是MIMO通信系统中的3大部分：空时编码、系统容量、信道估计
06-09
主要研究的是MIMO通信系统中的3大部分：空时编码、系统容量、信道估计 代码很多、很杂，留给后来人学习、参考。 各个文件的相应使用方法在具体的文件夹里都有 MIMO技术是一项十分重要的技术，能在不增加带宽的前提下...
mimo信道容量及和速率的仿真代码
11-18
在无线通信领域，MIMO（Multiple-Input Multiple-Output）技术是一种重要的传输方式，通过利用空间多样性和信号处理技术来显著提升系统容量和数据传输速率。本资源提供了一个基于MATLAB的MIMO信道容量及和速率的仿真...
如何在笔记本上建立无线代码
yy2065913784的专栏
10-19 177
下面我来教大家怎么在笔记本上建立无线代码，首页把以下w
代码随想录算法训练营第十七天（二叉树 五）
b1ankwall的博客
08-16 779
参加代码随想录第43期算法训练营day17|654.最大二叉树，617.合并二叉树，700.二叉搜索树中的搜索，98.验证二叉搜索树
【Python】生成二维迷宫的算法
HarryXYC的博客
08-16 332
基于Pygame和PIL的迷宫生成器
【算法】马踏棋盘（骑士周游）问题回溯算法实现以及使用贪心算法优化
2302_78914800的博客
08-18 538
将马儿随机放在国际象棋的 8*8 棋盘的某个方格中，马儿按照“马走日”进行移动，要求每个方格只进入一次，走遍棋盘上全部 64 个方格。注意：马儿不同的走法会得到不同的结果，效率也会有影响。
Leetcode 219. 存在重复元素 II
最新发布
铁马冰河入梦来
08-20 322
哈希表解决存在重复元素二问题
noip2001数的划分
no_play_no_games的博客
08-18 188
将整数n分成k份，且每份不能为空，任意两份不能相同(不考虑顺序)。例如：n=7，k=3，下面三种分法被认为是相同的。1，1，5;1，5，1;5，1，1;问有多少种不同的分法。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值