无线通信代码搬运/复现系列(2) ： 可实现速率优化对于具有可重构智能表面的 MIMO 系统

N. S. Perović, L. -N. Tran, M. Di Renzo and M. F. Flanagan, “Achievable Rate Optimization for MIMO Systems With Reconfigurable Intelligent Surfaces,” in IEEE Transactions on Wireless Communications, vol. 20, no. 6, pp. 3865-3882, June 2021, doi: 10.1109/TWC.2021.3054121.

可重构智能表面（RIS）代表了一种新技术，可以塑造无线网络中的无线电波传播，并提供了多种可能的性能和实施收益。受此启发，我们研究了配备 RIS 的多流多输入多输出（MIMO）系统的可实现速率优化，并制定了传输信号的协方差矩阵和 RIS 元素的联合优化问题。为了解决这个问题，我们提出了一种基于投影梯度方法（PGM）的迭代优化算法。我们推导出保证所提算法收敛的步长，并定义了一个回溯线搜索来提高其收敛速度。此外，我们引入了间接和直接链路的总自由空间路径损耗（FSPL）比作为所考虑通信系统中 RIS 适用性的一阶度量。模拟结果显示，所提出的 PGM 实现了与最先进基准方案相同的可实现速率，但计算复杂度显著降低。 此外，我们证明 RIS 应用特别适合提高室内环境中的可实现速率，因为即使少量的 RIS 元素也可以提供显著的可实现速率增益。

信道矩阵代码

function [Hdir,H1,H2] = chan_mat_RIS_surf_univ_new(Nt,Nr,Nris,lt,lr,D,no_mat,K,f,dist_ris,varargin)

lambda = 3e8/f;     % 波长
dt = lambda/2;      % 发射天线间距
dr = lambda/2;      % 接收天线间距
dris = lambda/2;    % RIS元件间距
k = 2*pi/lambda;    % 波数

% 几何布局
% x, y 和 z 轴 
% 发射天线阵列
tx_arr(1,:) = zeros(1,Nt); 
tx_arr(2,:) = (sort(0:Nt-1,'descend')-(Nt-1)/2)*dt+lt; 
tx_arr(3,:) = zeros(1,Nt); 
% 接收天线阵列
rx_arr(1,:) = D*ones(1,Nr);
rx_arr(2,:) = (sort(0:Nr-1,'descend')-(Nr-1)/2)*dr+lr; 
rx_arr(3,:) = zeros(1,Nr);
% RIS 
center = [dist_ris 0]; % RIS 中心位置 
N1 = sqrt(Nris);
N2 = N1;                                    % RIS 在两个维度 N1 和 N2 中的元素数量
ris_pos = RISPosition(N1,N2,dris,center);   % RIS 元素的坐标 
a = repmat(ris_pos{1},N1,1);                % 将RIS元素放置在适当的坐标
ris_arr(1,:) = a(:)';        
ris_arr(2,:) = zeros(1,Nris);
ris_arr(3,:) = repmat(ris_pos{2},1,N2); 

if isempty(varargin)                        % 加载直连链路的自由空间路径损耗 (FSPL)
    alpha = 2;
else
    alpha = varargin{1};
end

% 直接 TX-RX 路径/信道矩阵
for i1 = 1:Nr                                        % 计算 TX 和 RX 天线之间的距离                                          
    for j1 = 1:Nt
        d(i1,j1) = norm(rx_arr(:,i1)-tx_arr(:,j1));
    end 
end
Hdir_los = exp(-1i*k*d);                             % 直接链路，视距 (LOS) 矩阵指数 
tx_rx_dist = sqrt(D^2+(lt-lr)^2);                    % TX-RX 距离   
FSPL_dir = (lambda/(4*pi))^2/tx_rx_dist^alpha(1);    % 直接链路的自由空间路径损耗 (FSPL) 的倒数 
Hdir = Rician(Hdir_los,sqrt(FSPL_dir),no_mat,K);     % 直接链路信道矩阵           

% 间接路径 (TX-RIS-RX)
for l1 = 1:Nris                                      % 计算 RIS 元素和 RX 天线之间的距离                                                            
    for r1 = 1:Nr  
        d2(r1,l1) = norm(rx_arr(:,r1)-ris_arr(:,l1)); 
    end
    for t1 = 1:Nt                                    % 计算 RIS 元素和 TX 天线之间的距离                  
        d1(l1,t1) = norm(tx_arr(:,t1)-ris_arr(:,l1));   
    end
end

tx_ris_dist = sqrt(dist_ris^2+lt^2);                 % TX-RIS 距离
ris_rx_dist = sqrt((D-dist_ris)^2+lr^2);             % RIS-RX 距离   

FSPLindir = lambda^4/(256*pi^2)*...                  % 间接链路的自由空间路径损耗 (FSPL) 的倒数 
           ((lt/tx_ris_dist+lr/ris_rx_dist)^2)*...
           1/(tx_ris_dist*ris_rx_dist)^2;

% TX-RIS 信道矩阵
H1_los = exp(-1i*k*d1);                             % TX-RIS 链路，视距 (LOS) 矩阵指数  
FSPL_1 = sqrt(FSPLindir);                           % 间接链路的 FSPL 嵌入在 TX-RIS 信道矩阵中 
H1 = Rician(H1_los,FSPL_1,no_mat,K);

% RIS-RX 信道矩阵
H2_los = exp(-1i*k*d2);                             % RIS-RX 链路，视距 (LOS) 矩阵指数 
FSPL_2 = 1;
H2 = Rician(H2_los,FSPL_2,no_mat,K);
end

function pos = RISPosition(N1,N2,dist,center)                 % 确定RIS元素的位置
d1 = (0:N1-1)-(N1-1)/2;
d2 = (0:N2-1)-(N2-1)/2;
pos{1} = center(1)+d1*dist;
pos{2} = center(2)+d2*dist;
end 

function Hout = Rician(Hlos,FSPL,no_mat,K)                     % 创建 Rician 信道矩阵
Hlos = repmat(Hlos,no_mat,1);
Hnlos = sqrt(1/2)*(randn(size(Hlos))+1i*randn(size(Hlos)));
Htot = FSPL/sqrt(K+1)*(Hlos*sqrt(K)+Hnlos);
dim = size(Hlos,1)/no_mat;
for ind = 1:no_mat
   Hout{ind} = Htot((ind-1)*dim+1:ind*dim,:); 
end
end

代码实现了一个基于投影梯度法（PGM）的迭代优化算法，用于联合优化发射协方差矩阵和再生智能表面（RIS）相移，从而最大化通信系统中的可达速率。代码还包括了与协方差矩阵和RIS相移相关的梯度计算、速率计算、单位圆投影以及水填充算法。

% 实现PGM迭代优化的函数

function [Cout,iter_time] = PGM_opt(Pt,Hdir,H1,H2,maxIter,Qinit,myomegaunitball,c)

% 加载初始协方差矩阵和RIS相移  
myomegaunitcirc = myomegaunitball.';
Q = Qinit; 

% 线搜索参数 
delta = 1e-5;       
rho = 0.5;

iIter = 0;         
stepsize = 10000;   % 初始步长 

% 初始化可达速率
Cout = [RIScap(Hdir,H2,H1,myomegaunitcirc,Q)]; % 初始可达速率 
Cprev = Cout;

iter_time = 0;
tic
while iIter<maxIter
    iIter = iIter+1;
    q_Q = gracov(Hdir,H2,H1,myomegaunitcirc,Q);     % 相对于Q的梯度
    g_RIS = gradRIS(Hdir,H2,H1,myomegaunitcirc,Q);  % 相对于RIS的梯度
   
    for iLineSearch = 0:30

        % 更新协方差 (Q) 矩阵, 算法1第3行 
        Qnew = Q + stepsize*q_Q;
        Qnew = cov_mat_proj_modified(Qnew,Pt*c^2);
        
        % 更新RIS相移, 算法1第4行
        yunitcirc = myomegaunitcirc + stepsize*g_RIS;
        myomegaunitcircnext = projectontounitcircle(yunitcirc)/c;
        
        % 新的可达速率 
        Cnew = RIScap(Hdir,H2,H1,myomegaunitcircnext,Qnew);
        
        % 线搜索过程
        if (((Cnew-Cprev) >= delta*(norm(Qnew-Q)^2+ ...
                norm(myomegaunitcircnext-myomegaunitcirc)^2)) ...
            || (stepsize<1e-4) )
            % 如果步长满足条件(35c)或步长过小：    
            % 获得的Q矩阵和RIS相移
            myomegaunitcirc = myomegaunitcircnext;
            Q = Qnew;
            Cprev = Cnew;
            break
        else 
            % 缩小步长 
            stepsize=stepsize*rho;
        end
    end
    % 新的可达速率
    Cout = [Cout Cnew];
    % 每次迭代的执行时间
    iter_time = [iter_time toc];
end
end

% 相对于Q的梯度计算
function y = gracov(Hdir,H2,H1,myomega,Q)
Z = Hdir+H2*diag(myomega)*H1;
Nr = size(H2,1);
y = Z'*inv(eye(Nr)+Z*Q*Z')*Z;
end

% 相对于RIS的梯度计算
function y = gradRIS(Hdir,H2,H1,myomega,Q)
Z = Hdir+H2*diag(myomega)*H1;
Nr = size(H2,1);
y = diag(H2'*inv(eye(Nr)+Z*Q*Z')*Z*Q*H1');
end

% 可达速率计算
function y = RIScap(Hdir,H2,H1,myomega,Q)
Z = Hdir+H2*diag(myomega)*H1;
Nr = size(H2,1);
y = real(log(det(eye(Nr)+Z*Q*Z')))/log(2);
end

% RIS投影
function y = projectontounitcircle(x)
y = x./abs(x);
end

% 协方差矩阵投影
function Qnew = cov_mat_proj_modified(Qold,Pt)
    [U,D] = eig(Qold);
    Dnew = water_fill(Pt,real(diag(D))).';
    Qnew = U*diag(Dnew)*U';
end 

% 水填充算法
function vect_out = water_fill(Pt,vect_in)
    vect_in = vect_in.';
    [sort_val,sort_idx] = sort(vect_in,'descend');
    for n = length(vect_in):-1:1
        water_level = (sum(sort_val(1:n))-Pt)/n;
        di = sort_val(1:n)-water_level;
        if di(:)>=0
            break
        end   
    end
    vect_out = zeros(1,length(vect_in));
    vect_out(sort_idx(1:n)) = di;
end 

主程序

clear; close all; clc;

Nt = 8;                 % 发射天线的数量
Nr = 4;                 % 接收天线的数量
Nris = 15^2;            % RIS 元素的数量
       
K = 1;                  % Rician 因子    
D = 500;                % 发射端到接收端的距离
dist_ris = 40;          % RIS 到发射端的距离
f = 2e9;                % 频率

lt = 20;                % 发射端的位置 
lr = 100;               % 接收端的位置 
Pt = 1;                 % 发射功率，单位：瓦特
N0 = -120;              % 噪声功率，单位：dB
SNR = db2pow(-N0);      % 信噪比 (SNR)            
no_mat = 10;            % 信道实现的数量 

no_iter = 500;          % 迭代次数 
alpha_dir = 3;          % 直接链路的自由空间路径损耗 (FSPL) 指数

% 生成信道矩阵（注意：此代码仅适用于 Nris 为平方数的情况）
[Hdirt,H1t,H2t] = chan_mat_RIS_surf_univ_new(Nt,Nr,Nris,lt,lr,D,no_mat,K,f,dist_ris,alpha_dir);

Cpgm = zeros(1,no_iter+1);

for i = 1:no_mat 
    Hdir = Hdirt{i}; H1 = H1t{i}; H2 = H2t{i};                  
    
    % 缩放因子
    c = sqrt(norm(Hdir)/norm(H2*H1))*max(sqrt(Pt),1)/sqrt(Pt)*10;     
    
    % 初始的 Q 矩阵和 RIS 相移    
    Qinit = eye(Nt)*(Pt/Nt);
    omega_init = ones(1,Nris);
    
    % PGM 迭代优化 
    [dCpgm,~] = PGM_opt(Pt,Hdir*sqrt(SNR)/c,H1*sqrt(SNR),H2,no_iter,Qinit*c^2,omega_init/c,c);
    Cpgm = Cpgm+dCpgm;
    
end

% 绘制可达速率与迭代次数的图像
semilogx(1:length(Cpgm),Cpgm/no_mat,'r','DisplayName','PGM');
xlabel('迭代次数'); ylabel('可达速率 [bit/s/Hz]');
xlim([0 no_iter]); 
legend('show','Location','SouthEast');
print('../results/Achievable_Rate', '-dpdf')