ToB企服应用市场:ToB评测及商务社交产业平台

标题: 智能反射面（IRS）在无线通讯安全领域应用的论文复现 [打印本页]
作者: 盛世宏图    时间: 2024-6-11 19:06
标题: 智能反射面（IRS）在无线通讯安全领域应用的论文复现
弁言

Zhang Rui老师的将IRS引入无线通讯安全的论文《Secure Wireless Communication via Intelligent Reflecting Surface》有较高的引用量，在此给出要论文的复现及代码。
主要问题

该论文的目的是引入IRS并联合优化基站的主动式波束和IRS的被动式波束，使得抑制窃听者信噪比的同时最大化合法用户处的信噪比。其场景如下：
      
       体系模型图        因此可以构造出以下的优化问题：
      
       优化问题        即在基站发射功率的约束下，优化基站和IRS的波束使得保密速率最大化。
给定IRS相位矩阵时，优化基站波束

可简单地将求绝对值的平方进行简单展开，令

将对数相减变更为真数相除，对数是单调递增函数，因此最大化对数，即是最大化真数即可。因此，可简化为以下的问题： 

该结构可以参考文献[1]直接给出解的形式如下： 

 其中
对应于矩阵
的最大特性值对于的归一化特性向量。
给定基站波束，优化IRS相位矩阵

该部分推导较为复杂，可以详细阅读论文，如果有不懂的地方，可以评论或私信互换。主要是利用了分式规划将其转化为一个半正定松弛问题，求解该问题然后利用高斯随机化过程进行求解，转化后的问题如下所示：

 仿真复现

仿真参数设置

  1. clc;clear;
  3. epsilon = 1e-3; % 收敛停止条件
  5. % 天线数
  6. M = 4; % AP天线数
  7. Nx = 8;
  8. Nz = 8;
  9. N = Nx*Nz; % IRS单元个数
  11. % 用户位置
  12. APloc = [0;0]; % AP位置
  13. userloc = [150;0]; % user位置
  14. edloc= [145;0]; % 窃听者位置
  15. IRSloc = [145;5]; % IRS位置
  17. C0 = db2pow(-30); % 参考距离时的路损
  18. D0 = 1; % 参考距离
  19. sigmaK2 = db2pow(-80); % 噪声功率，-80dBm
  21. L = @(d, alpha)C0*(d/D0)^(-alpha); % 路损模型
  23. % 路损参数
  24. alpha_AU = 3;
  25. alpha_AE = 3;
  26. alpha_AI = 2.2;
  27. alpha_IU = 3;
  28. alpha_IE = 3;
  30. % 莱斯因子
  31. K_AU = db2pow(1);
  32. K_AE = db2pow(1);
  34. K_AI = db2pow(1);
  35. K_IU = db2pow(1);
  36. K_IE = db2pow(1);
  39. R = 1000; % 信道实现数
  41. P_AP = db2pow(25); % 发射功率15dBm
复制代码
 产生信道

论文中说明信道都为独立的莱斯信道，论文中有些信道考虑的是具有空间相关性的莱斯信道，需要在NLoS部分前后乘以一个相关系数矩阵，具体内容可以参考论文[1]，为简化，在此没有考虑相关系数矩阵，则可以产生如下信道：
  1. dAE = norm(APloc-edloc);
  2. hAE = sqrt(L(dAE,alpha_AE)/sigmaK2)*(sqrt(1/(1+K_AE))*ones(M,1)'+sqrt(K_AE/(1+K_AE))*(randn(1,M)+1i*randn(1,M))/sqrt(2));
  4. dAU = norm(APloc-userloc);
  5. hAU = sqrt(L(dAU,alpha_AU)/sigmaK2)*(sqrt(1/(1+K_AU))*ones(M,1)'+sqrt(K_AU/(1+K_AU))*(randn(1,M)+1i*randn(1,M))/sqrt(2));
  7. dAI = norm(APloc-IRSloc);
  8. thetaIRS = atan(145/5);phi = 0; thetaAP = atan(5/145);
  9. HAI = sqrt(L(dAI,alpha_AI)/sigmaK2)*(sqrt(K_AI/(1+K_AI))*URA_sv(thetaIRS, phi,Nx,Nz)*ULA_sv(thetaAP,M)'+sqrt(1/(1+K_AI))*(randn(N,M)+1i*randn(N,M))/sqrt(2));
  11. dIU = norm(IRSloc-userloc);
  12. thetaIRS = -pi/4; phi = 0;
  13. hIU = sqrt(L(dIU,alpha_IU))*(sqrt(K_IU/(1+K_IU))*URA_sv(thetaIRS, phi,Nx,Nz)'+sqrt(1/(1+K_IU))*(randn(1,N)+1i*randn(1,N))/sqrt(2));
  15. dIE = norm(IRSloc-edloc);
  16. thetaIRS = 0; phi = 0;
  17. hIE = sqrt(L(dIE,alpha_IE))*(sqrt(K_IE/(1+K_IE))*URA_sv(thetaIRS, phi,Nx,Nz)'+sqrt(1/(1+K_IE))*(randn(1,N)+1i*randn(1,N))/sqrt(2));
复制代码
 迭代优化

  1. q = 2*pi*rand(1,N); % 随机初始化IRS的相位
  2. Q = diag(exp(1i*q));
  4. % 给定q优化W
  5. A = (hIU*Q*HAI+hAU)'*(hIU*Q*HAI+hAU); % 公式(9)
  6. B = (hIE*Q*HAI+hAE)'*(hIE*Q*HAI+hAE); % 公式(10)
  8. I_M = eye(M);
  9. C = (B+1/P_AP*I_M)\(A+1/P_AP*I_M);
  10. [V,D] = eig(C); % 特征值分解
  11. [d,ind] = sort(diag(D));
  12. u_max = V(:,ind(end))/norm(V(:,ind(end)));
  13. w_opt = sqrt(P_AP) * u_max;
  15. % 给定W优化q
  16. hU = conj(hAU)*conj(w_opt)*(w_opt.')*(hAU.');
  17. hE = conj(hAE)*conj(w_opt)*(w_opt.')*(hAE.');
  18. GU = [diag(conj(hIU))*conj(HAI)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU)  diag(conj(hIU))*conj(HAI)*conj(w_opt)*(w_opt.')*(hAU.');...
  19.     conj(hAU)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU) 0];
  20. GE = [diag(conj(hIE))*conj(HAI)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE)  diag(conj(hIE))*conj(HAI)*conj(w_opt)*(w_opt.')*(hAE.');...
  21.     conj(hAE)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE) 0]; % 公式(18)
  22. q = SDR(hU,hE,GU,GE,N);
  24. Q = diag(q);
  25. R = max(0, log2(1+abs((hIU*Q*HAI+hAU)*w_opt)^2)- log2(1+abs((hIE*Q*HAI+hAE)*w_opt)^2))
复制代码
 URA导向矢量函数

  1. function ura_sv = URA_sv(theta, phi,Nx,Ny)
  2.     m = 0:Nx-1;
  3.     a_az = exp(1i*pi*m*sin(theta)*cos(phi)).';
  4.     n = 0:Ny-1;
  5.     a_el = exp(1i*pi*n*sin(phi)).';
  6.     ura_sv = kron(a_az,a_el);
  7. end
复制代码
ULA导向矢量函数

  1. function ula_sv = ULA_sv(theta, M)
  2.     m = 0:M-1;
  3.     ula_sv = exp(1i*pi*m*sin(theta)).';
  4. end
复制代码
半正定松弛优化函数

SDR求解问题(22a)
  1. function [q,cvx_optval] = SDR(hU,hE,GU,GE,N)
  2.     L = 1000; % 高斯随随机化次数
  3.     cvx_begin sdp quiet
  4.         variable X(N+1,N+1) hermitian
  5.         variable mu1(1,1)
  6.         maximize(real(trace(GU*X)+mu1*(hU+1)))
  7.         subject to
  8.             real(trace(GE*X))+mu1*(hE+1)==1;
  9.             for i=1:N+1
  10.                 En = zeros(N+1,N+1);
  11.                 En(i,i)=1;
  12.                 real(trace(En*X)) == mu1;
  13.             end
  14.             X ==  hermitian_semidefinite(N+1);
  15.             mu1 >= 0;
  16.     cvx_end
  18.    % 高斯随机化过程
  19.     %% method 1
  20.     max_F = 0;
  21.     max_q = 0;
  22.     S = X / mu1;
  23.     [U, Sigma] = eig(S);
  24.     for l = 1 : L
  25.         r = sqrt(2) / 2 * (randn(N+1, 1) + 1j * randn(N+1, 1));
  26.         q = U * Sigma^(0.5) * r;
  27.         if q' * GU * q > max_F
  28.             max_q = q;
  29.             max_F = q' * GU * q;
  30.         end
  31.     end
  32.    
  33.     q = exp(1j * angle(max_q / max_q(end)));
  34.     q = q(1 : N);
  35. end
复制代码
以上步伐是给定发射功率的单点优化步伐，仿真随着发射功率变革的完备步伐以及对比算法如下：

  1. clc;clear;
  3. epsilon = 1e-3; % 收敛停止条件
  5. % 天线数
  6. M = 4; % AP天线数
  7. Nx = 8;
  8. Nz = 8;
  9. N = Nx*Nz; % IRS单元个数
  11. % 用户位置
  12. APloc = [0;0]; % AP位置
  13. userloc = [150;0]; % user位置
  14. edloc= [145;0]; % 窃听者位置
  15. IRSloc = [145;5]; % IRS位置
  17. C0 = db2pow(-30); % 参考距离时的路损
  18. D0 = 1; % 参考距离
  19. sigmaK2 = db2pow(-80); % 噪声功率，-80dBm
  21. L = @(d, alpha)C0*(d/D0)^(-alpha); % 路损模型
  23. % 路损参数
  24. alpha_AU = 3;
  25. alpha_AE = 3;
  26. alpha_AI = 2.2;
  27. % alpha_IU = 3;
  28. alpha_IU = 2.3;
  29. % alpha_IE = 3;
  30. alpha_IE = 2.5;
  32. % 莱斯因子
  33. K_AU = db2pow(1);
  34. K_AE = db2pow(1);
  36. % K_AI = db2pow(1);
  37. K_AI = db2pow(10);
  38. % K_IU = db2pow(1);
  39. K_IU = db2pow(10);
  40. % K_IE = db2pow(1);
  41. K_IE = db2pow(10);
  44. P = db2pow(-5:5:25); % 发射功率15dBm
  45. frame = 10;
  46. maxIter = 20;
  47. RSDR = zeros(length(P),1);
  48. RMRT = zeros(length(P),1);
  49. RWIRS = zeros(length(P),1);
  50. RUB = zeros(length(P),1);
  51. for p=1:length(P)
  52.     p
  53.     P_AP = P(p);
  54.     for fr = 1:frame
  55.         dAE = norm(APloc-edloc);
  56.         hAE = sqrt(L(dAE,alpha_AE)/sigmaK2)*(sqrt(K_AE/(1+K_AE))*ones(M,1)'+sqrt(K_AE/(1+K_AE))*(randn(1,M)+1i*randn(1,M))/sqrt(2));
  57.         
  58.         dAU = norm(APloc-userloc);
  59.         hAU = sqrt(L(dAU,alpha_AU)/sigmaK2)*(sqrt(K_AU/(1+K_AU))*ones(M,1)'+sqrt(K_AU/(1+K_AU))*(randn(1,M)+1i*randn(1,M))/sqrt(2));
  60.         
  61.         dAI = norm(APloc-IRSloc);
  62.         thetaIRS = atan(145/5);phi = 0; thetaAP = atan(5/145);
  63.         HAI = sqrt(L(dAI,alpha_AI)/sigmaK2)*(sqrt(K_AI/(1+K_AI))*URA_sv(thetaIRS, phi,Nx,Nz)*ULA_sv(thetaAP,M)'+sqrt(K_AI/(1+K_AI))*(randn(N,M)+1i*randn(N,M))/sqrt(2));
  64.         
  65.         dIU = norm(IRSloc-userloc);
  66.         thetaIRS = -pi/4; phi = 0;
  67.         hIU = sqrt(L(dIU,alpha_IU))*(sqrt(K_IU/(1+K_IU))*URA_sv(thetaIRS, phi,Nx,Nz)'+sqrt(K_IU/(1+K_IU))*(randn(1,N)+1i*randn(1,N))/sqrt(2));
  68.         
  69.         dIE = norm(IRSloc-edloc);
  70.         thetaIRS = 0; phi = 0;
  71.         hIE = sqrt(L(dIE,alpha_IE))*(sqrt(K_IE/(1+K_IE))*URA_sv(thetaIRS, phi,Nx,Nz)'+sqrt(K_IE/(1+K_IE))*(randn(1,N)+1i*randn(1,N))/sqrt(2));
  73.         q = 2*pi*rand(1,N); % 随机初始化IRS的相位
  74.         Q = diag(exp(1i*q));
  75.         
  76.         R_old = 0;
  77.         R_new = 10;
  78.         count = 0;
  79.         while(abs(R_old-R_new)/R_new > epsilon && count < maxIter)
  80.             count = count + 1;
  81.             % 给定q优化W
  82.             A = (hIU*Q*HAI+hAU)'*(hIU*Q*HAI+hAU); % 公式(9)
  83.             B = (hIE*Q*HAI+hAE)'*(hIE*Q*HAI+hAE); % 公式(10)
  84.             
  85.             I_M = eye(M);
  86.             C = (B+1/P_AP*I_M)\(A+1/P_AP*I_M);
  87.             [V,D] = eig(C); % 特征值分解
  88.             [d,ind] = sort(diag(D));
  89.             u_max = V(:,ind(end))/norm(V(:,ind(end)));
  90.             w_opt = sqrt(P_AP) * u_max;
  91.             
  92.             % 给定W优化q, SDR
  93.             hU = conj(hAU)*conj(w_opt)*(w_opt.')*(hAU.');
  94.             hE = conj(hAE)*conj(w_opt)*(w_opt.')*(hAE.');
  95.             GU = [diag(conj(hIU))*conj(HAI)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU)  diag(conj(hIU))*conj(HAI)*conj(w_opt)*(w_opt.')*(hAU.');...
  96.                 conj(hAU)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU) 0];
  97.             GE = [diag(conj(hIE))*conj(HAI)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE)  diag(conj(hIE))*conj(HAI)*conj(w_opt)*(w_opt.')*(hAE.');...
  98.                 conj(hAE)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE) 0]; % 公式(18)
  99.             [q,upper_bound] = SDR(hU,hE,GU,GE,N);
  100.             
  101.             Q = diag(q);
  102.             R_old = R_new;
  103.             R = max(0, log2(1+abs((hIU*Q*HAI+hAU)*w_opt)^2)- log2(1+abs((hIE*Q*HAI+hAE)*w_opt)^2));
  104.             R_new = R;
  105.         end
  106.         RSDR(p) = RSDR(p) + R;
  107.         RUB(p) = RUB(p) + log2(upper_bound);
  109.         % AP MRT with IRS
  110.         R_old = 0;
  111.         R_new = 10;
  112.         count = 0;
  113.         while(abs(R_old-R_new)/R_new > epsilon && count < maxIter)
  114.             count = count + 1;
  115.             % 给定q优化W
  116.             w_opt = sqrt(P_AP)*HAI(1,:)'/norm(HAI(1,:));
  117.             
  118.             % 给定W优化q, SDR
  119.             hU = conj(hAU)*conj(w_opt)*(w_opt.')*(hAU.');
  120.             hE = conj(hAE)*conj(w_opt)*(w_opt.')*(hAE.');
  121.             GU = [diag(conj(hIU))*conj(HAI)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU)  diag(conj(hIU))*conj(HAI)*conj(w_opt)*(w_opt.')*(hAU.');...
  122.                 conj(hAU)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU) 0];
  123.             GE = [diag(conj(hIE))*conj(HAI)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE)  diag(conj(hIE))*conj(HAI)*conj(w_opt)*(w_opt.')*(hAE.');...
  124.                 conj(hAE)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE) 0]; % 公式(18)
  125.             [q,~] = SDR(hU,hE,GU,GE,N);
  126.             
  127.             Q = diag(q);
  128.             R_old = R_new;
  129.             R = max(0, log2(1+abs((hIU*Q*HAI+hAU)*w_opt)^2)- log2(1+abs((hIE*Q*HAI+hAE)*w_opt)^2));
  130.             R_new = R;
  131.         end
  132.         RMRT(p) = RMRT(p) + R;
  134.         % without IRS
  135.         A = (hAU)'*(hAU); % 公式(9)
  136.         B = (hAE)'*(hAE); % 公式(10)
  137.         
  138.         I_M = eye(M);
  139.         C = (B+1/P_AP*I_M)\(A+1/P_AP*I_M);
  140.         [V,D] = eig(C); % 特征值分解
  141.         [d,ind] = sort(diag(D));
  142.         u_max = V(:,ind(end))/norm(V(:,ind(end)));
  143.         w_opt = sqrt(P_AP) * u_max;
  144.         R = max(0, log2(1+abs((hAU)*w_opt)^2)- log2(1+abs((hAE)*w_opt)^2));
  145.         RWIRS(p) =RWIRS(p) +R;
  147.         
  148.     end
  149. end
  151. plot(pow2db(P), RSDR/frame,'b-o','LineWidth',2)
  152. hold on
  153. plot(pow2db(P), RMRT/frame,'k-o','LineWidth',2)
  154. plot(pow2db(P), RWIRS/frame,'r-.d','LineWidth',2)
  155. plot(pow2db(P), RUB/frame,'m-.+','LineWidth',2)
  156. grid on
  157. xlabel('P_{AP} (dBm)')
  158. ylabel('Average Secrecy Rate (bps/Hz)')
  159. legend('Alternating Optimization with IRS','AP MRT with IRS','Optimal AP without IRS','Upper bound','Location','northwest')
复制代码
 仿真结果

      
       复现的结果图            
       论文原图       可以看出不同算法的趋势根本复现，数值上大概有些不同，大概还是信道建模部分以及反射面个数的问题，不影响对于算法整体的明白。
看到评论区和私信许多人问关于随着RIS单位数N变革的图，本身改写的步伐始终出问题，因为没有具体调试的代码，不清晰具体错在那里了，我本身改写了重新跑了1000frame也没有出现错误，图像也比较平滑，下面是我跑出来的结果：
   
     复现的结果图   
 

      
       论文原图       这里可以看出随N变革的趋势是一致的，不太一致的地方是我跑的图随着N变革MRT和AO的方法性能会接近，而原论文的图性能差距会变大，这里不是太清晰是不是由于参数设置的问题，或者是没有采用空间相关性信道。
代码链接为： 随N变革部分代码链接
参考文献

[1] A. Khisti and G. W. Wornell, “Secure transmission with multiple antennas I: the MISOME wiretap channel,”IEEE Trans. Inf. Theory, vol. 56, no. 7, pp. 3088-3104, Jul. 2010.

有任何不清晰的写错或步伐有误的地方，欢迎在评论区留言或私信互换！


免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。



欢迎光临 ToB企服应用市场:ToB评测及商务社交产业平台 (https://dis.qidao123.com/) Powered by Discuz! X3.4