clear all
%%%%%%initialize variables
tolerance=1e-3;
N=5; % number of users
Rate_DL_min=100e3;
Rate_UL_min= 100e3;
go = 1.42e-2;
Go=2.2846;
phi_min=0.2;
phi_max=1.2;
H_min=50;
H_max=200;
No = 1.2589254118e-20; 
DTx=8; %tx antennas, receiver antenna is 1
sigma_D= No*rand(1); % noise at drone BS
sigma_n= No*rand(N,1);
W=10e6;
Wk_old=0;
H_old=0;
phi_old=0;
%alpha = 0.5;



%%%%%%%%%%%% user-generate %%%%%%%%%%%%%%%%%%%%%%%

radius = 100;



% Data
% users generate
xc = randn;
yc = randn;
Qx=xc; % BS-drone initial location 
Qy=yc;

theta = rand(1,N)*(2*pi);
r = sqrt(rand(1,N))*radius;
X = xc + r.*cos(theta);
Y = yc + r.*sin(theta);
D=sqrt((X-xc).^2+(Y-yc).^2);
D_max= max(D);


P_D = 0.5:0.3:2.5;
%alpha=0.1:0.1:0.9;
alpha=0.5;
%P_D=1.2; % P_D 1.2 for multuple alphas
%P_D=2.5; %for multiple P_B
%P_B = 1e-5:0.5e-5:0.4e-4 ; % uplink user max power
P_B=1e-4; % uplink user max power_ used for all the previous data set
%P_D=2.5;
%loop control manually
for y=1:1:size(P_B,2) %multiple P_B
for u=1:1:size(alpha,2) %multiple alpha
for e=1:1:size(P_D,2) %multiple P_D
    
%%%%lagrange multipliers initial values
mu_n=0.1*ones(N,1);
nu_n=0.1*ones(N,1);
zeta_n=0.1*ones(N,1);
xi=1;

%optim param init
H=50;
phi=0.7;

p_min=1e-5; %1e-5 for PD
p_n=1e-5*ones(N,1);

%P_D= 1; % downlink max power
%P_B=1e-4; % uplink user max power_ used for all the previous data set

%%%%%%%%%%%%%%%%%%%channels generate%%
vari= No;
h_hat=  rand(DTx,1); % random variable for channel
 
%h_hat = 1/sqrt(2)*[randn(DTx,1) + j*randn(DTx,1)]; % Rayleigh channel

%%%SI channel on drone and user including CC%%
h_D= sqrt(vari)*(rand(DTx)); %%SIC on FD MIMO BS

for i=1:N
    for j=1:N
        if i==j
            g_mn(i,j)=sqrt(vari)*(0.1143); %%SIC on fd user %fixing rand(1)
            %g_mn(i,j)=0;
        else
            g_mn(i,j)= go/(sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2));     % channel between users
        end
    end
end

%%%%%%%%%CCI+SI on the ground users (sum)
for i=1:N
    for j=1:N
    CCI_SIU(i,j) = g_mn(i,j)^2 * p_n(i);
    end
end
CCI_SIU=sum(CCI_SIU,'all');

%%initializing tansmit beamforming k_n 
for(p=1:1:N)
k_n(:,p)= (P_D(e)/N)*rand(DTx,1); 
end

%%%overwrite all random parameters
%%% Following two lines loaded user locations and channels from a file for comparison purpose
%load('users_drone_loc_andetc.mat')
%load('channel_bf.mat')


%%%%%initliazing channel BS-user and transmit beamforming vector k_n 
t_H=2;
t_phi=2;

H_opt=0; phi_opt=0; 

while (t_H>0.5 || t_phi>0.001)

for(p=1:1:N)
h_n(:,p) = h_hat./((X(p)-Qx)^2+(Y(p)-Qy)^2+H^2); % (DTx,user) 
H_n(:,:,p) = h_n(:,p)*h_n(:,p)'; %% H= h*h^H

end

%%%%%%%%%receive beamforming w_n solution ZF-BF%%%%%%%

for i=1:N
    z=[zeros(i-1,1); 1; zeros(N-i,1)]';
    %w_n(:,i)=((inv(h_n'*h_n)*h_n'))/norm((inv(h_n'*h_n)*h_n'))
    w_n(:,i)=(z*(pinv(h_n'*h_n)*h_n'))'/norm(z*(pinv(h_n'*h_n)*h_n'))';
    
    W_n(:,:,i) = w_n(:,i)* w_n(:,i)' ; % Wn=wn*wn'
    


end




%%%%%%%%%%%%%%%%%%%% optimize H and Phi using matlab optimizer (fmincon)
%y(1)= H and y(2) = phi

fun=@(y)- sum(arrayfun(@(l)W*alpha(u)* log2(1+ ( (h_hat/((X(l)-Qx)^2+(Y(l)-Qy)^2+y(1)^2))' *w_n(:,l)*sqrt(p_n(l)))^2 / (W*sigma_D^2) )...
    +W*(1-alpha(u))* log2(1+ ( (Go/(y(2)^2)*( (h_hat/((X(l)-Qx)^2+(Y(l)-Qy)^2+y(1)^2))' * k_n(:,l))^2  )/(CCI_SIU+ W*sigma_n(l)^2))) ,1:N));

c= @(y) horzcat(arrayfun(@(l) ((X(l)-Qx)^2+(Y(l)-Qy)^2)- y(1)^2*(tan(y(2)))^2 ,1:N) , ...
    arrayfun(@(l) Rate_UL_min - W*alpha(u)* log2(1+ ((h_hat/((X(l)-Qx)^2+(Y(l)-Qy)^2+y(1)^2))'*w_n(:,l)*sqrt(p_n(l)))^2 / (W*sigma_D^2) ),1:N)...
    ,arrayfun(@(l) Rate_DL_min - W*(1-alpha(u))* log2(1+ ( (Go/y(2)^2)*((h_hat/((X(l)-Qx)^2+(Y(l)-Qy)^2+y(1)^2))' * k_n(:,l))^2  )/(CCI_SIU+ W*sigma_n(l)^2)),1:N));
%c2= @(y) arrayfun(@(l) Rate_UL_min - W*alpha(u)* log(1+ ((h_hat/((X(l)-Qx)^2+(Y(l)-Qy)^2+y(1)^2))'*w_n(:,i)*sqrt(p_n(i)))^2 / (W*sigma_D^2) )/log(2),1:N);
ceq = @(y)[];
nonl_con = @(y)deal(c(y),ceq(y));
gs = GlobalSearch('Display','off');
%opts = optimoptions('fmincon','MaxFunEvals', 10000,'TolX',1e-4,'TolCon',1e-4);
opts = optimoptions(@fmincon,'Display','none','Algorithm','interior-point');
%problem = createOptimProblem('fmincon','x0',[0.1],'objective',fun,'lb',[0.02],'ub',[0.9],'options',opts);
problem = createOptimProblem('fmincon','x0',[H_min phi_min],'objective',fun,'lb',[H_min phi_min],'ub',[H_max phi_max],'nonlcon',nonl_con,'options',opts);

H_and_phi=run(gs,problem);
H=H_and_phi(1);
phi=H_and_phi(2);

%%%%%%%%%%%%%%%%%% optimize Q and Pt using ALM%%%%%%%
%y(1)=Qx, y(2)=Qy and y(3) to y(l+3) p_n
t_Qx=2;
t_Qy=2;
t_p_n=ones(1,N);

Qx_opt=0; Qy_opt=0; p_n_opt =zeros(1,N);
while (t_Qx>0.5 || t_Qy>0.5 || all(t_p_n>1e-7))

    
    
    
%%%%%%%%%%%%%%%%%%%%%%%transmit beamforming k_n solution using SDP
    sum2=0;
cvx_begin sdp quiet
variable  K(DTx,DTx,N) 
 sum1=0;
    for i = 1:N
        %sum1 = sum1 + W*(1-alpha(u))*log(1+ ((Go/phi^2)*trace(H_n(:,:,i)*K(:,:,i)))/(CCI_SIU+W*sigma_n(i)^2))/log(2);
        sum1 = sum1 + trace(H_n(:,:,i)*K(:,:,i));
    end
maximize sum1
subject to
    %for i = 1:N
      %W* (1-alpha(u))*log(1+((Go/phi^2)*trace(H_n(:,:,i)*K(:,:,i)))/(CCI_SIU+W*sigma_n(i)^2))/log(2) >= Rate_DL_min; %CVX issue, taking always higher   
    %end

    
    for i = 1:N
        trace(K(:,:,i)) == P_D(e)/N;
        K(:,:,i)== semidefinite(DTx);
        %rank(K(:,:,i))==1;
    end
    for i=1:N
        trace(h_D*K(:,:,i)*h_D') == 0;
    end
       for i=1:N
           for j=1:N
               if j==i 
               else
                  trace(H_n(:,:,i)*K(:,:,j)) == 0;  
               end
           end  
       end
         
cvx_end


for i=1:N
    S=eig(K(:,:,i));
    %[A,B]=eigs(K(:,:,i),1);
    k_n(:,i)= S;
end


    
fun2=@(y)- (sum(arrayfun(@(l)W*alpha(u)* log2(1+ ( (h_hat/((X(l)-y(1))^2+(Y(l)-y(2))^2+H^2))' *w_n(:,l)*sqrt(y(l+2)))^2 / (W*sigma_D^2) )...
    +W*(1-alpha(u))* log2(1+ ( (Go/(phi^2)*( (h_hat/((X(l)-y(1))^2+(Y(l)-y(2))^2+H^2))' * k_n(:,l))^2  )/(sum(arrayfun(@(m) sum(arrayfun(@(l) g_mn(l,m)^2 * y(l+2),1:N)),1:N))+ W*sigma_n(l)^2))) ,1:N)) ...
    + (1/(2*xi))*(sum(arrayfun(@(l)(max(0, mu_n(l) - xi*( Rate_UL_min -(W*alpha(u)* log2(1+ ( (h_hat/((X(l)-y(1))^2+(Y(l)-y(2))^2+H^2))' *w_n(:,l)*sqrt(y(l+2)))^2 / (W*sigma_D^2) ))))^2 - mu_n(l)^2),1:N ) ) ...
    + sum(arrayfun(@(l)(max(0, nu_n(l) - xi*( Rate_DL_min-(W*(1-alpha(u))* log2(1+ ( (Go/(phi^2)*( (h_hat/((X(l)-y(1))^2+(Y(l)-y(2))^2+H^2))' * k_n(:,l))^2 )/( sum(arrayfun(@(m) sum(arrayfun(@(l) g_mn(l,m)^2 * y(l+2),1:N)),1:N)) + W*sigma_n(l)^2))))))^2 - nu_n(l)^2),1:N ) ) ...
    + sum(arrayfun(@(l) (max(0,zeta_n(l) - xi*( -(X(l)-y(1))^2+(Y(l)-y(2))^2+H^2*(tan(phi))^2))^2) - zeta_n(l)^2 ,1:N))));

 %+ sum(arrayfun(@(l) (max(0,iota_n(l)+xi*(p_n(l)-P_B(y)))^2-iota_n^2) ,1:N))

d = [];
gs = GlobalSearch('Display','off');
%opts = optimoptions('fmincon','MaxFunEvals', 10000,'TolX',1e-15,'TolCon',1e-15);
opts = optimoptions(@fmincon,'Display','none','Algorithm','interior-point');
problem = createOptimProblem('fmincon','x0',[1 1 p_min*ones(1,N)],'objective',fun2,'lb',[-inf -inf p_min*zeros(1,N)],'ub',[+inf +inf P_B(y)*ones(1,N)],'options',opts);
x=run(gs,problem);
Qx=x(1);
Qy=x(2);
p_n=x(3:N+2);

%%%updating CCI+SIU with p_n
for i=1:N
    for j=1:N
    CCI_SIU(i,j) = g_mn(i,j)^2 * p_n(i);
    end
end
CCI_SIU=sum(CCI_SIU,'all');

%%%%%%%update lagrange multipliers
for i=1:N
mu_n(i) = max(0, mu_n(i) - xi*(Rate_UL_min - (W*alpha(u)* log2(1+ ( (h_hat/((X(i)-Qx)^2+(Y(i)-Qy)^2+H^2))' *w_n(:,i)*sqrt(p_n(i)))^2 / (W*sigma_D^2) ))));
nu_n(i) = max(0, nu_n(i) - xi* (Rate_DL_min-(W*(1-alpha(u))* log2(1+ ( (Go/(phi^2)*( (h_hat/((X(i)-Qx)^2+(Y(i)-Qy)^2+H^2))' * k_n(:,i))^2 )/( CCI_SIU + W*sigma_n(i)^2))))));
zeta_n(i)=max(0,zeta_n(i)- xi*( -(X(i)-Qx)^2+(Y(i)-Qy)^2+ H^2*(tan(phi))^2));
%iota_n(i)= max(0, iota_n(i) + xi*(p_n(i)-P_B(y)) ); %matlab optimizer handles

end
xi=xi+1;
t_Qx = abs(abs(Qx)-abs(Qx_opt));
t_Qy = abs(abs(Qy)-abs(Qy_opt));
t_p_n=abs(p_n-p_n_opt);

Qx_opt=Qx; Qy_opt=Qy; p_n_opt=p_n;

end

t_H=abs(H-H_opt);
t_phi=abs(phi-phi_opt);

H_opt=H; phi_opt=phi;

end

for i=1:N
    FD_Rate_n(i)= W*alpha(u)* log2(1+ ( (h_hat/((X(i)-Qx)^2+(Y(i)-Qy)^2+H^2))' *w_n(:,i)*sqrt(p_n(i)))^2 / (W*sigma_D^2) )...
    +W*(1-alpha(u))* log2(1+ ( (Go/(phi^2)*( (h_hat/((X(i)-Qx)^2+(Y(i)-Qy)^2+H^2))' * k_n(:,i))^2  )/(CCI_SIU+ W*sigma_n(i)^2)));
    HD_DL_Rate_n(i)= W*(1-alpha(u))* log2(1+ ( (Go/(phi^2)*( (h_hat/((X(i)-Qx)^2+(Y(i)-Qy)^2+H^2))' * k_n(:,i))^2  )/(CCI_SIU+ W*sigma_n(i)^2)));
    HD_UL_Rate_n(i)= W*alpha(u)* log2(1+ ( (h_hat/((X(i)-Qx)^2+(Y(i)-Qy)^2+H^2))' *w_n(:,i)*sqrt(p_n(i)))^2 / (W*sigma_D^2) );
end


Sum_UL_power(e)=sum(p_n); %select between y, u and e to variate P_B,alpha, P_D, respectively.
Sum_FD_Rate(e)= sum(FD_Rate_n)
Sum_HD_DL_Rate(e)= sum(HD_DL_Rate_n);
Sum_HD_UL_Rate(e)= sum(HD_UL_Rate_n);

end
end
end
 
