In the first two iterations, there is no error. However, there is a error in the third iteration.
i have prove this problem is convex, but i have know why the error occurs. Can you help me?
clc,clear;
close all;
%%
K=2;
Nt = 4;
M = 10;
iter_num = 100;
inner_iter_num = 1000;
h_irs_u = rand(M,K)+1i*rand(M,K);
h_ai = rand(M,Nt)+1i*rand(M,Nt);
W = zeros(Nt,K,iter_num+1);
val = zeros(1,iter_num+1);
mu = zeros(K,iter_num+1);
v = zeros(K,iter_num+1);
Phi1= zeros(M,M,iter_num+1);
Phi1(:,:,1) = sqrt(0.5)*diag(exp(1i*2*pi*rand(M,1)));
Phi2= zeros(M,M,iter_num+1);
Phi2(:,:,1) = sqrt(0.5)*diag(exp(1i*2*pi*rand(M,1)));
Pmax = 3;
sigma = 1;
delta = 1e-3;
% initialization
W(:,:,1) = rand(Nt,K)+1i*rand(Nt,K);
mu(:,1) = rand(K,1);
v(:,1) = rand(K,1);
rho = 1.5;
for iter = 1:iter_num
% optimize mu and v
% mu_k = zeros(K,1);
% v_k = zeros(K,1);
H = [h_irs_u(:,1)'*Phi1(:,:,iter)*h_ai;
h_irs_u(:,2)'*Phi2(:,:,iter)*h_ai]';
for kk=1:K
rho = real(conj(v(kk,iter))*H(:,kk)'*W(:,kk,iter));
mu(kk,iter+1)= rho/2*(rho+sqrt(rho^2+4));
v(kk,iter+1) = sqrt(1+mu(kk,iter+1))*H(:,kk)'*W(:,kk,iter)/(norm(H(:,kk)'*W(:,:,iter))^2+sigma);
end
% optimize W
cvx_begin quiet
variable W0(Nt,K) complex
maximize (2*sqrt(1+mu(1,iter+1))*real(v(1,iter+1)'*H(:,1)'*W0(:,1))-abs(v(1,iter+1))^2*(square_pos(norm(H(:,1)'*W0))+sigma)...
+2*sqrt(1+mu(2,iter+1))*real(v(2,iter+1)'*H(:,2)'*W0(:,2))-abs(v(2,iter+1))^2*(square_pos(norm(H(:,2)'*W0))+sigma))
subject to
norm(W0)<=Pmax;
cvx_end
W(:,:,iter+1)=W0;
% optimize phase shift at IRS
Theta1 = zeros(M,inner_iter_num+1);
Theta1(:,1) = diag(Phi1(:,:,iter));
Theta2 = zeros(M,inner_iter_num+1);
Theta2(:,1) = diag(Phi2(:,:,iter));
psi1=zeros(M,inner_iter_num+1);
psi1(:,1) = exp(1i*2*pi*rand(M,1));
psi2= zeros(M,inner_iter_num+1);
psi2(:,1) = exp(1i*2*pi*rand(M,1));
u1 = zeros(M,inner_iter_num+1);
u1(:,1) = rand(M,1);
u2 = zeros(M,inner_iter_num+1);
u2(:,1) = rand(M,1);
x = zeros(1,inner_iter_num+1);
x(1) = sqrt(0.5);
y = zeros(1,inner_iter_num+1);
y(1)=sqrt(0.5);
for inner_iter=1:inner_iter_num
% optimize theta1 and theta2
cvx_begin quiet
variable theta1(M,1) complex
variable theta2(M,1) complex
H1 = [h_irs_u(:,1)'*diag(theta1)*h_ai;
h_irs_u(:,2)'*diag(theta2)*h_ai]';
maximize(2*sqrt(1+mu(1,iter+1))*real(v(1,iter+1)'*H1(:,1)'*W0(:,1))-abs(v(1,iter+1))^2*(square_pos(norm(H1(:,1)'*W0))+sigma)...
+2*sqrt(1+mu(2,iter+1))*real(v(2,iter+1)'*H1(:,2)'*W0(:,2))-abs(v(2,iter+1))^2*(square_pos(norm(H1(:,2)'*W0))+sigma)...
-rho/2*(square_pos(norm(theta1-x(inner_iter)*psi1(:,inner_iter)+u1(:,inner_iter)))+square_pos(norm(theta2-y(inner_iter)*psi2(:,inner_iter)+u2(:,inner_iter)))))
cvx_end
Theta1(:,inner_iter+1) = theta1;
Theta2(:,inner_iter+1) = theta2;
% optimze psi1 and psi2
psi1(:,inner_iter+1) = exp(1i*angle(theta1+u1(:,inner_iter)));
psi2(:,inner_iter+1) = exp(1i*angle(theta2+u2(:,inner_iter)));
% optimize x and y
a = 2*real((Theta1(:,inner_iter+1)+u1(:,inner_iter))'*psi1(:,inner_iter+1));
c = 2*real((Theta2(:,inner_iter+1)+u2(:,inner_iter))'*psi2(:,inner_iter+1));
d = -(M+norm(Theta1(:,inner_iter+1)+u1(:,inner_iter))^2+norm(Theta2(:,inner_iter+1)+u2(:,inner_iter))^2);
x(inner_iter+1) = a/sqrt(a^2+c^2);
y(inner_iter+1) = c/sqrt(a^2+c^2);
% optimze u1 and u2
u1(:,inner_iter+1) = u1(:,inner_iter) + Theta1(:,inner_iter+1) -x(inner_iter+1)*psi1(:,inner_iter+1);
u2(:,inner_iter+1) = u2(:,inner_iter) + Theta2(:,inner_iter+1) -y(inner_iter+1)*psi2(:,inner_iter+1);
% disp([norm(Theta1(:,inner_iter+1)-x(inner_iter+1)*psi1(:,inner_iter+1)), norm(Theta2(:,inner_iter+1)-y(inner_iter+1)*psi2(:,inner_iter+1))]);
err1 = norm(Theta1(:,inner_iter+1)-x(inner_iter+1)*psi1(:,inner_iter+1));
err2 = norm(Theta2(:,inner_iter+1)-y(inner_iter+1)*psi2(:,inner_iter+1));
if err1<=delta && err2<=delta
break
end
end
Phi1(:,:,iter+1) = diag(Theta1(:,inner_iter+1));
Phi2(:,:,iter+1) = diag(Theta2(:,inner_iter+1));
H = [h_irs_u(:,1)'*Phi1(:,:,iter+1)*h_ai;
h_irs_u(:,2)'*Phi2(:,:,iter+1)*h_ai]';
val(iter+1) = log2(1+abs(H(:,1)'*W0(:,1))^2/(abs(H(:,1)'*W0(:,2))^2+sigma)) + log2(1+abs(H(:,2)'*W0(:,2))^2/(abs(H(:,2)'*W0(:,1))^2+sigma));
disp([iter,val(iter+1) ])
if abs(val(iter+1) - val(iter))<=delta
break
end
end
figure
plot(val(2:iter+1),'b-','linewidth',1.5);
xlabel('iteration index');
ylabel('Objective value');
grid on