%% ECE 602 Course Project: Joint Trajectory and Communication Design for Multi-UAV Enabled Wireless Networks
% 
% *IEEE Transactions on Wireless Communications*
%
% *Team members* : Abdallah Mohamed Hussein(20779025), Ahmed Hussein Salamah (20800291)
%
% This report replicates, analyzes and produces results from the reseach paper, 
% titled: *"Joint Trajectory and Communication Design for Multi-UAV Enabled Wireless Networks"*.

%% Introduction
% Wireless communication shows a sustainable interest to use Unmanned Aerial Vehicle (UAVs) 
% to provide data coverage. This is mainly due to their flexible deployment, high maneuverability, 
% and significant low cost. 
% In this paper, the authors consider multi-UAVs to be used as base stations (BSs) 
% to provide coverage to a group of users from a certain altitude. The
% issue at hand is selecting the optimum passing trajectory of the UAVs and
% power and resource allocation/user association strategy for providing the best service to the
% users. In order to maintain fairness between the different users the authors 
% chose their maximization target to be the minimum throughput. 
% The UAVs cooperate to solve this joint problem by optimizing the scheduling 
% and association jointly with the UAVs’ trajectory and 
% power control to achieve more efficient communications using the different design 
% dimensions. The intially formulated optimization problem is non-convex.
%
% To find an  feasible solution for the given problem efficiently, the authors relaxed the binary 
% variables for user scheduling and association to be continuous variables and they 
% have used block coordinate decent method as an iterative algorithm to solve the resulting 
% problem.  The entire optimization problem is splitted into three partitions namely 
% the resourece allocation and user association, UAV trajectory, and transmitted 
% power control. Each variable is optimized in each iteration while fixing the other two 
% partitions. However, the iterative method did not solve the problem completely where 
% when fixed user scheduling and association, both UAV power control problem with fixed 
% trajectory case and its reverse, are still difficult to solve as due to their 
% non-convexity. Thus, the successive convex optimization technique is used to solve 
% this problem the underestimator taylor expansion of the function is used instead
% Additionaly, the author proposed a trajectory initialization scheme based 
% on the simple circular trajectory and the circle packing scheme. This scheme compared 
% to other benchmark schemes with heuristic UAV and conventional static UAV trajectories 
% showing a significant throughput gain and also proportionally with increasing 
% the UAV trajectory design period. 
% 

%% Problem Formulation
% let the user $k$ association with UAV $m$ at time instance $n$ be denoted as:
% $\mathbf{A} = \left \{\alpha_{k,m}[n],\forall k,m,n\right \}$, the
% trajectory of UAV $m$ at time instance $n$ be:
% $\mathbf{Q} = \left \{\mathbf{q}_{k,m}[n],\forall m,n\right \}$ and the power 
% transmitted by a UAV $m$ at time instance $n$ be: 
% $\mathbf{P} = \left \{\mathbf{p}_{k,m}[n],\forall m,n\right \}$ . 
% 
% The Goal is to maximize the minimum average rate among all users by jointly 
% optimizing the user scheduling and association (i.e., $A$), UAV trajectory 
% (i.e., $Q$), and transmit power (i.e., $P$) over all time slots.
%
%  
% $$\max_{\eta, \mathbf{A},\mathbf{Q},\mathbf{P}} \ \eta$$
%
% $\frac{1}{N}\sum_{n=1}^{N}\sum_{m=1}^{M}
% \alpha_{n,m}[n]\log_{2}(1+\gamma_{k,m}[n]])\geq \eta,\forall k$
%
% $\sum_{k=1}^{K} \alpha_{n,m}[n]\leq 1,\forall m,n$
%
% $\sum_{m=1}^{M} \alpha_{n,m}[n]\leq 1,\forall k,n$
%
% $\alpha_{n,m}[n]\in{0,1},\forall k,m,n$
%
% $\left \| q_{m}[n+1]-q_{m}[n]\right \|^2\leq S_{max}^{2},n=1,...,N-1$
%
% $q_{m}[1]=q_{m}[N],\forall m$
%
% $\left \| q_{m}[n]-q_{j}[n]\right \|^2\leq d_{min}^{2} \forall n,m,j\neq
% m$
% 
% $0\leq p_{m}[n]\leq P_{max}\forall m,n$
% 
% The first constraint is the received rate per each user which should be
% greater than the minimum rate $\eta$. The following 3 constraints
% describe the user association mechanisim. Only one user is associated at
% time per UAV so at any time instance $n$ both the summation over $m$ and
% $k$ of the association binary variable should be less than or equal 1. 
% The following constraints indicate that a UAV can take a maximum step of
% Smax which is a function of the max speed of the UAV. Also, the
% trajectory of any UAV should be periodic i.e. the start location is the
% same as the final point. Finally, the UAVs should always be at least
% separated by a distance $d_{min}$ at any time instance. The last 
% constraint is for the power control which simply states that the UAV
% power at any time should be between 0 and $P_{max}$. 

%% System parameters
% Defining the most important system parameters
%
%
% <html>
% <table border=1><tr><td>Parameter</td><td>Value</td></tr>
% <tr><td>Number of users (K)</td><td>6 users</td></tr>
% <tr><td>Channel power gain from UAV $m$ to user k</td><td>Free-Space Path Loss Model</td></tr>
% <tr><td>Noise power (sigma2)</td><td>-110 dBm</td></tr>
% <tr><td>UAVs Maximum Speed (V_max)</td><td>50 m/s</td></tr>
% <tr><td>Channel Power gain at distance 1m (rho)</td><td>-60 dB</td></tr>
% <tr><td>Maximum Power (Pmax)</td><td>0.1 W</td></tr>
% <tr><td>Minimum inter-UAV distance in m to ensure collision avoidance (dmin)</td><td>100 m</td></tr>
% <tr><td>Users distribution within 2D Area</td><td>2x2 km2</td></tr>
% <tr><td>Epsilon Maximum (NOt Given)</td><td>0.1</td></tr>
% <tr><td>Epsilon Thersould</td><td>10^(-4)</td></tr>
% <tr><td>Fixed altitude (H)</td><td>100 m</td></tr></table>
% </html>
%
% 
% The epsilon maximum $\varepsilon_{max}$ selection is based on several trials
% after which this value is capable of providing a good trade-off between
% results accuracy compared to results presneted by the authors and pc run-time.  

%% System Model
% The realization used in the paper for $K$ users  position $w_{max}=[x_{k}(t),y_{k}(t)]$
% is approximated to acheive the same results. All UAVs are assumed to fly 
% at a fixed altitude $H$ above ground and the time-varying horizontal 
% coordinate of UAV  time instant $t$ is denoted by
% $q_{m}(t)=[x_{m}(t),y_{m}(t)]^{T}$ with $0\leq t \leq T$
%
% Maximum Horizontal Distance
%
% $S_{max}=\varepsilon_{max} H$
%
% the minimum number of time slots required for achieving the accuracy 
% with a given $\varepsilon_{max}$ can be obtained as 
%
% $N\geq\frac{V_{max}T}{H\varepsilon_{max}}$
%
% Free-Space path loss model:
%
% $h_{k,m}[n]=\rho_{0} d_{k,m}^{-2}[n]= \frac{\rho_{0}}{H^{2}+ \left \|
% q_{m}[n]-w_{k}[n]\right \|^{2}}$
%
% Signal-to-interference-plus-noise ratio (SINR):
% 
% $\gamma_{k,m}[n] = \frac{p_{m}[n]H_{k,m}[n]}{\sum_{j=1,j\neq
% m}^{M}p_{j}[n]h_{k,j}[n]+\sigma^{2}}$

% Defining the most important system parameters

K = 6 ;       % Number of users
M = 1 ;       % Number of UAVs
Nr = 1  ;     % Number of realizations

L = 2e3 ;     % Coverage area length (m)
W = 2e3 ;     % Coverage area Width (m)
H = 100 ;     % UAV flight height (m)

sigma2 = 10^(-110/10)/1e3 ;% Noise power -110 dBm
do = 1 ;                   % refrence distance
rho = 10^(-60/10);         % channel power at refrence distance (-60 dBm)

Pmax = 0.1 ;               % UAV max power
Vmax = 50  ;               % UAV max speed (m/s)
dmin = 100 ;               % UAVs min separation distance

tau = 100;                 % time-slot factor
epislon_o = 1e-4;          % algorithim stopping criteria
epislon_max = 0.01;
Smax = H * epislon_max ;   % the maximum horizontal distance that the UAV can travel in each time slot

R_upperBound = (1/K)*log2(1+(rho*Pmax)/((H^2)*sigma2)); %max-min Rate for each user 


%% Proposed Solution
% *A. User Scheduling and Association Optimization*
% 
% For any given UAV trajectory and transmit power $\left\{\mathbf{Q},\mathbf{P}\right\}$ 
% ,the user scheduling and association of problem can be optimized by solving the following problem
%
%
% <<User Scheduling and Association.PNG>>
%
%
% *B. UAV Trajectory Optimization*
% 
% For any given user scheduling and association as well as UAV transmit power $\left\{\mathbf{A},\mathbf{P}\right\}$
% ,the UAV trajectory of problem can be optimized by solving the following problem
% 
% 
% <<UAV Trajectory Optimization.PNG>>
%
% where $\hat{R}_{k,m}^{lb}[n]$ is a convex function and considered as globally lower-bounded 
% by its first-order Taylor expansion at any point
%
% <<UAV Trajectory Optimization1.PNG>>
%
%
% <<UAV Trajectory Optimization2.PNG>>
% 
% *C. UAV Transmit Power Control*
%
% For any given user scheduling and association as well as UAV trajectory 
%  $\left\{\mathbf{A},\mathbf{Q}\right\}$ ,the UAV transmit power of problem can be 
% optimized by solving the following problem
%
% <<UAV Transmit Power Control.PNG>>
% 
% <<UAV Transmit Power Control1.PNG>>
%
% *D. Overall Algorithm and Convergence*
%
% As the problem is paritioned in to three problems that is solved iteratively 
% by fixing two partition and optimizing the remained. From
% the shown Algorithm, in step 3 the optimal solution will obtain
% $\eta(A^r,Q^r,P^r) \leq \eta(A^{r+1},Q^r,P^r)$. Then Algorithm will end 
% reaching $\eta(A^r,Q^r,P^r) \leq \eta(A^{r+1},Q^{r+1},P^{r+1})$ which
% indicates that the objected function is non-decearsing after each
% iteration which will be shown the analysis section. The stopping
% cariteria will be depend reaching a slight improving in the rate
% compared to a certain threshold value $\varepsilon_{max}$. 
%
% *Algorithm 1*
%
% <<Algo.PNG>>
%

%% *How to intiate the Code and Control the Algorithm for M = 1 & 2*
% In this section, we designed a way for the user to control the required
% case depending on the raised flags. As for M = 2, the power optimzation is only the case at 
% M=2 that's control the interference added on the other UAVs. So if we
% setted the power optimization to zero, it means that the algorithm will
% work with full power transmission. 
%
% Both F_run_M1 and F_run_M2 flags, they control the optimzation Algorithm
% to start as whole block. We have done this because it takes a long time 
% to converage, so we have saved the figures and data and we loaded them 
% in the code to save time and we seted both flags to be zero. we have 
% some test cases that are fully functional but it will take a long time to
% sweep all over muliple t values.
%
% *Flags to Control Code and UAV Motions and Select the optimization partition*
%

% Test Case "F_run_M1 = 1" and "F_ProposedTraj_UAV_M1 = 1"

% For M = 1
F_run_M1                = 0; %Flag to run the enable the next flags
F_static_UAV_M1         = 0; %Static
F_circular_UAV_M1       = 0; %Circular trajectory
F_ProposedTraj_UAV_M1   = 1; %Proposed trajectory

% Test Case "F_run_M2 = 1", "F_circular_UAVs  = 1" and 'F_PowerOpt = 1"
% Test Case "F_run_M2 = 1", "F_static_UAVs  = 1"   and 'F_PowerOpt = 1"

% For M = 2
F_run_M2                = 0; %Flag to run the enable the next flags
F_static_UAVs           = 0; %Static
F_circular_UAVs         = 1; %Circular trajectory
F_ProposedTraj_UAVs     = 0; %Proposed trajectory

F_PowerOpt              = 0; %Power Optimization
F_All                   = 0; % All Partitions (Whole Algorithm for M = 2)
%% Initialization for 1 UAV (M=1)
%

T = [30 60 120 210 300 400 800] ;
eta_A_M1 =[];
eta_A_t1 =[];
eta_T_M1 =[];
eta_T_t1 =[];
eta_All_M1=[];
eta_t1 = [];

for t = T
    N = ceil(Vmax*t/(Smax)) ;
    % generating the user's location realization approximated from the
    % paper
    w = [200 200 ; 600 1150;-400 400 ; ...
        -1100 600 ; -700 -600; 800 -400]; % realization used in the paper
    
    % Initialization (M=1)
    Po  = ones(N,1)*Pmax ;         % Initial power
    theta_n = 2*pi*([0:N-2 ,0])./(N-1) ;
    
    C = mean(w) ;                  % center for the trajectory
    R_max = Vmax*t/(2*pi) ;        % Maximum radius of the trajectory
    R_cg = max(sqrt(sum(((kron(C,ones(K,1))- w).^2),2))) ;
    R = min(R_cg,R_max) ;
    Qo = [ C(1)+R*cos(theta_n) ;  C(2)+R*sin(theta_n) ]' ;
    
    
    eta_o = 0 ;
    delta_eta = 1;
    
    P = Po ;
    r = 0 ;
    % Static position UAV
    if(F_static_UAV_M1 == 1)
        Q = [C(1)*ones(size(Qo,1),1) C(2)*ones(size(Qo,1),1)];
    else
        % Initial Circular Trajactory
        Q = Qo ;
    end
    %% Plot the initial trajectory for 1 UAV
%     figure(1)
%     plot(w(:,1) , w(:,2),'b o')
%     grid on
%     xlabel('x(m)')
%     ylabel('y(m)')
%     hold on
%     plot(Qo(:,1,1),Qo(:,2,1),'r -- ')
%     plot(C(1),C(2),'r x')
%     plot(Q(:,1,1),Q(:,2,1),'k --  ') 
    %% User Allocation Optimization Partition
    while (F_run_M1)
        [eta_Alloc_M1(r+1) , A, gamma] = optimize_Allocation_M1(Q,P,w,rho,H,sigma2,M,K,N);    % initialization
        
        delta_eta = eta_Alloc_M1(r+1)- eta_o;
        eta_o = eta_Alloc_M1(r+1);
        if delta_eta < epislon_o
            eta_A_M1 = [eta_A_M1;eta_Alloc_M1(r+1)];
            eta_A_t1 = [eta_A_t1;t];
            eta_All_M1 = [eta_All_M1;eta_Alloc_M1(r+1)];
            eta_t1 = [eta_t1;t];
            break
        end
    %% Trajectory Optimization Partition
         if (F_ProposedTraj_UAV_M1==1)
             [eta_Traj_underEst_M1(r+1),Q] = optimize_Trajectory_M_1(A,Q,P,N,K,sigma2,rho,H,w,Smax);
             eta_traj_M1(r+1) = min(Rate_Calculation(SINR_calculation(Q,P,w,rho,H,sigma2,M,K,N),A,K,M,N));
             delta_eta = eta_traj_M1(r+1)- eta_o;
             eta_o = eta_traj_M1(r+1);
             if delta_eta < epislon_o
                 eta_T_M1 = [eta_T_M1;eta_traj_M1(r+1)];
                 eta_T_t1 = [eta_T_t1;t];
                 eta_All_M1 = [eta_All_M1;eta_traj_M1(r+1)];
                 eta_t1 = [eta_t1;t];
                 break
             end
         end
     r = r + 1
    end
    clear Q Qo % Remove the generated Trajactory to prevent overwrite
end

if(F_run_M1)
    figure(1)
    hold on
    plot(Q(:,1,1),Q(:,2,1),'k --')
    
    figure(2)
    title('Max-min rate versus period T for a one UAV system with different optimization schemes')
    hold on
    plot(eta_t1,eta_All_M1,'-')
    if(F_circular_UAV_M1)
        legend('Circular trajectory_Allocation_Phase')
    elseif (F_static_UAV_M1)
        legend('Static UAV_Allocation_Phase')
    else
        legend('Proposed Trajectory')
    end
    % This figure show that the alogrithm converage to a higher rate
    % or it saturate at a certain value.
    figure(3)
    title('Compare Between Allocation and Trajectory Optization')
    subplot(1,2,1)
    plot(eta_A_t1,eta_A_M1,'-- r')
    legend('Allocation Opt')
    subplot(1,2,2)
    plot(eta_T_t1,eta_T_M1,'-- b')
    legend('Trajectory Opt')
end

    %% Upper Bound Cal. for M = 1
%
% figure(2)
%
% hold on
% 
% plot(T,R_upperBound*ones(T,1),'-- r') % the upper bound 
% 
% legend('Upper Bound')

%% Allocation and Trajectory Optimization
%
% In this part, we plot the implementation of the whole alogrithm for 
% previous time instances. It is shown that the algorithm works as 
% described in the paper  as it depends on the shown trajeactory as a
% result of the given time. For example, for T = 30 the UAV does not have time to
% cover the whole users. As seen at T = 120, the UAV cares for outer users as the inner 
% user will be explicitly covered and for T = 210 the UAV have sufficient time to cover
% all users. We have run two experment, the next shown is the first
% experiemnt.
%
%
% <<paperRef1.PNG>>
%
%
%

%
% plot_results_M1() %this function is used to plot the saved .mat files for the
% down mentioned figures. 
%


% from figure 4 to 10 (Paper Figure 3)
uiopen('fig3_30.fig',1)
uiopen('fig3_60.fig',1)
uiopen('fig3_120.fig',1)
uiopen('fig3_210.fig',1)
uiopen('fig3_300.fig',1)
uiopen('fig3_400.fig',1)
uiopen('fig3_800.fig',1)

%% The Second Experiment trajectory at M = 1

uiopen('fig_T1_M_1.fig',1)


%% Analysis and Conclusions for Single UAV
%
% In this figure, we compare the average max-min rate between three schemes, 
% which are :
% 
% * Proposed trajectory, which is the result of the proposed Algorithm.
% * Circular trajectory, which is the result of the initialization scheme.
% * Static UAV, which the UAV is located at the center without any motion
%  
% All schemes were optimizaed using Algorithm 1 given the previous
% trajectories., which shows that the max-min rate of the static UAV is
% invariant witt increasing time since there UAV is in static position
% without any mobility. In construct, the proposed and circular trajectory 
% schemes increase the max-min rate with time and ultimately becomes 
% saturated when T becomes sufficiently large to almost near to the upper bound.
% This behavior is expected since with the UAV mobility, a larger value 
% of T provides the UAV more time to fly near to the users to be served, 
% which will improve the max-min rate. Both figures for different
% experiement that reveal for the same behaviour. 

uiopen('Fig5.fig',1) % figure 11 (Paper Figure 5)
uiopen('fig_M_1.fig',1)
%% Initialization of Circular trajectory for 2 UAVs (M=2)
% 
% In the shown figure, we show the Circular trajectory for 2 UAVs serving different areas, each
% circle depends on the ground servered users position and its radius depends on 
% the radius of the larger circle that will cover the last user. The small
% circle indicates that the UAV assigned a short time to cover the servered
% users in its region. Then, by scaling up the coverage circle area, which mean 
% the UAV will have longer time to cover the surounded users. This Circular
% trajectory will be selected initial to start the algorithm as done before in M =
% 1 just to start the allocation optimization, then the trajectory
% optimization start after. If the  trajectory optimization is not
% selected as a flag the the circular trajectory will be always used.


uiopen('2UAVs_Motion.fig',1) % figure 12


%%
T = [5 30 60 90 200 300 400 500 600] ;
eta_All_M2 = [];
eta_t2 = [];

eta_A = [];
eta_A_t2 = [];
eta_T2 = [];
eta_T_t2 = [];
eta_P2 = [];
eta_P_t2 = [];


% 
% In this part the T is sweeped all over different values of t to generate
% different trajectories as shown before and then comapre between the
% reproduced rates. 
%


% Defining M for 2 UAVs for another users position compared to M =1. We 
% relaized an approzimated positions for the user as much as possible form
% the paper.
M = 2 ;       % Number of UAVs

for t = T
    N = ceil(Vmax*t/(H*epislon_max)) ;
    delta_t = t/N;
%     Smax = delta_t * Vmax ;
%     Smax = epislon_max * H;
%     if ((Smax/H)>epislon_max)
%         Smax = epislon_max * H;
%     end
    

    % generating the user's location realization
    w = [-1200 570 ; -450 400;-700 -600 ; ...
        600 1100 ; 150 200;750 -400]; % realization used in the paper

    % Initialization (M=1)
    Po  = ones(N,M)*Pmax ;         % Initial power
    theta_n = 2*pi*([0:N-2 ,0])./(N-1) ;

    C = mean(w) ;                  % center for the trajectory
    R_max = Vmax*t/(2*pi) ;        % Maximum radius of the trajectory
    R_cg = max(sqrt(sum(((kron(C,ones(K,1))- w).^2),2))) ;
    R = min(R_cg,R_max) ;
    Qo(:,:,1)  = [ C(1)+R*cos(theta_n) ;  C(2)+R*sin(theta_n) ]' ;
    Qo(:,:,2)  = [ C(1)+R_max*cos(theta_n) ;  C(2)+R_max*sin(theta_n) ]' ;
    C_1 = mean(w(4:6,:));
%     C_1 = [200 210];
    R_1 = R/2;
    Qo(:,:,3)  = [ C_1(1)+R_1*cos(theta_n) ;  C_1(2)+R_1*sin(theta_n) ]' ;
    Qo(:,:,4)  = [ C_1(1)+(R_1/2)*cos(theta_n) ;  C_1(2)+(R_1/2)*sin(theta_n) ]' ;
    
    C_2 = mean(w(1:3,:));
%     C_2 = [-600 180];
    R_2 = R/2;
    Qo(:,:,5)  = [ C_2(1)+R_2*cos(theta_n) ;  C_2(2)+R_2*sin(theta_n) ]' ;
    Qo(:,:,6)  = [ C_2(1)+(R_2/2)*cos(theta_n) ;  C_2(2)+(R_2/2)*sin(theta_n) ]' ;

    eta_o = 0 ;
    delta_eta = 1;
    P = Po ;
    r = 0 ;

    % Static position for 2 UAVs
    if(F_static_UAVs == 1)
         Q(:,:,1) = [ones(size(Qo(:,:,1),1),1)*C_1(1) ones(size(Qo(:,:,1),1),1)*C_1(2)];
         Q(:,:,2) = [ones(size(Qo(:,:,1),1),1)*C_2(1) ones(size(Qo(:,:,1),1),1)*C_2(2)];
    else
    % Initial Circular Trajactory of 2 UAVs
        Q(:,:,1) = Qo(:,:,4);
        Q(:,:,2) = Qo(:,:,6);
    end

%% Plot the initial trajectory for Both UAVs  
%     figure(12)
%     title('Users and UAVs Trajectory')
%     grid on
%     hold on
%     plot(w(:,1) , w(:,2),'b o')
%     xlabel('x(m)')
%     ylabel('y(m)')
%     plot(Qo(:,1,1),Qo(:,2,1),'r -- ')
%     plot(Qo(:,1,2),Qo(:,2,2),'r - ')  
%     plot(Qo(:,1,3),Qo(:,2,3),'b -- ')
%     plot(Qo(:,1,4),Qo(:,2,4),'b - ')
%     plot(Qo(:,1,5),Qo(:,2,5),'-- ', 'color',[0.21,0.30,0.08])
%     plot(Qo(:,1,6),Qo(:,2,6),' - ', 'color',[0.21,0.30,0.08])
%     C = mean(w);
%     plot(C(1),C(2),'r x')
%     plot(C_1(1),C_1(2),'b x')
%     plot(C_2(1),C_2(2),'x', 'color',[0.21,0.30,0.08])
%     plot([C_1(1) C_2(1)],[C_1(2) C_2(2)],'b --')
%     % Used Trajectory in the Current Run
%     figure(13)
%     hold on
%     grid on
%     xlabel('x(m)')
%     ylabel('y(m)')
%     title('Current Runing Trajectory')
%     plot(Q(:,1,1),Q(:,2,1),'r * ')
%     plot(Q(:,1,2),Q(:,2,2),'b *')



    %% User Allocation Optimization Partition
    % In the Users Allocation allocation part the resulted rate is the
    % actual value where we did not used any type of approxiation to
    % manipulate with the non-convexity of the function 
    while (F_run_M2)
        [eta_Alloc(r+1) , A, gamma] = optimize_Allocation_M2(Q,P,w,rho,H,sigma2,M,K,N);    % initialization
         delta_eta = eta_Alloc(r+1)- eta_o;
         eta_o = eta_Alloc(r+1);
         if delta_eta < epislon_o
             eta_A = [eta_A;eta_Alloc(r+1)];
             eta_A_t2 = [eta_A_t2;t]
             eta_All_M2 = [eta_All_M2;eta_Alloc(r+1)];
             eta_t2 = [eta_t2;t];
             break
         end
         
    %% Trajectory Optimization Partition
    %
    % The CVX could not solve this problem as the auothers mentioned. The
    % reason that $\hat{R}_{k,m}^lb[n]$ is concave and we are subtracting $\log
    % \left ( \sum_{j=1,j\neq m}^{M} \frac{p_{j}[n]  \rho_{0}}{H^{2}+S_{k,j}[n]} + \sigma^{2}  \right )$
    % ,which is log of a convex function not a concave that prevents CVX from to find its maximum.
    % We have used another solver built in the matlab instead of CVX, but it 
    % does not improve the throughput and not relevant to the authors results. 
    % In the down mentioned function it is written for the CVX, the other
    % solver is mention at the end of the script.
         if (F_ProposedTraj_UAVs ==1 || F_All)
             % eta_Traj_underEst is the approximated Rate after Maximiizing
             % but not the actual rate. The resulted Q will be used for
             % calucating the actual rate.
             [eta_Traj_underEst(r+1),Q] = optimize_Trajectory_M_2(A,Q,P,N,K,M,sigma2,rho,H,w,Smax,dmin);
              eta_traj(r+1) = min(Rate_Calculation(SINR_calculation(Q,P,w,rho,H,sigma2,M,K,N),A,K,M,N))
              delta_eta = eta_traj(r+1)- eta_o;
              eta_o = eta_traj(r+1);
              if delta_eta < epislon_o
                  eta_T2 = [eta_T2;eta_traj(r+1)];
                  eta_T_t2= [eta_T_t2;t];
                  eta_All_M2 = [eta_All_M2;eta_traj(r+1)];
                  eta_t2 = [eta_t2;t];
                  break
              end
         end
    %% Users Power Optimization Partition         
        if (F_PowerOpt ==1 || F_All)
            [eta_Pow_underEst(r+1),P_r1] = power_optimization(Pmax,A,Q,P,M,N,K,rho,w,H,sigma2);
            eta_pow(r+1) = min(Rate_Calculation(SINR_calculation(Q,P,w,rho,H,sigma2,M,K,N),A,K,M,N))
            delta_eta = eta_pow(r+1)- eta_o;
            eta_o = eta_pow(r+1);
            if delta_eta < epislon_o
                eta_P2 = [eta_P2;eta_pow(r+1)];
                eta_P_t2 = [eta_P_t2;t];
                eta_All_M2 = [eta_All_M2;eta_pow(r+1)];
                eta_t2 = [eta_t2;t];
                break
            end
        end
         r = r + 1
    end
        clear Qo Q % Remove the generated Trajactory to prevent overwrite
end
if (F_run_M2)
    figure(15)
    title('Max-min rate versus period T for a two-UAV system with different optimization schemes')
    hold on
    plot(eta_t2,eta_All_M2,'-')
    if (F_circular_UAVs || F_static_UAVs)
        if(F_circular_UAVs)
            legend('Circular trajectory/Scheme III')
        elseif (F_static_UAVs)
            legend('Static UAV')
        end
        figure(16)
        title('Max-min rate versus period T for a two-UAV system after Allocation Optimization')
        hold on
        plot(eta_A_t2,eta_A,'-')
        if(F_circular_UAVs)
            legend('Circular trajectory/Scheme III')
        elseif (F_static_UAVs)
            legend('Static UAV')
        end
    end
    
    
    if (F_ProposedTraj_UAVs || F_All)
        figure(13)
        hold on
        plot(Q(:,1,1),Q(:,2,1),'r - ')
        plot(Q(:,1,2),Q(:,2,2),'b - ')
        figure()
        hold on
        plot(eta_T_t2,eta_T2,'-')
        legend('Proposed Trajectory')
    end
    
    if (F_PowerOpt || F_All)
        figure(17)
        title('UAV transmit power versus time for a two-UAV system after Power Optimization')
        plot(eta_P_t2, eta_P2,'-')
        
    end
end


%% Analysis and Conclusions for Multiple UAVs
% 
% In figure 8, we implement the 3rd scheme in the paper which only depend
% on the user allocation and transmitting full power. The other two schemes 
% mentioned as they depend on the trajectory optimization optimization. If
% we tried to optimimize the user allocation and power the and fixing the
% circular trajector, the result rate from the power stage will be the same
% as the allocation rate, which shows the importance of trajectory optimzation.
% This proof that the optimzation algorithm converage in the right 
% direction as mentioned before to be $\eta(A^r,Q^r,P^r) \leq \eta(A^{r+1},Q^r,P^r)$. 
% From figure 8, we also proofed that max-min rates of the 3rd scheme increase with the period T
% becomes larger. 
% 
% In figure 11, the average max-min rate achieved for 2 UAVs by two trajectory, which
% are the circular trajectory and static UAV, where it is fixed at its own
% center. It is shown that without any trajectory optimzation there will be
% no imporvement in the average rate. As shown, when the UAV has enough
% time to travel the optimal average rate is almost 1.1 bps and then it got
% saturated as the interference between the two UAVs can be further be
% optimizated. 

uiopen('Fig8.fig',1) % (Paper Figure 8)
uiopen('Fig11.fig',1)% (Paper Figure 11)

%% Reference
% 
% # Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication
% design for multi-uav enabled wireless networks,” IEEE Transactions
% on Wireless Communications, Early access, 2018.
% # Packings of Equal Circles in Fixed-Sized Containers With Maximum
% Packing Density. Accessed: Apr. 5, 2017. [Online]. Available:
% <http://www.packomania.com>

%% Solving the Trajectory partition using other solvers
% In order to solve the given problem for trajectory optimization we have
% tried all the below solvers with no luck: 
%%
% 
% * *CVX solver:* we have overrided the DCP rule that did not allow a
% subtraction of a concave and a convex term which results from the LHS log
% term in the inequality constraint. This is done using the method
% described in <http://cvxr.com/cvx/doc/advanced.html#new-functions-via-partially-specified-problems>
% where we defined a trivial maximization problem that will result in the
% required concave expression. This expression is regarded as affine by the
% CVX and thus it allowed the CVX to continue. Afterwards, we tried both
% solver offered by the CVX free edition namely, _SeDuMi_ and _SDPT3_ sadly
% without any success "the output of the function was an infeasible which
% was the origin point at all times. This solution is infeasible because of
% the constraint of the minimum distance separation $d_{min}$.
% solution. 
% * *MATLAB Solvers* when CVX failed to solve the given problem we resorted
% to solving it via MATLAB solvers where the function was implemented again
% from scratch and passed to the global search solver. despite claiming
% that the solver is capable of providing a global minimizer/maximizer
% value for the given function. It failed to deliver any solution other
% than the initial trajectory providied as an initial point. which was
% better that what CVX could do. 
% 
% The functions used for implementing this functionality is provided as well in the below section  

% Matlab Solver
uiopen('FigM2_Matlab_Solver.fig',1)

uiopen('M_2_30.fig',1)

uiopen('M_2_60.fig',1)

uiopen('M_2_90.fig',1)



type('optimize_Trajectory_M_2_Using_MATLAB_Solver')

type('constraints_for_MATLAB')
%%
% CVX

type('main_M_2')

type('optimize_Trajectory_M_2_1')

type('rate_convex_trajectory_M_2_1')

type('partial_definition_cvx')

%% Print all the functions used in the project
% Common functions between M = 1 and 2 cases
type('SINR_calculation')

type('Rate_Calculation')

% *Used Function for M = 1*

% Allocation Functions
type('optimize_Allocation_M1')

% Trajectory Functions

type('optimize_Trajectory_M_1')

type('rate_convex_trajectory_M1')


% *Used Functions for M = 2*

% Allocation Functions
type('optimize_Allocation_M2')

% Trajectory Functions
type('optimize_Trajectory_M_2')

type('rate_convex_trajectory_M_2')

% Power Optimziation

type('power_optimization')

type('Rate_function_approx_for_power_alloc')



