In mostly channel realizations, the following optimization problem can obtain one solution.

But it appears “infeasible” at the beginning in some channel realizations,

Complete program:

Initialized point

```
Nt = 4 , K = 4 ;
error_radius = 1*ones(K,K);
noise_power = 0.1*ones(K,1);
power_max = 10^(1/10);
epsilon = 0.01 ;
for seed_index =1:50
randn('state',seed_index);
h_bar = 1 / sqrt(2) * (randn(Nt,K,K) + i * randn(Nt,K,K));
wb = zeros(Nt,K) ;
for k=1:K
wb(:,k) = (randn(Nt,1) + i * randn(Nt,1)) ;
norm_wb(1,k) = norm(wb(:,k)) ;
wb(:,k) = sqrt(power_max)*wb(:,k)/norm(1,k) ; % normal distribution = power_max
end
outloop = 0 ;
while(outloop == 0) % *Until Stopping criterion holds*
pre_sum_rate = 0 ;
if (iteration == 0)
s = zeros(K,K) ;
for i = 1:K
h = 0 ;
s(i,i) = max(abs(h_bar(:,i,i)'*w_b(:,i)) - error_radius(i,i)*norm(w_b(:,i)),0.01) ;
for j = 1:K
if (j ~= i)
s(j,i) = abs(h_bar(:,j,i)'*w_b(:,j)) + error_radius(j,i)*norm(w_b(:,j)) ;
h = h + s(j,i)^2 ;
end
end
u_old(1,i) = h + noise_power(i,1) ;
t_old(1,i) = 1 + (s(i,i)^2)/(h + noise_power(i,1)) + epsilon ;
pre_sum_rate = pre_sum_rate + weight(i,1)*log2(t_old(1,i)) ;
end
else
pre_sum_rate = post_sum_rate ;
end
%% CVX for updating
cvx_begin quiet
variables t(1,K) u(1,K) z(1,K-1) lambda(K,K) lambda_hat(1,K) beta1(K,K)
variable W(Nt,Nt,K) complex
variable wb(Nt,K) complex
expression f_ni(1,K)
expression Phi(Nt+1,Nt+1,K,K)
for i = 1:K
f_ni(1,i) = sqrt((t_old(1,i)-1)*u_old(1,i)) + sqrt((t_old(1,i)-1)/u_old(1,i))*(u(1,i)-u_old(1,i))/2 + sqrt(u_old(1,i)/(t_old(1,i)-1))*(t(1,i)-t_old(1,i))/2 ;
for j=1:K
if j~= i
Phi(:,:,j,i) = [eye(Nt) ; h_bar(:,j,i)']*(-W(:,:,j))*([eye(Nt) ; h_bar(:,j,i)']')+ [lambda(j,1)*eye(Nt),zeros(Nt,1) ;zeros(1,Nt),beta1(j,i)-lambda(j,i)*(error_radius(j,i)^2)] ;
end
end
end
maximize (z(1,3))
subject to
norm([2*z(1,1),t(1,1)-t(1,2)]') <= t(1,1)+t(1,2) ; % SOC expressions
norm([2*z(1,2),t(1,3)-t(1,4)]') <= t(1,3)+t(1,4) ;
norm([2*z(1,3),z(1,1)-z(1,2)]') <= z(1,1)+z(1,2) ;
t(1,:) >= ones(1,K) + epsilon ; % t >= (1 + SINR)^1
z(1,1) >= 0 ;
z(1,2) >= 0 ;
for i = 1:K
real(h_bar(:,i,i)'*wb(:,i)) - lambda_hat(1,i)*delta_(i,i) >= f_ni(1,i) ;
norm(wb(:,i)) <= lambda_hat(1,i) ; % SOC
noise_power(i,1) + sum(beta1([1:i-1 i+1:K],i)) <= u(1,i) ; % slack variable
for j=1:K
if j ~= i
lambda(j,i) >= 0 ;
Phi(:,:,j,i) == hermitian_semidefinite(Nt+1) ; % S-lemma
[1,wb(:,j)';wb(:,j),W(:,:,j)] == hermitian_semidefinite(Nt+1) ; % Schur complement
end
end
square_pos(norm(wb(:,i))) <= power_max ; % BS power constraint
end
cvx_end
post_sum_rate = sum(log2(t)) ;
w_b = wb ;
t_old = t ;
u_old = u ;
%% Stopping criterion
diff = post_sum_rate - pre_sum_rate ;
if ( diff > epsilon)
outloop = 0 ;
else
outloop = 1 ;
end
end
end
```