Why is there an error in reproducing this alternate optimization implemented using the SCA algorithm and penalty function？

hugema750 · February 27, 2024, 5:52am

I would like to implement his Algorithm 1 and Algorithm 2 according to the method in the paper. i think i followed his steps strictly in the simulation, but why does the second algorithm appear to be FAILED in the simulation? mark, can you take a look at my code for me, thank you very much!

passage:
W. Lv, J. Bai, Q. Yan and H. M. Wang, “RIS-Assisted Green Secure Communications: Active RIS or Passive RIS?,” in IEEE Wireless Communications Letters , vol. 12, no. 2, pp. 237-241, Feb. 2023, doi: 10.1109/LWC.2022.3221609. keywords: {Artificial intelligence;Fading channels;Wireless communication;Reflection;Optimization;Minimization;Base stations;Reconfigurable intelligent surface;active RIS;double fading;secrecy rate;transmit power},

hugema750 · February 27, 2024, 5:56am

clear all;
clc;
N = 36; % IRS 反射元件数量
M = 12; % 基站天线数量
P = 10; % 基站发射功率
R = 10; %初始化保密率
epsilon=1e-3; %convergence threshold
Lmax=200; %maximum iteration number
yta = 1e15; %罚函数的系数
%% Channel
C0 = db2pow(-30); % 参考距离时的路损
D0 = 1;
sigma = db2pow(-90); % 噪声功率
PL = @(d, alpha)C0*(d/D0)^(-alpha);

Alice=[0,0]; %BS location
IRS=[60,30]; %IRS location
Bob=[70,20]; %Bob location
Eve=[100,20]; %Eve location

ahpahAI = 2.2;
ahpahIB = 3;
ahpahIE = 3;
ahpahAB = 3;
ahpahAE = 3;
kr=db2pow(1); %rician factor

HAI=sqrt(PL(norm(Alice-IRS),ahpahAI))(ones(N,M)+sqrt(1/(kr+1))(sqrt(1/2)*randn(N,M)+sqrt(-1/2)randn(N,M)));
hIB=sqrt(PL(norm(IRS-Bob),ahpahIB))(sqrt(kr/(kr+1))ones(N,1)+sqrt(1/(kr+1))(sqrt(1/2)*randn(N,1)+sqrt(-1/2)randn(N,1)));
hIE=sqrt(PL(norm(IRS-Eve),ahpahIE))(sqrt(kr/(kr+1))ones(N,1)+sqrt(1/(kr+1))(sqrt(1/2)*randn(N,1)+sqrt(-1/2)randn(N,1)));
hAB=sqrt(PL(norm(Alice-Bob),ahpahAB))(sqrt(1/2)*randn(M,1)+sqrt(-1/2)randn(M,1));
hAE=sqrt(PL(norm(Alice-Eve),ahpahAE))(sqrt(1/2)*randn(M,1)+sqrt(-1/2)*randn(M,1));

%%
ite=0; %initialize iteration number
obj=zeros(1,Lmax); %store objective value in each iteration

%================================固定Q优化W================================
q=ones(N,1); %initialize RIS coefficients
Q=diag(q’);
hB = HAI’QhIB+hAB;
hE = HAI’QhIE+hAE;
A = 1/(sigma*(1+norm(hIB’Q,2)^2));
B = 1/(sigma(1+norm(hIE’*Q,2)^2));
Wt=zeros(M,M);
[V, D] = eig(Wt); % 计算矩阵A的特征值和特征向量
% 找到最大特征值及其对应的特征向量
[lamda, idx] = max(diag(D));
v = V(:, idx)/norm(V(:, idx));
z=exp(R)-1;
f=P-(norm(Q,2)^2)*sigma;
while(1)

cvx_begin sdp
    variable W(M,M) hermitian semidefinite
    minimize (real(trace(W)+yta*(trace(W)-lamda-trace(v*v'*(W-Wt)))))
    subject to
        real(A*trace(hB*hB'*W)-B*trace(hE*hE'*W))>=z;
        real(trace(Q*HAI*W*HAI'*Q'))<=f; 
cvx_end

aresult=real(A*trace(hB*hB'*W)-B*trace(hE*hE'*W))-z
aaresult=real(trace(Q*HAI*W*HAI'*Q'))-f
[V, D] = eig(W); % 计算矩阵A的特征值和特征向量
[lamda, idx] = max(diag(D)); % 找到最大特征值及其对应的特征向量
v = V(:, idx)/norm(V(:, idx));
Wt=W;
if(abs(real(trace(Wt)-lamda))<epsilon)
    atest=abs(real(trace(Wt)-lamda))-epsilon
    break;
end
atw=real(trace(Wt));

end
r=rank(Wt)
w=sqrt(lamda)v;%找到了离最优的Wt最近的w
aaaaaaaaaaaaaaa=ww’-Wt
%================================固定W优化Q================================
GB=diag(hIB’)HAI;
GE=diag(hIE’)HAI;
G=[diag(HAIw)diag(HAIw)'+sigmaeye(N) zeros(N,1);zeros(1,N) 0];
GIB=[sigmadiag(hIB’)diag(hIB’)’ zeros(N,1);zeros(1,N) sigma];
GIE=[sigmadiag(hIE’)diag(hIE’)’ zeros(N,1);zeros(1,N) sigma];
GAB=[GBww’GB’+sigmadiag(hIB’)diag(hIB’)’ GBww’hAB;hAB’ww’GB’ sigma+hAB’ww’hAB];
GAE=[GEww’GE’+sigmadiag(hIE’)diag(hIE’)’ GEww’hAE;hAE’ww’GE’ sigma+hAE’ww’hAE];
%进行一些初始化变量的设置
%设置迭代起始条件（表示了一个r=1的矩阵）
Ut=[1ones(N,1);1][1ones(N,1);1]';
[V, D] = eig(Ut); % 计算矩阵A的特征值和特征向量
% 找到最大特征值及其对应的特征向量
[lamda, idx] = max(diag(D));
v = V(:, idx)/norm(V(:, idx));
%wcount记录迭代次数
wcount=0;
%记录一些特定的变量
CIB=GIB/real(trace(GIBUt));
CAE=GAE/real(trace(GAE*Ut));
while(1)

cvx_begin sdp
    variable U(N+1,N+1) hermitian semidefinite
    variable deta(1,1);
    maximize (real(deta-yta*(trace(U)-lamda-trace(v*v'*(U-Ut)))));
    subject to
        deta>=0;
        %（1）约束
        real(-rel_entr(1,real(trace(GAB*U)))-rel_entr(1,real(trace(GIE*U)))+rel_entr(1,real(trace(GIB*Ut)))+rel_entr(1,real(trace(GAE*Ut))) ...
            -real(trace(CIB*(U-Ut)))-real(trace(CAE*(U-Ut)))-deta)>=R;
        real(trace(G*U))<=P;
        real(U(N+1,N+1))==1;  
cvx_end

% result=  real(log(real(trace(GAB*U)))+log(real(trace(GIE*U)))-log(real(trace(GIB*Ut)))-log(real(trace(GAE*Ut))) ...
%             -real(trace(CIB*(U-Ut)))-real(trace(CAE*(U-Ut)))-deta);
result=real(-rel_entr(1,real(trace(GAB*U)))-rel_entr(1,real(trace(GIE*U)))+rel_entr(1,real(trace(GIB*Ut)))+rel_entr(1,real(trace(GAE*Ut))) ...
            -real(trace(CIB*(U-Ut)))-real(trace(CAE*(U-Ut)))-deta)
[V, D] = eig(U); % 计算矩阵A的特征值和特征向量
[lamda, idx] = max(diag(D)); % 找到最大特征值及其对应的特征向量
v = V(:, idx)/norm(V(:, idx));  

if(abs(real(trace(G*U)-trace(G*Ut)))<=epsilon)
 break;
end

%更新数值继续迭代
Ut=U;
CIB=GIB/real(trace(GIB*Ut));
CAE=GAE/real(trace(GAE*Ut));

end
%？？？？？
%？？？？？？疑问：yta设置小的话，固定W优化Q可能不会满足r=1，为什么最后可以得解但是result记录的值小于R，这并不满足上方的（1）约束啊
%？？？？？？但是如果yta设置太大，求解固定w优化Q会出现failed
%？？？？？？
%？？？？？

%经过以上运算得到的Ut是最佳的，进而推出最优的Q
v1=v*sqrt(lamda);
Q=diag(v1(1:N));

Mark_L_Stone · February 27, 2024, 12:06pm

Getting SCA to run nicely is out of scope of this forum. SCA is unreliable, and many of the algorithms are too complicated for volunteers with limited time to be able to spend adequate effort to figure out things.

hugema750 · February 28, 2024, 2:58pm

Thank you very much for your answer and have a nice life!

alemayehu_kebede · March 23, 2024, 2:43pm

hello ,@hugema 750 if you correct cvx opt. IRS-phase error would you message me via this platform. thank you