When i use cvx to solve this convex problem, something wrong, can you help me?

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``````

Check the numerical values of the coefficients used in the objective function on the iteration in which the error occurs.

You can read many posts by me on this forum of the hazards of crude unsafeguarded SCA. What you have experienced is perhaps one such manifestation.

ok, thank you very much