引言
Zhang Rui老師的將IRS引入無線通信安全的論文《Secure Wireless Communication via Intelligent Reflecting Surface》有較高的引用量，在此給出要論文的復(fù)現(xiàn)及代碼。
主要問題
該論文的目的是引入IRS并聯(lián)合優(yōu)化基站的主動(dòng)式波束和IRS的被動(dòng)式波束，使得抑制竊聽者信噪比的同時(shí)最大化合法用戶處的信噪比。其場景如下：
?因此可以構(gòu)造出以下的優(yōu)化問題：
?即在基站發(fā)射功率的約束下，優(yōu)化基站和IRS的波束使得保密速率最大化。
給定IRS相位矩陣時(shí)，優(yōu)化基站波束
可簡單地將求絕對(duì)值的平方進(jìn)行簡單展開，令
將對(duì)數(shù)相減變換為真數(shù)相除，對(duì)數(shù)是單調(diào)遞增函數(shù)，因此最大化對(duì)數(shù)，即是最大化真數(shù)即可。因此，可簡化為以下的問題：?
該結(jié)構(gòu)可以參考文獻(xiàn)[1]直接給出解的形式如下：?
?其中對(duì)應(yīng)于矩陣的最大特征值對(duì)于的歸一化特征向量。
給定基站波束，優(yōu)化IRS相位矩陣
該部分推導(dǎo)較為復(fù)雜，可以詳細(xì)閱讀論文，如果有不懂的地方，可以評(píng)論或私信交流。主要是利用了分式規(guī)劃將其轉(zhuǎn)化為一個(gè)半正定松弛問題，求解該問題然后利用高斯隨機(jī)化過程進(jìn)行求解，轉(zhuǎn)化后的問題如下所示：
?仿真復(fù)現(xiàn)
仿真參數(shù)設(shè)置
clc;clear;
epsilon = 1e-3; % 收斂停止條件
% 天線數(shù)
M = 4; % AP天線數(shù)
Nx = 8;
Nz = 8;
N = Nx*Nz; % IRS單元個(gè)數(shù)
% 用戶位置
APloc = [0;0]; % AP位置
userloc = [150;0]; % user位置
edloc= [145;0]; % 竊聽者位置
IRSloc = [145;5]; % IRS位置
C0 = db2pow(-30); % 參考距離時(shí)的路損
D0 = 1; % 參考距離
sigmaK2 = db2pow(-80); % 噪聲功率，-80dBm
L = @(d, alpha)C0*(d/D0)^(-alpha); % 路損模型
% 路損參數(shù)
alpha_AU = 3;
alpha_AE = 3;
alpha_AI = 2.2;
alpha_IU = 3;
alpha_IE = 3;
% 萊斯因子
K_AU = db2pow(1);
K_AE = db2pow(1);
K_AI = db2pow(1);
K_IU = db2pow(1);
K_IE = db2pow(1);
R = 1000; % 信道實(shí)現(xiàn)數(shù)
P_AP = db2pow(25); % 發(fā)射功率15dBm?產(chǎn)生信道
論文中說明信道都為獨(dú)立的萊斯信道，論文中有些信道考慮的是具有空間相關(guān)性的萊斯信道，需要在NLoS部分前后乘以一個(gè)相關(guān)系數(shù)矩陣，具體內(nèi)容可以參考論文[1]，為簡化，在此沒有考慮相關(guān)系數(shù)矩陣，則可以產(chǎn)生如下信道：
dAE = norm(APloc-edloc);
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));
dAU = norm(APloc-userloc);
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));
dAI = norm(APloc-IRSloc);
thetaIRS = atan(145/5);phi = 0; thetaAP = atan(5/145);
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));
dIU = norm(IRSloc-userloc);
thetaIRS = -pi/4; phi = 0;
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));
dIE = norm(IRSloc-edloc);
thetaIRS = 0; phi = 0;
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));?迭代優(yōu)化
q = 2*pi*rand(1,N); % 隨機(jī)初始化IRS的相位
Q = diag(exp(1i*q));
% 給定q優(yōu)化W
A = (hIU*Q*HAI+hAU)'*(hIU*Q*HAI+hAU); % 公式(9)
B = (hIE*Q*HAI+hAE)'*(hIE*Q*HAI+hAE); % 公式(10)
I_M = eye(M);
C = (B+1/P_AP*I_M)\(A+1/P_AP*I_M);
[V,D] = eig(C); % 特征值分解
[d,ind] = sort(diag(D));
u_max = V(:,ind(end))/norm(V(:,ind(end)));
w_opt = sqrt(P_AP) * u_max;
% 給定W優(yōu)化q
hU = conj(hAU)*conj(w_opt)*(w_opt.')*(hAU.');
hE = conj(hAE)*conj(w_opt)*(w_opt.')*(hAE.');
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.');...
conj(hAU)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU) 0];
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.');...
conj(hAE)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE) 0]; % 公式(18)
q = SDR(hU,hE,GU,GE,N);
Q = diag(q);
R = max(0, log2(1+abs((hIU*Q*HAI+hAU)*w_opt)^2)- log2(1+abs((hIE*Q*HAI+hAE)*w_opt)^2))?URA導(dǎo)向矢量函數(shù)
function ura_sv = URA_sv(theta, phi,Nx,Ny)
m = 0:Nx-1;
a_az = exp(1i*pi*m*sin(theta)*cos(phi)).';
n = 0:Ny-1;
a_el = exp(1i*pi*n*sin(phi)).';
ura_sv = kron(a_az,a_el);
endULA導(dǎo)向矢量函數(shù)
function ula_sv = ULA_sv(theta, M)
m = 0:M-1;
ula_sv = exp(1i*pi*m*sin(theta)).';
end半正定松弛優(yōu)化函數(shù)
SDR求解問題(22a)
function [q,cvx_optval] = SDR(hU,hE,GU,GE,N)
L = 1000; % 高斯隨隨機(jī)化次數(shù)
cvx_begin sdp quiet
variable X(N+1,N+1) hermitian
variable mu1(1,1)
maximize(real(trace(GU*X)+mu1*(hU+1)))
subject to
real(trace(GE*X))+mu1*(hE+1)==1;
for i=1:N+1
En = zeros(N+1,N+1);
En(i,i)=1;
real(trace(En*X)) == mu1;
end
X == hermitian_semidefinite(N+1);
mu1 >= 0;
cvx_end
% 高斯隨機(jī)化過程
%% method 1
max_F = 0;
max_q = 0;
S = X / mu1;
[U, Sigma] = eig(S);
for l = 1 : L
r = sqrt(2) / 2 * (randn(N+1, 1) + 1j * randn(N+1, 1));
q = U * Sigma^(0.5) * r;
if q' * GU * q > max_F
max_q = q;
max_F = q' * GU * q;
end
end
q = exp(1j * angle(max_q / max_q(end)));
q = q(1 : N);
end以上程序是給定發(fā)射功率的單點(diǎn)優(yōu)化程序，仿真隨著發(fā)射功率變化的完整程序以及對(duì)比算法如下：
clc;clear;
epsilon = 1e-3; % 收斂停止條件
% 天線數(shù)
M = 4; % AP天線數(shù)
Nx = 8;
Nz = 8;
N = Nx*Nz; % IRS單元個(gè)數(shù)
% 用戶位置
APloc = [0;0]; % AP位置
userloc = [150;0]; % user位置
edloc= [145;0]; % 竊聽者位置
IRSloc = [145;5]; % IRS位置
C0 = db2pow(-30); % 參考距離時(shí)的路損
D0 = 1; % 參考距離
sigmaK2 = db2pow(-80); % 噪聲功率，-80dBm
L = @(d, alpha)C0*(d/D0)^(-alpha); % 路損模型
% 路損參數(shù)
alpha_AU = 3;
alpha_AE = 3;
alpha_AI = 2.2;
% alpha_IU = 3;
alpha_IU = 2.3;
% alpha_IE = 3;
alpha_IE = 2.5;
% 萊斯因子
K_AU = db2pow(1);
K_AE = db2pow(1);
% K_AI = db2pow(1);
K_AI = db2pow(10);
% K_IU = db2pow(1);
K_IU = db2pow(10);
% K_IE = db2pow(1);
K_IE = db2pow(10);
P = db2pow(-5:5:25); % 發(fā)射功率15dBm
frame = 10;
maxIter = 20;
RSDR = zeros(length(P),1);
RMRT = zeros(length(P),1);
RWIRS = zeros(length(P),1);
RUB = zeros(length(P),1);
for p=1:length(P)
p
P_AP = P(p);
for fr = 1:frame
dAE = norm(APloc-edloc);
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));
dAU = norm(APloc-userloc);
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));
dAI = norm(APloc-IRSloc);
thetaIRS = atan(145/5);phi = 0; thetaAP = atan(5/145);
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));
dIU = norm(IRSloc-userloc);
thetaIRS = -pi/4; phi = 0;
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));
dIE = norm(IRSloc-edloc);
thetaIRS = 0; phi = 0;
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));
q = 2*pi*rand(1,N); % 隨機(jī)初始化IRS的相位
Q = diag(exp(1i*q));
R_old = 0;
R_new = 10;
count = 0;
while(abs(R_old-R_new)/R_new > epsilon && count < maxIter)
count = count + 1;
% 給定q優(yōu)化W
A = (hIU*Q*HAI+hAU)'*(hIU*Q*HAI+hAU); % 公式(9)
B = (hIE*Q*HAI+hAE)'*(hIE*Q*HAI+hAE); % 公式(10)
I_M = eye(M);
C = (B+1/P_AP*I_M)\(A+1/P_AP*I_M);
[V,D] = eig(C); % 特征值分解
[d,ind] = sort(diag(D));
u_max = V(:,ind(end))/norm(V(:,ind(end)));
w_opt = sqrt(P_AP) * u_max;
% 給定W優(yōu)化q, SDR
hU = conj(hAU)*conj(w_opt)*(w_opt.')*(hAU.');
hE = conj(hAE)*conj(w_opt)*(w_opt.')*(hAE.');
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.');...
conj(hAU)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU) 0];
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.');...
conj(hAE)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE) 0]; % 公式(18)
[q,upper_bound] = SDR(hU,hE,GU,GE,N);
Q = diag(q);
R_old = R_new;
R = max(0, log2(1+abs((hIU*Q*HAI+hAU)*w_opt)^2)- log2(1+abs((hIE*Q*HAI+hAE)*w_opt)^2));
R_new = R;
end
RSDR(p) = RSDR(p) + R;
RUB(p) = RUB(p) + log2(upper_bound);
% AP MRT with IRS
R_old = 0;
R_new = 10;
count = 0;
while(abs(R_old-R_new)/R_new > epsilon && count < maxIter)
count = count + 1;
% 給定q優(yōu)化W
w_opt = sqrt(P_AP)*HAI(1,:)'/norm(HAI(1,:));
% 給定W優(yōu)化q, SDR
hU = conj(hAU)*conj(w_opt)*(w_opt.')*(hAU.');
hE = conj(hAE)*conj(w_opt)*(w_opt.')*(hAE.');
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.');...
conj(hAU)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIU) 0];
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.');...
conj(hAE)*conj(w_opt)*(w_opt.')*(HAI.')*diag(hIE) 0]; % 公式(18)
[q,~] = SDR(hU,hE,GU,GE,N);
Q = diag(q);
R_old = R_new;
R = max(0, log2(1+abs((hIU*Q*HAI+hAU)*w_opt)^2)- log2(1+abs((hIE*Q*HAI+hAE)*w_opt)^2));
R_new = R;
end
RMRT(p) = RMRT(p) + R;
% without IRS
A = (hAU)'*(hAU); % 公式(9)
B = (hAE)'*(hAE); % 公式(10)
I_M = eye(M);
C = (B+1/P_AP*I_M)\(A+1/P_AP*I_M);
[V,D] = eig(C); % 特征值分解
[d,ind] = sort(diag(D));
u_max = V(:,ind(end))/norm(V(:,ind(end)));
w_opt = sqrt(P_AP) * u_max;
R = max(0, log2(1+abs((hAU)*w_opt)^2)- log2(1+abs((hAE)*w_opt)^2));
RWIRS(p) =RWIRS(p) +R;
end
end
plot(pow2db(P), RSDR/frame,'b-o','LineWidth',2)
hold on
plot(pow2db(P), RMRT/frame,'k-o','LineWidth',2)
plot(pow2db(P), RWIRS/frame,'r-.d','LineWidth',2)
plot(pow2db(P), RUB/frame,'m-.+','LineWidth',2)
grid on
xlabel('P_{AP} (dBm)')
ylabel('Average Secrecy Rate (bps/Hz)')
legend('Alternating Optimization with IRS','AP MRT with IRS','Optimal AP without IRS','Upper bound','Location','northwest')
?仿真結(jié)果
可以看出不同算法的趨勢(shì)基本復(fù)現(xiàn)，數(shù)值上可能有些不同，可能還是信道建模部分以及反射面?zhèn)€數(shù)的問題，不影響對(duì)于算法整體的理解。
看到評(píng)論區(qū)和私信很多人問關(guān)于隨著RIS單元數(shù)N變化的圖，自己改寫的程序始終出問題，因?yàn)闆]有具體調(diào)試的代碼，不清楚具體錯(cuò)在哪里了，我自己改寫了重新跑了1000frame也沒有出現(xiàn)錯(cuò)誤，圖像也比較平滑，下面是我跑出來的結(jié)果：
?
這里可以看出隨N變化的趨勢(shì)是一致的，不太一致的地方是我跑的圖隨著N變化MRT和AO的方法性能會(huì)接近，而原論文的圖性能差距會(huì)變大，這里不是太清楚是不是由于參數(shù)設(shè)置的問題，或者是沒有采用空間相關(guān)性信道。
代碼鏈接為：?隨N變化部分代碼鏈接
參考文獻(xià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.文章來源：http://www.zghlxwxcb.cn/news/detail-414060.html
有任何不清楚的寫錯(cuò)或程序有誤的地方，歡迎在評(píng)論區(qū)留言或私信交流！文章來源地址http://www.zghlxwxcb.cn/news/detail-414060.html
到了這里，關(guān)于智能反射面（IRS）在無線通信安全領(lǐng)域應(yīng)用的論文復(fù)現(xiàn)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
