Below provided code is taking many days (even week) to run completely, even though using CVX Version 2.2, and Mosek solver (Version 9.1.9) is selected. How to handle such big codes, is there setup/solver correction need to be done? (even with other solvers also taking very high time than this)
% This is a code package related to the following conference paper:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Y. Mao, B. Clerckx and V. O. K. Li, “Rate-splitting multiple access for
% downlink communication systems: bridging, generalizing, and outperforming
% SDMA and NOMA.” EURASIP Journal on Wireless Communications and Networking
% 2018.1 (2018): 133.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The code is implemented in Matlab environment with CVX toolbox assisted.
%
% Fig. 6(b) of the above paper will be reproduced by running this Matlab
% script. By changing the variables ‘bias’ (channel gain difference between
% the users), ‘NT’(number of transmit antenna), Fig. 5–Fig. 6 can be
% reproduced.
clc;clear all; clf
tic
%number of transmit antenna
NT=4;
%channel gain difference
bias=0.3;
%SNR in dB
SNRdB=20;
%accuracy of convergence
tolerance = 10^-6;
%tolerance = 10^-2;
%tolerance = 10^-1;
%number of random channel
N_channel=100;
%N_channel=25;
%N_channel=3;
%user weights
weight=[-3 -1:0.05:1 3];
%weight=[-3 -1:1:1 3];
%weight=[-3 0 3];
u2=10.^weight;
u1=ones(1,length(u2));
capacity_DPC_UE1_average =zeros(length(u1),1);
capacity_DPC_UE2_average =zeros(length(u1),1);
capacity_MULP_UE1_average=zeros(length(u1),1);
capacity_MULP_UE2_average=zeros(length(u1),1);
capacity_NOMA_order1_UE1_average=zeros(length(u1),1);
capacity_NOMA_order1_UE2_average=zeros(length(u1),1);
capacity_NOMA_order2_UE1_average=zeros(length(u1),1);
capacity_NOMA_order2_UE2_average=zeros(length(u1),1);
capacity_RS_order1_UE1_average=zeros(length(u1),1);
capacity_RS_order1_UE2_average=zeros(length(u1),1);
capacity_RS_order2_UE1_average=zeros(length(u1),1);
capacity_RS_order2_UE2_average=zeros(length(u1),1);
for i_weight=1:length(u1)
for i_channel=1:N_channel
weights=[u1(i_weight),u2(i_weight)];
[i_weight, i_channel]
randn('seed',i_channel*4) %Correct (Authors given)
%randn('state',i_channel*4) %not correct
H_BC(:,:,1)=1/sqrt(2)*(randn(1,NT)+1i*randn(1,NT));
%H_BC(:,:,1)=1/sqrt(2)*(randn(1,NT)+1i*randn(1,NT))*bias; %not correct
H_BC(:,:,2)=1/sqrt(2)*(randn(1,NT)+1i*randn(1,NT))*bias;
H_MAC(:,:,1)=H_BC(:,:,1)';
H_MAC(:,:,2)=H_BC(:,:,2)';
Capacity_DPC = DPC_rateRegion(weights,H_MAC,SNRdB,tolerance);
[Capacity_RS_order1, Capacity_RS_order2, P_common1, P_common2] = RS_rateRegion(weights,H_BC,SNRdB,tolerance);
[Capacity_NOMA_order1, Capacity_NOMA_order2 ]= NOMA_rateRegion(weights,H_BC,SNRdB,tolerance);
Capacity_MULP = MULP_rateRegion(weights,H_BC,SNRdB,tolerance);
capacity_DPC_UE1(i_channel)=Capacity_DPC(1);
capacity_DPC_UE2(i_channel)=Capacity_DPC(2);
capacity_MULP_UE1(i_channel)=Capacity_MULP(1);
capacity_MULP_UE2(i_channel)=Capacity_MULP(2);
capacity_NOMA_order1_UE1(i_channel)=Capacity_NOMA_order1(1);
capacity_NOMA_order1_UE2(i_channel)=Capacity_NOMA_order1(2);
capacity_NOMA_order2_UE1(i_channel)=Capacity_NOMA_order2(1);
capacity_NOMA_order2_UE2(i_channel)=Capacity_NOMA_order2(2);
capacity_RS_order1_UE1(i_channel) = Capacity_RS_order1(1);
capacity_RS_order1_UE2(i_channel) = Capacity_RS_order1(2);
capacity_RS_order2_UE1(i_channel) = Capacity_RS_order2(1);
capacity_RS_order2_UE2(i_channel) = Capacity_RS_order2(2);
end %end i_channel
capacity_DPC_UE1_average(i_weight)=mean(capacity_DPC_UE1);
capacity_DPC_UE2_average(i_weight)=mean(capacity_DPC_UE2);
capacity_MULP_UE1_average(i_weight)=mean(capacity_MULP_UE1);
capacity_MULP_UE2_average(i_weight)=mean(capacity_MULP_UE2);
capacity_NOMA_order1_UE1_average(i_weight)=mean(capacity_NOMA_order1_UE1);
capacity_NOMA_order1_UE2_average(i_weight)=mean(capacity_NOMA_order1_UE2);
capacity_NOMA_order2_UE1_average(i_weight)=mean(capacity_NOMA_order2_UE1);
capacity_NOMA_order2_UE2_average(i_weight)=mean(capacity_NOMA_order2_UE2);
capacity_RS_order1_UE1_average(i_weight)=mean(capacity_RS_order1_UE1);
capacity_RS_order1_UE2_average(i_weight)=mean(capacity_RS_order1_UE2);
capacity_RS_order2_UE1_average(i_weight)=mean(capacity_RS_order2_UE1);
capacity_RS_order2_UE2_average(i_weight)=mean(capacity_RS_order2_UE2);
clear capacity_DPC_UE1
clear capacity_DPC_UE2
clear capacity_MULP_UE1
clear capacity_MULP_UE2
clear capacity_NOMA_order1_UE1
clear capacity_NOMA_order1_UE2
clear capacity_NOMA_order2_UE1
clear capacity_NOMA_order2_UE2
clear capacity_RS_order1_UE1
clear capacity_RS_order1_UE2
clear capacity_RS_order2_UE1
clear capacity_RS_order2_UE2
end %end for i_weight
toc
figure(1)
plot(reshape(capacity_DPC_UE1_average,[length(u1),1]),reshape(capacity_DPC_UE2_average,[length(u2),1]),‘-o’,‘LineWidth’,2,‘Color’, [ 0 0.4470 0.7410])
hold on,grid on
x=[capacity_RS_order1_UE1_average capacity_RS_order2_UE1_average];
y=[capacity_RS_order1_UE2_average capacity_RS_order2_UE2_average];
k=convhull(x,y);
x1 = x(k);
y1 = y(k);
xx=floor(x1.*10^(5))./(10^(5));
ind=find(xx==0);
[~,ind_ini]=max(x1);
plot(x1(ind_ini(1):ind),y1(ind_ini(1):ind),‘-+’,‘LineWidth’,2,‘Color’, [0.8500 0.3250 0.0980])
hold on,grid on
x=[capacity_NOMA_order1_UE1_average capacity_NOMA_order2_UE1_average];
y=[capacity_NOMA_order1_UE2_average capacity_NOMA_order2_UE2_average];
k=convhull(x,y);
x1 = x(k);
y1 = y(k);
xx=floor(x1.*10^(5))./(10^(5));
ind=find(xx==0);
[~,ind_ini]=max(x1);
ind=ind(find(ind>ind_ini));
plot(x1(ind_ini(1):ind(1)),y1(ind_ini(1):ind(1)),‘–’,‘LineWidth’,2,‘Color’,[0.9290 0.6940 0.1250 ])
hold on,grid on
x=capacity_MULP_UE1_average;
y=capacity_MULP_UE2_average;
k=convhull(x,y);
x1 = x(k);
y1 = y(k);
xx=floor(x1.*10^(5))./(10^(5));
ind=find(xx==0);
[~,ind_ini]=max(x1);
plot(x1(ind_ini(1):ind(1)),y1(ind_ini(1):ind(1)),‘-.’,‘LineWidth’,2,‘Color’, [0.4660 0.6740 0.1880])
legend(‘DPC’,‘RS’,‘SC-SIC’,‘MU-LP’)
xlabel(‘R_1 (bit/s/Hz)’)
ylabel(‘R_2 (bit/s/Hz)’)
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% Functions defining
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [ Capacity1, Capacity2, P_common1, P_common2] = RS_rateRegion(weights,H,SNRdB,tolerance)
[Capacity1, P_common1 ] = RS_rateRegion_order1(weights,H,SNRdB,tolerance);
H1(:,:,1)=H(:,:,2);
H1(:,:,2)=H(:,:,1);
weights1(1)=weights(2);
weights1(2)=weights(1);
[Capacity22, P_common2] = RS_rateRegion_order1(weights1,H1,SNRdB,tolerance);
Capacity2(2)=Capacity22(1);
Capacity2(1)=Capacity22(2);
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [ Capacity,P_common] = RS_rateRegion_order1(weights,H,SNRdB,tolerance)
SNR = 10^(SNRdB/10);
[A_U,NT,N_user] = size(H) ;
h1=H(:,:,1);
h2=H(:,:,2);
H_tot=[h1’ h2’];
%Step 1: Initialization–MRC+SVD
P_common=SNR*0.9;
P_private_k=(SNR-P_common)/N_user;
[U2,~,~]=svd(H_tot);
hat_p_c=U2(:,1);
p_1=h1’/norm(h1)*sqrt(P_private_k);
p_2=h2’/norm(h2)sqrt(P_private_k);
p_c=hat_p_csqrt(P_common);
loop=1;
cap_tot_past=0;
count=0;
while (loop)
%T_ck and T_k calculation
T_1=abs(h1*p_1)^2+abs(h1*p_2)^2+1;
T_2=abs(h2*p_1)^2+abs(h2*p_2)^2+1;
T_c_1=abs(h1*p_c)^2+T_1;
T_c_2=abs(h2*p_c)^2+T_2;
I_1=T_1-abs(h1*p_1)^2;
I_2=T_2-abs(h2*p_2)^2;
%Step 2: MMSE Common and private Receiver update at each user
g_c_1=inv(T_c_1)*p_c'*h1';
g_c_2=inv(T_c_2)*p_c'*h2';
g_1=inv(T_1)*p_1'*h1';
g_2=inv(T_2)*p_2'*h2';
%Step 3: MSE Common and private Weight update at each user
E_c_1=inv(T_c_1)*T_1;
E_c_2=inv(T_c_2)*T_2;
E_1=inv(T_1)*I_1;
E_2=inv(T_2)*I_2;
W_c_1=inv(E_c_1);
W_c_2=inv(E_c_2);
W_1=inv(E_1);
W_2=inv(E_2);
%Step 4: Update P
[cap_tot, p_1, p_2, p_c]=RS_update_P_order1(weights,H,SNR,g_c_1,g_c_2,g_1,g_2,W_c_1,W_c_2,W_1,W_2);
if abs(cap_tot-cap_tot_past)<=tolerance
loop=0;
else
cap_tot_past=cap_tot;
count=count+1;
end
if count>=2000
break;
end
end
P_common=real(trace(p_c*p_c’));
%Calculate the capacity of each UE
T_1=abs(h1p_1)^2+abs(h1p_2)^2+1;
T_2=abs(h2p_1)^2+abs(h2p_2)^2+1;
T_c_1=abs(h1p_c)^2+T_1;
T_c_2=abs(h2p_c)^2+T_2;
I_1=T_1-abs(h1p_1)^2;
I_2=T_2-abs(h2p_2)^2;
%Common Rate
R_c_1=real(log2(T_c_1/T_1));
R_c_2=real(log2(T_c_2/T_2));
%Private Rate
R_1=real(log2(T_1/I_1));
R_2=real(log2(T_2/I_2));
%ro_c is not a optimization variable in timesharing scenario
Capacity(1)=R_1+min(R_c_1,R_c_2);
Capacity(2)=R_2;
Capacity_common(1)=R_c_1;
Capacity_common(2)=R_c_1;
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [cap_tot, p_1, p_2, p_c ]=RS_update_P_order1(weights,H,SNR,g_c_1,g_c_2,g_1,g_2,W_c_1,W_c_2,W_1,W_2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Then we try to use CVX method to solve the problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u1=weights(1);
u2=weights(2);
[A_U,NT,N_user] = size(H);
h1=H(:,:,1);
h2=H(:,:,2);
cvx_begin quiet
variable p_1(NT,A_U) complex
variable p_2(NT,A_U) complex
variable p_c(NT,A_U) complex
% variable ro_c(N_user)
expression constraints;
%objective function
%common rate part:
T_c_1=square_abs(h1*p_c)+square_abs(h1*p_1)+square_abs(h1*p_2)+1;
T_c_2=square_abs(h2*p_c)+square_abs(h2*p_1)+square_abs(h2*p_2)+1;
E_c_1=square_abs(g_c_1)*T_c_1-2*real(g_c_1*h1*p_c)+1;
E_c_2=square_abs(g_c_2)*T_c_2-2*real(g_c_2*h2*p_c)+1;
c_1=(W_c_1*E_c_1-log2(W_c_1));
c_2=(W_c_2*E_c_2-log2(W_c_2));
c=max(c_1,c_2);
%we assume common rate is for user 1, therefore, c_2=0
object_common=u1*c;
%Private rate part
T_1=square_abs(h1*p_1)+square_abs(h1*p_2)+1;
T_2=square_abs(h2*p_1)+square_abs(h2*p_2)+1;
E_1=abs(g_1)^2*T_1-2*real(g_1*h1*p_1)+1;
E_2=abs(g_2)^2*T_2-2*real(g_2*h2*p_2)+1;
object_private=u1*(W_1*E_1-log2(W_1))+u2*(W_2*E_2-log2(W_2));
object_func=object_common+object_private;
minimize(object_func)
%constraints
constraints=trace(p_1'*p_1)+trace(p_2'*p_2)+trace(p_c'*p_c)-SNR;
subject to
constraints<=zeros(size(constraints))
cvx_end
cap_tot=object_func;
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [Capacity1 Capacity2 ] = NOMA_rateRegion(weights,H,SNRdB,tolerance)
[Capacity1] = NOMA_rateRegion_order1(weights,H,SNRdB,tolerance);
%switch the order of UE1 and UE2
H1(:,:,1)=H(:,:,2);
H1(:,:,2)=H(:,:,1);
weights1(1)=weights(2);
weights1(2)=weights(1);
[Capacity22] = NOMA_rateRegion_order1(weights1,H1,SNRdB,tolerance);
Capacity2(1)=Capacity22(2);
Capacity2(2)=Capacity22(1);
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [Capacity]=NOMA_rateRegion_order1(weights,H,SNRdB,tolerance)
% Here assumed user-1 is Far, and user-2 is Near user.
SNR = 10^(SNRdB/10);
[A_U,NT,N_user] = size(H) ;
h1=H(:,:,1);
h2=H(:,:,2);
%Step 1: Initialization
P_k=(SNR)/N_user;
p_1=h1’/norm(h1)*sqrt(P_k);
p_2=h2’/norm(h2)*sqrt(P_k);
loop=1;
cap_tot_past=0;
count=0;
while (loop)
%T_ck and T_k calculation
T_1_1=abs(h1*p_1)^2+abs(h1*p_2)^2+1;
T_2_1=abs(h2*p_1)^2+abs(h2*p_2)^2+1; %to decode user-1 signal at user-2
T_2_2=abs(h2*p_2)^2+1;
I_1_1=abs(h1*p_2)^2+1;
I_2_1=abs(h2*p_2)^2+1;
I_2_2=1;
%Step 2: MMSE Common and private Receiver update at each user
g_1_1=inv(T_1_1)*p_1'*h1';
g_2_1=inv(T_2_1)*p_1'*h2';
g_2_2=inv(T_2_2)*p_2'*h2';
%Step 3: MSE Common and private Weight update at each user
E_1_1=inv(T_1_1)*I_1_1;
E_2_1=inv(T_2_1)*I_2_1;
E_2_2=inv(T_2_2)*I_2_2;
W_1_1=inv(E_1_1);
W_2_1=inv(E_2_1);
W_2_2=inv(E_2_2);
%optimization
[cap_tot, p_1, p_2 ]=NOMA_update_P_order1(weights,H,SNR,g_1_1,g_2_1,g_2_2,W_1_1,W_2_1,W_2_2);
if abs(cap_tot-cap_tot_past)<=tolerance
loop=0;
else
cap_tot_past=cap_tot;
count=count+1;
end
if count>=2000
break;
end
end
%Calculate the capacity of each UE
T_2_2=abs(h2p_2)^2+1;
I_2_2=T_2_2-abs(h2p_2)^2;
T_1_1=abs(h1p_1)^2+abs(h1p_2)^2+1;
T_2_1=abs(h2p_1)^2+abs(h2p_2)^2+1;
I_1_1=abs(h1p_2)^2+1;
I_2_1=abs(h2p_2)^2+1;
R_1_1=real(log2(T_1_1/I_1_1));
R_2_1=real(log2(T_2_1/I_2_1));
%Private Rate
R_2=real(log2(T_2_2/I_2_2));
Capacity(1)=min(R_1_1,R_2_1);
Capacity(2)=R_2;
% Capacity
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [cap_tot, p_1, p_2 ]=NOMA_update_P_order1(weights,H,SNR,g_1_1,g_2_1,g_2_2,W_1_1,W_2_1,W_2_2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Then we try to use CVX to solve the problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u1=weights(1);
u2=weights(2);
[A_U,NT,N_user] = size(H);
h1=H(:,:,1);
h2=H(:,:,2);
cvx_begin quiet
variable p_2(NT,A_U) complex
variable p_1(NT,A_U) complex
expression constraints(1,2*N_user);
T_1_1=square_abs(h1*p_1)+square_abs(h1*p_2)+1;
T_2_1=square_abs(h2*p_1)+square_abs(h2*p_2)+1;
E_1_1=square_abs(g_1_1)*T_1_1-2*real(g_1_1*h1*p_1)+1;
E_2_1=square_abs(g_2_1)*T_2_1-2*real(g_2_1*h2*p_1)+1;
r_1_1=(W_1_1*E_1_1-log2(W_1_1));
r_2_1 =(W_2_1*E_2_1-log2(W_2_1));
%objective function
object_common=u1*max(r_1_1,r_2_1);
T_2_2=square_abs(h2*p_2)+1;
E_2_2=abs(g_2_2)^2*T_2_2-2*real(g_2_2*h2*p_2)+1;
object_private=u2*(W_2_2*E_2_2-log2(W_2_2));
object_func=object_common+object_private;
minimize(object_func)
%constraints
constraints=trace(p_2'*p_2)+trace(p_1'*p_1)-SNR;
subject to
constraints<=zeros(size(constraints))
cvx_end
cap_tot=object_func;
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% MULP_rateRegion
function Capacity = MULP_rateRegion(weights,H,SNRdB,tolerance)
%function Capacity = WMMSE_rateRegion(weights,H,SNRdB,tolerance)
%------------------Algorithm WMMSE Algorithm in Paper: [Weighted
%Sum-Rate Maximization using Weighted MMSE for MIMO-BC
%Beamforming Design]-------------------------------------------
SNR = 10^(SNRdB/10);
[A_U,NT,N_user] = size(H);
H_tot=;
for i_user = 1:N_user
H_tot = cat(1,H_tot,H(:,:,i_user)); %size: (A_U*N_user)*NT
end
%Step 1: initialization B
P_k=SNR/N_user;
for i_user=1:N_user
h1=H(:,:,i_user);
B(:,:,i_user)=h1’/norm(h1)*sqrt(P_k);
end
loop=1;
cap_tot_past=100000;
count=0;
while(loop)
%Step 2: MMSE Receiver Calculation
for i_user =1 : N_user
% B
R(:,:,i_user)=eye(A_U);
for j_user=1:N_user
if j_user~=i_user
R(:,:,i_user)=R(:,:,i_user)+H(:,:,i_user)*B(:,:,j_user)*B(:,:,j_user)‘*H(:,:,i_user)’;
end
end
A_MMSE{i_user}=B(:,:,i_user)‘*H(:,:,i_user)’*inv(H(:,:,i_user)B(:,:,i_user)B(:,:,i_user)'H(:,:,i_user)'+R(:,:,i_user));
end
A_MMSE_blkdiag = blkdiag(A_MMSE{:}); %size: (A_UN_user)(A_UN_user)
%Step 3: MSE Weight
for i_user=1:N_user
E(:,:,i_user)=inv(eye(A_U)+B(:,:,i_user)'*H(:,:,i_user)'*inv(R(:,:,i_user))*H(:,:,i_user)*B(:,:,i_user));
W{i_user}=weights(i_user)*inv(E(:,:,i_user));
end
W_blkdiag = blkdiag(W{:}); %size: (A_U*N_user)*(A_U*N_user)
%Step 4: WMMSE transmit filter
B_bar=inv(H_tot'*A_MMSE_blkdiag'*W_blkdiag*A_MMSE_blkdiag*H_tot+trace(W_blkdiag*A_MMSE_blkdiag*A_MMSE_blkdiag')*eye(NT)/SNR)*H_tot'*A_MMSE_blkdiag'*W_blkdiag;
b=sqrt(SNR/trace(B_bar*B_bar'));
B_MMSE=b*B_bar; %size: NT*(A_U*N_user)
for i_user=1:N_user
B(:,:,i_user)=B_MMSE(:,(i_user-1)*A_U+1:i_user*A_U);
end
%loop determination
for i_user=1:N_user
R(:,:,i_user)=eye(A_U);
for j_user=1:N_user
if j_user~=i_user
R(:,:,i_user)=R(:,:,i_user)+H(:,:,i_user)*B(:,:,j_user)*B(:,:,j_user)'*H(:,:,i_user)';
end
end
E(:,:,i_user)=inv(eye(A_U)+B(:,:,i_user)'*H(:,:,i_user)'*inv(R(:,:,i_user))*H(:,:,i_user)*B(:,:,i_user));
Capacity_loop(i_user)=real(log2(det(inv(E(:,:,i_user)))));
end
cap_tot=sum(Capacity_loop);
if abs(cap_tot-cap_tot_past)<tolerance
loop=0;
else
cap_tot_past=cap_tot;
count=count+1;
end
if count>=2000
break;
end
end
%Calculate the capacity of each UE
% trace(B(:,:,1)*B(:,:,1)‘)+trace(B(:,:,2)*B(:,:,2)’)
% SNR
for i_user=1:N_user
R(:,:,i_user)=eye(A_U);
for j_user=1:N_user
if j_user~=i_user
R(:,:,i_user)=R(:,:,i_user)+H(:,:,i_user)*B(:,:,j_user)*B(:,:,j_user)‘*H(:,:,i_user)’;
end
end
E(:,:,i_user)=inv(eye(A_U)+B(:,:,i_user)‘*H(:,:,i_user)’*inv(R(:,:,i_user))*H(:,:,i_user)*B(:,:,i_user));
Capacity(i_user)=real(log2(det(inv(E(:,:,i_user)))));
end
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function Capacity = DPC_rateRegion(weights,H,SNRdB,tolerance)
%------------------Algorithm proposed in Paper: Downlink Capacity
%Evaluation of Cellular Networks With Know-Interference Cancellation
%(Algorithm 1)
SNR = 10^(SNRdB/10);
[NT,A_U,N_user] = size(H) ;
%The algorithm is based on descend weights sequence of UEs
weights=weights(:);
[sortedWeights, index_UE]=sort(weights,‘descend’);
H = H(:,:,index_UE); %sorted H according to descend weights
%Step 1: Q initialization
Q = zeros(A_U,A_U,N_user);
gradient = zeros(A_U,A_U,N_user); %initialized gradient
loop=1;
cap_tot_past=0;
count=0;
%The algorithm iterates “Iterations” times to make Q converges
while(loop)
%covariance of each UE
for i_user=1:N_user
HQH(:,:,i_user)=H(:,:,i_user)*Q(:,:,i_user)*H(:,:,i_user)';
end
% Calculate the gradient of each UE on Q_i
for i_user=1:N_user
gradient_i_part = zeros(A_U,A_U);
for j_user=i_user:N_user-1
gradient_i_part = gradient_i_part+(sortedWeights(j_user)-sortedWeights(j_user+1))*...
(H(:,:,i_user)'*inv( eye(NT)+sum(HQH(:,:,1:j_user),3))*H(:,:,i_user));
end
gradient(:,:,i_user) = gradient_i_part + sortedWeights(N_user)*H(:,:,i_user)'*inv(eye(NT)+sum(HQH,3))*H(:,:,i_user);
end
%Step 2 & Step 3: principal eigenvector and the corresponding principal
%eigenvalues of the gradient
for i_user=1:N_user
[V,D]=eig(gradient(:,:,i_user));
[eigenValue(i_user) index_principal] = max(diag(D));
eigenVector(:,i_user)=V(:,index_principal);
end
%Step 4: optimal user selection
[~, i_opt_UE]=max(eigenValue);
% Step 5: calculating t, bisection or Matlab optimization toolbox
% Using MATLAB optimization toolbox, reference from 'J. Lee and N. Jindal,
%"Symmetric capacity of a downlink MIMO channel," IEEE ISIT 2006.'
obj_t = @(t) -obj_vvh(t,i_opt_UE,Q,sortedWeights,H,SNR,eigenVector(:,i_opt_UE));
t_opt=fminbnd(obj_t,0,1);
% Update Q
for i_user=1:N_user
if (i_user==i_opt_UE)
Q(:,:,i_user)=t_opt*Q(:,:,i_user)+(1-t_opt)*SNR*eigenVector(:,i_user)*eigenVector(:,i_user)';
else
Q(:,:,i_user)=t_opt*Q(:,:,i_user);
end
end
for i_user=1:N_user
info_UE(:,:,i_user) = eye(NT);
for j_user=1:i_user
info_UE(:,:,i_user) = info_UE(:,:,i_user) + H(:,:,j_user)*Q(:,:,j_user)*H(:,:,j_user)';
end
end
Capacity = zeros(N_user,1);
Capacity(index_UE(1)) = log2(real(det(info_UE(:,:,1))));
for i_user=2:N_user
Capacity(index_UE(i_user)) = log2(real(det(info_UE(:,:,i_user)))/real(det(info_UE(:,:,i_user-1))));
end
cap_tot_1=sum(Capacity);
if abs(cap_tot_1-cap_tot_past)<=tolerance
loop=0;
else
cap_tot_past=cap_tot_1;
count=count+1;
end
if count>=200
break;
end
end
% Calculate R
for i_user=1:N_user
info_UE(:,:,i_user) = eye(NT);
for j_user=1:i_user
info_UE(:,:,i_user) = info_UE(:,:,i_user) + H(:,:,j_user)*Q(:,:,j_user)*H(:,:,j_user)';
end
end
Capacity = zeros(N_user,1);
Capacity(index_UE(1)) = log2(real(det(info_UE(:,:,1))));
for i_user=2:N_user
Capacity(index_UE(i_user)) = log2(real(det(info_UE(:,:,i_user)))/real(det(info_UE(:,:,i_user-1))));
end
%
% Sub-function
%
function f=obj_vvh(t,i_opt_UE,Q,weights,H,SNR,v);
[NT,A_U,N_user] = size(H);
for i_user=1:N_user
if(i_user==i_opt_UE)
Q(:,:,i_user)=tQ(:,:,i_user)+(1-t)SNRvv’;
else
Q(:,:,i_user)=t*Q(:,:,i_user);
end;
end;
for i_user=1:N_user
HQH(:,:,i_user)=H(:,:,i_user)*Q(:,:,i_user)*H(:,:,i_user)';
end;
cap_part1=0;
for i_user=1:N_user-1
cap_part1=cap_part1+(weights(i_user)-weights(i_user+1))*log2(det(eye(NT)+sum(HQH(:,:,1:i_user),3)));
end;
cap_part2=weights(N_user)*log2(det(eye(NT)+sum(HQH(:,:,1:N_user),3)));
f=cap_part1+cap_part2;
end
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------