Relay

clear all
close all
clc


tic

global A_new B N_relay K xb yb v_max_mine delta_t h_b aggregators H F C1 C2 B_c x_cp y_cp init


%% 1. kmeans clustering - Define aggregators positions
% *************************************************************************
% *************************************************************************
% *************************************************************************
% *************************************************************************

% Define the number of sensors and their positions
LEN_X = 100;
LEN_Y = LEN_X;
K = 5;
x_aggr = -LEN_X/2 + rand(1, K) * LEN_X;
y_aggr = -LEN_Y/2 + rand(1, K) * LEN_Y;
aggregators = [x_aggr; y_aggr];

points = [x_aggr; y_aggr];
colors1 = ["red", "green", "blue", "magenta", "cyan", "yellow", "black"];
colors2 = zeros(K, 3);
for k = 1 : K
    colors2(k, :) = [rand, rand, rand];
end

colors = [];
if K <= length(colors1)
    colors = colors1;
else
    colors = colors2;
end

% aggregators = [20 10 0 -20 5;15 5 10 -30 -30];


%% 2. Fixed data about position and frequency
% *************************************************************************
% *************************************************************************
% *************************************************************************
% *************************************************************************

% BS === CP
xb = 1 * LEN_X / 2;    yb = 0;    h_b = 15;    qb = [xb, yb];
x_cp = xb;    y_cp = yb;    q_cp = [x_cp, y_cp];

% Data EM waves
f = 2.4e9;    c = 3e8;    lambda = c / f;    B = 1e6;
Gt = 1;    Gr = 1;    G = Gt*Gr;
Gt_dB = natural_to_dB(Gt);    Gr_dB = natural_to_dB(Gr);
P_tr = 0.001;                   % 0 dBm
P_ant = 0.001;                  % 0 dBm

% Formulas for noise
sigma_dBm = -174 + 10 * log10(B);
sigma_dB = sigma_dBm - 30;
P_noise = dB_to_natural(sigma_dB);

% Calculations
Gt_dB = natural_to_dB(Gt);    Gr_dB = natural_to_dB(Gr);
% FSPL
d0 = 1;
C0_1m_dB = - (20 * log10(d0 * f) + 20 * log10(4*pi/c) - Gt_dB - Gr_dB);
C0_1m = dB_to_natural(C0_1m_dB);

% Data of UAV-RIS
H = 30;
d = lambda / 10;
delta_t = 1;
v_max_kmh = 62;
v_max = v_max_kmh / 3.6;                          % 16.67 m/s

% SNR Threshold
SNR_thr = 1;
kappa_dB = 100;
kappa = 10^(kappa_dB/10);




%% 3. UAV Components - Find flight time
% *************************************************************************
% *************************************************************************
% *************************************************************************
% *************************************************************************

% UAV Frame = 50cm x 50cm
desired_side = 0.5;                 % RIS Frame Side
% BATTERY
name = "Tattu 30.000 mAh";
B_c = 180 * 3600;               % Wh  ---> W * s = Joule
B_w = 1.35;                     % kgr
% WEIGHTS
UAV_w = 3.25;                   % kgr
Ant_w = 8 / 1000;               % WiFi antenna = 8 gr weight
T_max = 20;                     % kgr
% Relay
a_PA = 1.875;
P_R = 1;
P_RC = 1.5;



% b) M_opt
max_logos_max = 1;
logos_max = 0.65;                % (v / vmax)_max
v_max_mine = v_max * logos_max;


rho_air = 1.225;                % kgr/m^3
v_a = 2.5;                      % 2.5 m/s
C_d = 0.005;
g = 9.81;                       % m/s^2
A_UAV = desired_side^2;
D_w = (rho_air * v_a^2 * C_d * A_UAV) / (2*g);


% a) M_max_relay
M_max_relay = find_M_max_relay(logos_max, T_max, B_w, D_w, UAV_w, Ant_w);
M_max_relay = floor(M_max_relay);




% Max antennas for the UAV to move with a standard logos_max (M_max)
fprintf('0 <= a <= a_max\n<strong>0 <= a <= %.2f\n</strong>', max_logos_max);
fprintf('<strong>Select: a = %.2f\n\n</strong>', logos_max);

display('**********************************************************');
fprintf('<strong>Max relay antennas ---> T_max [kgr]</strong>\n');
display(' ');


[Ltmin_minutes_relay, P_tot_max_relay, W_relay] = find_Lt_minutes_relay(M_max_relay, logos_max, lambda, UAV_w, Ant_w, T_max, rho_air, v_a, C_d, g, v_max, B_w, B_c, desired_side, a_PA, P_R, P_RC);
Ltmin_sec_relay = floor(Ltmin_minutes_relay * 60);
N_min_relay = floor(Ltmin_sec_relay / delta_t);
pretty_Lt(M_max_relay, logos_max, P_tot_max_relay, Ltmin_minutes_relay, N_min_relay, delta_t, W_relay, T_max);
display('**********************************************************');
display(' ');







% M selection
[M_relay, sep] = select_M_relay(desired_side, lambda);
M_rec = M_relay / 2;
M_tr = M_relay / 2;

display('**********************************************************');
fprintf('For a %dcm x %dcm frame, we need:\n', 100 * desired_side, 100 * desired_side);
fprintf('<strong>M = %d antennas with separation = %.2fcm</strong>\n', M_relay, 100*sep);
fprintf('located across the diagonal of frame\n');
display('**********************************************************');
fprintf('\n\n\n');

m_list = 1 : M_relay;
D_w = (rho_air * v_a^2 * C_d * A_UAV) / (2*g);

[Lt_minutes_relay, P_tot_relay, W_relay] = find_Lt_minutes_relay(M_relay, logos_max, lambda, UAV_w, Ant_w, T_max, rho_air, v_a, C_d, g, v_max, B_w, B_c, desired_side, a_PA, P_R, P_RC);
Lt_sec_relay = floor(Lt_minutes_relay * 60);
N_relay = floor(Lt_sec_relay / delta_t);



% CONSTRAINTS FOR 'N'
% 1. N must be even (for the spiral to be made)
% 2. N must be in form: N = K * (sth), where (sth) will be
% the time duration of benchmarks communication

display('**********************************************************');
fprintf('For M=%d and a=%.2f, we have:\n\n', M_relay, logos_max);
fprintf('N = %d, but:\n', N_relay);
for N_new = N_relay : -1 : 0
    if mod(N_new, 2) == 0 && mod(N_new, K) == 0
        N_relay = N_new;
        fprintf('1. N must be dividable by 2 (spiral)\n');
        fprintf('2. N must be dividable by %d (benchmarks)\n', K);
        fprintf('<strong>N = %d\n</strong>', N_relay);
        break;
    end
end
display('**********************************************************');
fprintf('\n\n\n');

display('**********************************************************');
fprintf('<strong>I select this value of M:</strong>\n\n');
pretty_Lt(M_relay, logos_max, P_tot_relay, Lt_minutes_relay, N_relay, delta_t, W_relay, T_max);
display('**********************************************************');
display(' ');








%% 4. Technical RIS Stats - Weights and Dimensions
fprintf('\n\n\n\n');
display('**********************************************************');
a = logos_max;
W_total = REL_info(a, B_w, D_w, T_max, M_relay, Ant_w, UAV_w, f, B, lambda, N_relay, sep);
display('**********************************************************');
display(' ');








%% 5. Initialize A, P, Q, Theta
% *************************************************************************
% *************************************************************************
% *************************************************************************
% *************************************************************************

n_list = 1 : N_relay;

% -------- VARIABLES ---------
% TDMA  = 1 x N ----> A = K x N
% POWER = 1 x N ----> P = K x N
% Q = 2 x N
% Theta = M x N, total_gain_optimized = 1 x N
% SNR = 1 x N, SNR_av_optimized = scalar
% ----------------------------

strofes = 3;
alpha = 10;
beta = (LEN_X/2 - alpha) / (strofes*2*pi);

iter = 0;
Q = traj_ellipsis(q_cp, LEN_X, LEN_Y, N_relay);
Q = traj_romvos(q_cp, LEN_X, LEN_Y, N_relay);
Q = traj_speira2(alpha, beta, N_relay, strofes, LEN_X, q_cp);
Q = traj_traveling_salesman(q_cp, N_relay, aggregators, K);
Q = traj_romvos(q_cp, LEN_X, LEN_Y, N_relay);
% Q = traj_straight_line(q_cp, LEN_X, LEN_Y, N);
[TDMA, schedule, timestamps] = TDMA_fair_benchmark(Q, K, N_relay, aggregators);
A = TDMA_to_A(TDMA, K, N_relay);
POWER = P_tr * ones(1, N_relay);
P = POWER_to_P(POWER, K, N_relay, TDMA);








%% 6. Plots, velocity, acceleration
min_sumrates_Mbps_list_relay = [];
sum_sumrates_Mbps_list_relay = [];
velocity = calculate_velocity(Q, delta_t, N_relay);
plot_trajectory_TDMA(Q, LEN_X, LEN_Y, aggregators, ...
    x_cp, y_cp, xb, yb, K, colors, iter, schedule, timestamps);
plot(Q(1, :), Q(2, :));




%% 7. Iteration 0 - Statistics
% *************************************************************************
% *************************************************************************
% *************************************************************************
% *************************************************************************

P_rel = M_relay*P_RC + (a_PA+1)*P_R;

display(' ');
display(' ');
display(' ');
display(' ');
display('**********************************************************');
display('**********************************************************');
display(' ');
display('MAIN FUNCTION PARAMETERS');
display(' ');
fprintf('<strong>a = %.2f\nM = %d antennas\nN = %d timeslots\nK = %d sensors\n</strong>\n\n\n', logos_max, M_relay, N_relay, K);


fprintf('<strong>-----------------------------</strong>\n');
fprintf('<strong>-------- Iteration %d --------</strong>\n', iter);
fprintf('<strong>-----------------------------</strong>\n\n');

sumrates_relay = find_sumrates_relay(Q, TDMA, K, N_relay, B, aggregators, qb, H, h_b, G, P_tr, P_noise, C0_1m, M_rec, M_tr, P_ant, kappa);
[min_sumrate_relay, min_agg_relay] = min(sumrates_relay);
min_sumrates_Mbps_list_relay = [min_sumrates_Mbps_list_relay, min_sumrate_relay/1e6];
SUM_SUMRATES_MBPS_RELAY = sum(sumrates_relay) / 1e6;
sum_sumrates_Mbps_list_relay = [sum_sumrates_Mbps_list_relay, SUM_SUMRATES_MBPS_RELAY];


fprintf('Schedule of Fair TDMA: \n');
SCHEDULE = "";
for k = 1 : K
    if k < K
        SCHEDULE = SCHEDULE + num2str(schedule(k)) + " ---> ";
    else
        SCHEDULE = SCHEDULE + num2str(schedule(k));
    end
end
fprintf('<strong>%s</strong>\n\n\n', SCHEDULE);


display('******** QoS ********');
sumrates_Mbps = sumrates_relay / 1e6;
fprintf('Sumrates [Mbps] = \n[');
fprintf('%g, ', sumrates_Mbps(1:end-1));
fprintf('%g]\n\n', sumrates_Mbps(end));
fprintf('<strong>Min (sumrate) = %.2f Mbps\nfor agg = %d\n</strong>', min_sumrate_relay/1e6, min_agg_relay);
display(' ');








%% 8. Alternate Algorithm - Iterations - Improve DF
% *************************************************************************
% *************************************************************************
% *************************************************************************
% *************************************************************************

max_iter = 1;   % In last iteration, I optimize only 'A', not 'Q'
improve_perc = Inf;
improve_perc_thr = 3/100;




while iter < max_iter


    %% 8a. Basics

    iter = iter + 1;
    initial_A = A;
    initial_Q = Q;
    initial_Theta = Theta;
    initial_total_gain_optimized = total_gain_optimized;
    initial_SNR = SNR;
    fprintf('<strong>-----------------------------</strong>\n');
    fprintf('<strong>-------- Iteration %d --------</strong>\n', iter);
    fprintf('<strong>-----------------------------</strong>\n\n');

    global F C1 C2 B_c;
    F = (P_tr/P_noise) * Gt*Gr * C0_1m^2 * M_relay^2;
    C1 = UAV_w + B_w + M_relay*Ant_w;
    C2 = (T_max - C1) / v_max;




    %% 8b. Optimization of 'A' - Find optimal TDMA
    fprintf('<strong>Subproblem 1 - Optimize A\n</strong>');

    % Create the d1, d2 list using the OLD TRAJECTORY Q
    d1_list = zeros(N_relay, K);
    d2_list = zeros(N_relay, K);
    for n = 1 : N_relay
        UAV = initial_Q(:, n);
        if iter>1
            UAV = [qq(n),qq(N_relay+n)];
        end
        for k = 1 : K
            agg = aggregators(:, k);
            d1_list(n, k) = sqrt( (UAV(1)-agg(1))^2 + (UAV(2)-agg(2))^2 + H^2 );
            d2_list(n, k) = sqrt( (UAV(1)-qb(1))^2 + (UAV(2)-qb(2))^2 + (H-h_b)^2 );
        end
    end

    % Solution VECTOR with: length = NK + 1
    % z = [a_1[1], ..., a_K[1],
    %      a_1[2], ..., a_K[2],
    %      ...................,
    %      a_1[N], ..., a_K[N],
    %      t]



    % 1a. Objective function 'f'
    f = zeros(N_relay*K+1, 1)';
    f(N_relay*K+1) = -1;                  % min (-t)
    % f = @(z) (-z(end));

    % 1b. Binary variables a_k[n] -------> (C3)
    intcon = 1 : N_relay*K;
    lb = zeros(N_relay*K+1, 1);
    ub = ones(N_relay*K+1, 1);
    lb(N_relay*K+1) = -Inf;
    ub(N_relay*K+1) = Inf;

    % 1c. Equality Constraint -----------> (C4)
    % For every n:   sum a_k[n] = 1
    Aeq = zeros(N_relay, N_relay*K+1);
    for n = 1 : N_relay
        % indeces = (n-1)*K + 1:K;
        Aeq(n, (n-1)*K+1:(n-1)*K+K) = 1;
    end
    beq = ones(N_relay, 1);



    % 1d. Inequality Constraints --------> (C7)
    % For every k: -sum(....) + t <= 0
    Aineq = zeros(K, N_relay*K+1);
    for k = 1 : K
        for n = 1 : N_relay
            d1 = d1_list(n, k);
            d2 = d2_list(n, k);
            Aineq(k, k+(n-1)*K) = -log2(1 + F/d1^2/d2^2);
        end
    end
    Aineq(:, N_relay*K+1) = ones(K, 1);           % for 't'
    bineq = zeros(K, 1);



    fprintf('* (C4) has N=%d equality constraints:\n', N_relay);
    fprintf('<strong> Aeq  = %dx%d</strong>, ', size(Aeq));
    fprintf('<strong> beq  = %dx%d</strong>\n', size(beq));
    fprintf('* (C7) has K=%d inequality constraints:\n', K');
    fprintf('<strong>Aineq = %dx%d</strong>, ', size(Aineq));
    fprintf('<strong>bineq = %dx%d</strong>\n\n', size(bineq));





    % 1e. Options (MaxTime = 60 seconds)
    options = optimoptions(@intlinprog, 'MaxTime', 60);
    [z, opt_val, exitflag, output] = ...
    intlinprog(f, intcon, Aineq, bineq, Aeq, beq, lb, ub,[],options);
    % This optimal value is in bps/Hz, so I have to multiply by
    % B = 1 MHz to find the value in bps or Mbps
    t_optimal_A = - opt_val * B;
    t_optimal_A_Mbps = t_optimal_A / 1e6;




    % 1f. Check if 'z' vector is alright
    SUM_Z = sum(z) - z(end)

    % 1g. Transform the vector 'z' into the 'A' matrix or 'TDMA' vector
    TDMA_new = zeros(1, N_relay);
    for n = 1 : N_relay
        subvector = z( (n-1)*K+1 : (n-1)*K+K );
        agg_comm = find(abs(subvector-1) <= 1e-4);
        TDMA_new(n) = agg_comm;
    end
    A_new = TDMA_to_A(TDMA_new, K, N_relay);
    global A_new;


    % 1h. Print statistics and improvement
    [Theta_A, total_gain_A, SNR_A, SNR_av_dB_A] = ...
    optimize_phase_shifts(N_relay, M_relay, TDMA_new, P_tr, initial_Q, aggregators, H, h_b, qb, ...
    lambda, d, C0_1m, P_noise);

    sumrates_A = find_sumrates(initial_Q, TDMA_new, K, N_relay, B, aggregators, qb, H, h_b, G, P_tr, P_noise, C0_1m, M_relay);
    [min_sumrate_A, min_agg_A] = min(sumrates_A);
    sumrates_A_Mbps = sumrates_A / 1e6;
    min_sumrates_Mbps_list_relay = [min_sumrates_Mbps_list_relay, min_sumrate_A/1e6];

    SUM_SUMRATES_MBPS_A = sum(sumrates_A_Mbps);
    sum_sumrates_Mbps_list_relay = [sum_sumrates_Mbps_list_relay, SUM_SUMRATES_MBPS_A];

    display('******** QoS ********');
    disp("SNR_av_dB = " + SNR_av_dB_A + " dB");
    fprintf('Sumrates [Mbps] = \n[');
    fprintf('%g, ', sumrates_A_Mbps(1:end-1));
    fprintf('%g]\n\n', sumrates_A_Mbps(end));
    fprintf('<strong>Optimization  = %.2f Mbps</strong>\n', t_optimal_A_Mbps);
    fprintf('<strong>Min (sumrate) = %.2f Mbps\nfor agg = %d\n</strong>', min_sumrate_A/1e6, min_agg_A);
    fprintf('\n\n');

    substep = " - Optimization TDMA - A";

    plot_trajectory_TDMA2(initial_Q, LEN_X, LEN_Y, aggregators, ...
    x_cp, y_cp, xb, yb, K, colors, iter, substep, TDMA_new, N_relay);
    plot(initial_Q(1, :), initial_Q(2, :));
    OK = 1000;









    %% 8c. Optimization of 'Q' - Find optimal trajectory

    if iter <= max_iter-1


        fprintf('<strong>\nSubproblem 2 - Optimize Q\n</strong>');

        % Solution vector in this form:
        % qq = [.... xn .... | .... yn .... |
        %       .... skn ... | .... tkn ... | t]
        % with:
        % skn = [s1(1), ..., sK(1),  s1(2), ..., sK(2), ..., s1(N), ..., sK(N)]
        % tkn = [t1(1), ..., tK(1),  t1(2), ..., tK(2), ..., t1(N), ..., tK(N)]



        % Objective function and options

       exitflag =0;
       if iter ==1
           qq0 = [1*initial_Q(1, :) 1*initial_Q(2, :) 30*rand(1, 2*K*N_relay) t_optimal_A_Mbps];
       else
           qq0 = [1*initial_Q(1, :) 1*initial_Q(2, :) qq(2*N_relay+1:end)];
       end
       iterr = 0;
       global init
       obj = @(qq)(-1*qq(end));
    %    exitflag==0 &&
       while (iterr<2)
        init = qq0/1.1;
        iterr= iterr +1;
        lb = -LEN_X/2 * ones(2*N_relay+2*K*N_relay+1, 1);
        ub = LEN_X/2 * ones(2*N_relay+2*K*N_relay+1, 1);
        lb(2*N_relay+1 : 2*N_relay+2*K*N_relay) = 0;
        ub(2*N_relay+1 : 2*N_relay+2*K*N_relay) = 1000;
        lb(2*N_relay+2*K*N_relay+1) = 0;
        ub(2*N_relay+2*K*N_relay+1) = 1000;
        options = optimoptions(@fmincon,'Algorithm','interior-point',...
                      'MaxIter',1000,'MaxFunEvals',3000, 'Display', 'off', 'ConstraintTolerance', 0.01);
        if iterr==1
            [qq, fval, exitflag] = fmincon(@(qq)(0), qq0, [], [], [], [], lb, [], @nonlcon2_nikos, options);
            qq0=qq;
            [qq, fval, exitflag] = fmincon(obj, qq0, [], [], [], [], lb, [], @nonlcon2_nikos, options);
            exitflag
        end
        if iterr>1
           [qq, fval, exitflag] = fmincon(obj, qq0, [], [], [], [], lb, [], @nonlcon2_nikos, options);
           exitflag
        end
         if exitflag ==0 || exitflag==1
            qq0=qq;
        end
       end

        % Transform vector 'qq' into new trajectory Q_new
        Q_new = [qq(1:N_relay); qq(N_relay+1:2*N_relay)];
        Q=Q_new;
        t_optimal_Q = fval;
        t_optimal_Q_Mbps = t_optimal_Q / 1e6;


        % Print statistics and improvement
        [Theta_Q, total_gain_Q, SNR_Q, SNR_av_dB_Q] = ...
        optimize_phase_shifts(N_relay, M_relay, TDMA_new, P_tr, Q_new, aggregators, H, h_b, qb, ...
        lambda, d, C0_1m, P_noise);

        sumrates_Q = find_sumrates(Q_new, TDMA_new, K, N_relay, B, aggregators, qb, H, h_b, G, P_tr, P_noise, C0_1m, M_relay);
        [min_sumrate_Q, min_agg_Q] = min(sumrates_Q);
        sumrates_Q_Mbps = sumrates_Q / 1e6;
        min_sumrates_Mbps_list_relay = [min_sumrates_Mbps_list_relay, min_sumrate_Q/1e6];

        SUM_SUMRATES_MBPS_Q = sum(sumrates_Q_Mbps);
        sum_sumrates_Mbps_list_relay = [sum_sumrates_Mbps_list_relay, SUM_SUMRATES_MBPS_Q];

        display('******** QoS ********');
        disp("SNR_av_dB = " + SNR_av_dB_Q + " dB");
        fprintf('Sumrates [Mbps] = \n[');
        fprintf('%g, ', sumrates_Q_Mbps(1:end-1));
        fprintf('%g]\n\n', sumrates_Q_Mbps(end));
        fprintf('<strong>Optimization  = %.2f Mbps</strong>\n', t_optimal_Q_Mbps);
        fprintf('<strong>Min (sumrate) = %.2f Mbps\nfor agg = %d\n</strong>', min_sumrate_Q/1e6, min_agg_Q);
        fprintf('\n\n');

        substep = " - Optimization TDMA - Q";

        plot_trajectory_TDMA2(Q_new, LEN_X, LEN_Y, aggregators, ...
        x_cp, y_cp, xb, yb, K, colors, iter, substep, TDMA_new, N_relay);
        plot(Q_new(1, :), Q_new(2, :));
        OK = 1000;

    end

end








%% 9. Plots of sumrates (min and sum)
min_sumrates_Mbps_list_relay
sum_sumrates_Mbps_list_relay

names = ["Initialization"];
for iter = 1 : max_iter-1
    str1 = "Iter." + num2str(iter) + " - A";
    str2 = "Iter." + num2str(iter) + " - Q";
    names = [names, str1, str2];
end
str3 = "Iter." + num2str(max_iter) + " - A";
names = [names str3];
X = categorical(names);

figure();
bar(min_sumrates_Mbps_list_relay);
xticklabels(X);
ylabel('Min sumrate [Mbps]');
title('Min sumrate [Mbps] vs Steps');

figure();
bar(sum_sumrates_Mbps_list_relay);
xticklabels(X);
ylabel('Sum of sumrates [Mbps]');
title('Sum of sumrates [Mbps] vs Steps');

toc