%% ECE 602 Course Project: Throughput Maximization in Wireless Powered Communication Networks
% *Team member*: Huaqing Wu (20704063), Haixia Peng (20535999)

%% Problem formulation
% A “harvest-then-transmit” protocol is proposed for wireless powered communication networks (WPCNs), where all users first harvest the 
% wireless energy broadcast by the hybrid access point (H-AP) in the downlink (DL) and then send their independent information to the 
% H-AP in the uplink (UL) by time-division-multiple-access (TDMA). Authors were interested in maximizing the UL throughput of by optimally
% allocating the time for the DL wireless energy transfer (WET) by the H-AP and the UL wireless information transmissions (WITs) by 
% different users. With the proposed protocol, the authors first maximize the sum-throughput of the WPCN (referred to as *sum-throughput 
% maximization problem*) by jointly optimizing the time allocated to the DL WET and the UL WITs given a total time constraint, 
% based on the users' DL and UL channels as well as their average harvested energy amount.
% However, the sum-sumthroughput maximization solution is shown to allocate substantially more time to the near
% users than the far users, thus resulting in unfair achievable rates among different users, reffered to as “doubly near-far” phenomenon.
% To overcome this phenomenon, by proposing a new performance metric referred to as common-throughput with the additional
% constraint that all users should be allocated with an equal rate in their UL WITs regardless of their distances to the H-AP, a 
% *common-throughput maximization problem* is formulated. 
% 
% The sum-throughput maximization problem is formulated as
%
% $$(P1): \max_{\tau} R_{sum}(\tau)$
%
% $$s.t. \sum_{i=0}^K \tau_i \leq 1,$
%
% $$\tau_i \geq 0, i = 0,1,\ldots, K$
%
% where, $$K$ is the number of users, $$\tau_i$ represent the time portions in each block allocated to the H-AP ($$i = 0$) and users 
% $$U_1, \ldots, U_K$, $$R_{sum}(\tau)$ is the sum-throughput of all users and given by $$R_{sum}(\tau) = \sum_{i=1}^K R_i(\tau)$. 
% $$R_i(\tau)$ is the achievable UL throughput of user $$U_i$ in bits/second/Hz (bps/Hz), and can be expressed as 
% $$R_i(\tau) = \tau_i log_2 (1 + \frac{g_iP_i}{\Gamma\sigma^2}) = \tau_i
% log_2 (1 + \gamma_i\frac{\tau_0}{\tau_i})$, where $$P_i$ denotes the
% average transmit power at $$U_i$, $$g_i$ is the channel power gains of the
% reversed UL channel to user $$U_i$, $$\Gamma$ represents the signal-to-noise
% ratio (SNR) gap from the additive white Gaussian noise (AWGN) channel capacity due to a practical modulation and
% coding scheme (MCS) used, $$\gamma_i = \frac{\zeta h_i g_i
% P_A}{\Gamma\sigma^2}, i=1,\ldots,K$. In which, $$P_A$ denotes the transmit power at the
% H-AP, $$h_i$ is the channel power gains of the DL channel from the H-AP to
% user $$U_i$, and $$\zeta$ is the energy harvesting efficiency at each
% receiver (it is assumed to be same for each receiver).
%
% The common-throughput maximization problem is formulated as
%
% $$(P2): \max_{\overline{R},\tau} \overline{R}$
%
% $$s.t. R_i(\tau) \geq \overline{R},i = 1,\ldots, K$
%
% $$\tau \in \mathcal{D}$
%
% where, $$\overline{R}$ denotes the common-throughput and $$\mathcal{D}$ is
% the feasible set of $$\tau$ specified by $$\sum_{i=0}^K \tau_i \leq 1$ and $$\tau_i \geq 0$ $$(i = 0,1,\ldots, K)$.
%
% Figure 3 in the paper shows the throughput given in (P1) for the special case of one single user in the network, i.e., K = 1, versus
% the time allocated to the DL WET, $$\tau_0$, with $$\gamma_1 = 10dB$,
% assuming that $$\sum_{i=0}^K \tau_i \leq 1$ holds with equality, i.e., for the UL WIT $$\tau_1 =1-\tau_0$.

gamma1 = 10;
tau0 = 0 : 0.001 : 1;
R = zeros(1,length(tau0)); R_max = 0;
for i = 1:length(tau0)
    R(i) = (1-tau0(i))*log2(1 + gamma1*tau0(i)/(1-tau0(i)));
    if R(i) > R_max
        R_max = R(i);
        tau_opt = tau0(i);
    end
end
figure;
plot(tau0,R,'b-','linewidth',2)
grid on;
xlabel('\tau_0','fontsize',12);
ylabel('Throughput (bps/Hz)','fontsize',12);
title('Fig. 3 in the Paper');
disp(['the optimal value of tau_0 is:', num2str(tau_opt)])

%%
% Observing the figure, we can find that there exist an
% optimal value of $$\tau_0$, which means that the time allocated to the DL
% wireless energy transfer should not be too large or too small.

%% Solution for Problem P1
% For the sum-throughput maximization problem (P1), since $$R_i(\tau)$ is a
% concave function of $$\tau$ for any given $$i \in \{1, \ldots, K\}$, which
% follows that $$R_{sum}(\tau)$ is also a concave function of $$\tau$ since it is the summation of $$R_i(\tau)$'s. Therefore, (P1) is 
% a convex optimization problem, and thus can be solved by convex
% optimization techniques. On the other hand, lemma 3.2 is defined as: given $$A > 0$, there exists a unique $$z^\ast > 1$ that is the solution of 
% $$f(z) = A$, where $$f(z) = zlnz - z + 1$, $$z \geq 0$. *Fig.4* shows that f(z) is a convex function over $$z\geq 0$ where the minimum
% is attained at $$z = 1$ with $$f(1) = 0$. Therefore, given $$0 < A \leq 1$, there are two different solutions for $$f(z) = A$, among
% which one is smaller than $$1$ and the other is larger than $$1$, i.e., $$z^\ast > 1$. On the other hand, if $$A>1$, there is only one
% solution for $$f(z)=A$, which is larger than $$1$, i.e., $$z^\ast > 1$.
%
z = 0 : 0.001 : 5;
f = z.*log(z) - z +1;
figure;
plot(z,f,'b-','linewidth',2)
grid on;
xlabel('z','fontsize',12);
ylabel('f(z)','fontsize',12);
title('Fig. 4 in the Paper');
%%
% Thus, based on Proposition 3.1 in the paper, the optimal time allocation solution for (P1), denoted by
% $$\tau^\ast=[\tau_0^\ast \ \tau_1^\ast \ \ldots \ \tau_K^\ast]$, is given by
% 
% $$\tau_0^\ast = \frac{z^\ast -1}{A + z^\ast -1}$,
%
% $$\tau_i^\ast = \frac{\gamma_i}{A + z^\ast -1}, i = 1,\ldots, K$
%
% where $$A = \sum_{i=i}^K \gamma_i > 1$ and $$z^\ast > 1$ is the
% corresponding solution of f (z) = A as given by Lemma 3.2.  
% 
% The optimal time allocation solution for (P1) follows that $$exists$ $$\tau_i^\ast \propto D_i^{-(\alpha_d + \alpha_u)},i = 1,\ldots, K$, which results in an unfair time and throughput
% allocation among users in the WPCN, i.e., the "doubly near-far" phenomenon. Here, $$\alpha_d \geq 2$ and $$\alpha_u \geq 2$
% denote the channel pathloss exponents in the DL and UL,respectively, and $$D_i$ denotes the distance between the H-AP
% and $$U_i$. *Fig.5* is used to further illustrate the "doubly near-far" phenomenon.
K = 2;
gamma_dB = [22, 10]; % the values of gamma in dB
gamma = 10.^(gamma_dB ./ 10);

% plot Fig 5 in the paper
tau_temp = 0.001:0.001:1;
[tau_1,tau_2] = meshgrid(0.001:0.001:1);
N_5 = length(tau_1);
R1_5= zeros(size(tau_1));
R2_5 = zeros(size(tau_1));
R_sum_5 = zeros(size(tau_1));
for i = 1 : N_5
    for j = 1 : N_5
        if tau_1(i,j) + tau_2(i,j) < 1
            R1_5(i,j) = tau_1(i,j)*log2(1+gamma(1)*(1-tau_1(i,j)-tau_2(i,j))/tau_1(i,j));
            R2_5(i,j) = tau_2(i,j)*log2(1+gamma(2)*(1-tau_1(i,j)-tau_2(i,j))/tau_2(i,j));
            R_sum_5(i,j) = R1_5(i,j) + R2_5(i,j);
        end
    end
end
figure;
contour(tau_1,tau_2,R_sum_5,10);
hold on
grid on;
xlabel('\tau_1 (allocated to U_1)','fontsize',12);
ylabel('\tau_2 (allocated to U_2)','fontsize',12);
title('Fig. 5 in the Paper');

% find the maximum value of the throughput and the corresponding tau
[m,n] = find(R_sum_5 == max(max(R_sum_5)));
tau_opt = [tau_1(m,n), tau_2(m,n)];
R1_opt_5 = R1_5(m,n);
R2_opt_5 = R2_5(m,n);
plot(tau_opt(1), tau_opt(2),'*')
disp(['the optimal value of tau is: [tau_0, tau_1, tau_2] = ', num2str([1-sum(tau_opt), tau_opt])])

% when calculating tau by proposition 3.1
tau_opt_5 = Optimal_tau_func_pro3(K, gamma);
disp(['the optimal value of tau based on Proposition 3.1 is: [tau_0, tau_1, tau_2] = ', num2str([tau_opt_5])])
%%
% Fig.5 shows the sum-throughput $$R_{sum}(\tau)$ versus UL WIT time allocation 
% $$\tau_1$ and $$\tau_2$ for a two-user network with $$K = 2$ and $$D_1 = \frac{1}{2} D_2$. It is assumed that the channel reciprocity 
% holds for the DL and UL and thus $$h_i = g_i, i = 1, 2$, with $$\alpha_d=
% \alpha_u= 2$, and $$\gamma_1 = 22dB$, $$\gamma_2 = 10dB$, with
% $$\gamma_1/\gamma_2 = (D_2/D_1)^{\alpha_d + \alpha_u}$. We can also see
% from the results that, the optimal solution found by brute force is the
% same as that obtained by using Proposition 3.1, which validates the
% correctness of the proposition.
% 
% With the same setup as for Fig. 5, Fig. 6 shows the optimal throughput of $$U_2$, $$R_2(\tau^\ast)$, normalized by
% that of $$U_1$, $$R_1(\tau^\ast)$, in a WPCN versus that in a conventional TDMA network for different values of the pathloss exponent
% $$\alpha$, with $$\alpha_d= \alpha_u= \alpha$. For the WPCN, $$\gamma_1$ is set to be fixed as $$\gamma_1 = 22dB$ for $$U_1$, while 
% $$\gamma_2 = 10,7,4,1$,and $$-2dB$ for $$\alpha= 2, 2.5, 3, 3.5$, and $$4$, respectively, since $$\gamma_1/\gamma_2 = (D_2/D_1)^{2\alpha}$.
% For the conventional TDMA network, $$\overline{E}$ for both $$U_1$ and
% $$U_2$ are assumed to be $$\overline{E} = \frac{1}{2}(E_1 + E_2)$ with
% $$E_1$ and $$E_2$ denoting the average harvested energy at $$U_1$
% and $$U_2$ in the WPCN under comparison, respectively; it then follows that for the conventional TDMA network, $$\gamma_1 = 13dB$ 
% for $$U_1$, and $$\gamma_2 = 7.1, 5.5, 4.0, 2.4$,and $$0.9dB$ for $$\alpha= 2, 2.5, 3, 3.5$, and $$4$, respectively.
alpha = 2 : 0.5 : 4;
N = length(alpha);

% For WPCN, the values of gamma are set as follows according to the paper
gamma1_wpcn = 10^(2.2);
gamma2_wpcn = [10^(1), 10^(0.7), 10^(0.4), 10^(0.1), 10^(-0.2)];
R2_R1_wpcn = zeros(1, N);
for i = 1 : N
    gamma_w = [gamma1_wpcn,gamma2_wpcn(i)];
    tau_opt_w = Optimal_tau_func_pro3(K, gamma_w);
    R1 =  tau_opt_w(2) * log2(1 + gamma_w(1)*tau_opt_w(1) / tau_opt_w(2));
    R2 =  tau_opt_w(3) * log2(1 + gamma_w(2)*tau_opt_w(1) / tau_opt_w(3));
    R2_R1_wpcn(i) = R2 / R1;
end

% For conventional TDMA network
gamma1_t = 10^(1.3);
gamma2_t = [10^(0.71) 10^(0.55) 10^(0.4) 10^(0.24) 10^(0.09)];
R2_R1_tdma = zeros(1, N);
for i = 1 : N
    gamma_t = [gamma1_t,gamma2_t(i)];
    tau_opt_t = Optimal_tau_func_pro3(K, gamma_t);
    R1 =  tau_opt_t(2) * log2(1 + gamma_t(1)*tau_opt_t(1) / tau_opt_t(2));
    R2 =  tau_opt_t(3) * log2(1 + gamma_t(2)*tau_opt_t(1) / tau_opt_t(3));
    R2_R1_tdma(i) = R2 / R1;
end
% plot Fig 6 in the paper
figure
semilogy(alpha,R2_R1_tdma, 'r-o ', alpha,R2_R1_wpcn, 'b-d ', 'linewidth',2)
xlabel('Pathloss Exponent, \alpha_d = \alpha_u = \alpha','fontsize',12);
ylabel('R_2 (\tau^*) / R_1 (\tau^*)','fontsize',12);
legend('Conventional TDMA network','WPCN');
grid on
title('Fig. 6 in the Paper');

%% Solution for Problem P2
% To solve (P2), given any $$\overline{R}>0$, the following feasiblity problem is first considered,
%
% $$Find \  \tau$
%
% $$s.t. R_i(\tau) \geq \overline{R},i = 1,\ldots, K$
%
% $$\tau \in \mathcal{D}    \,\,\,\,\,\,\,\,\,\,\,\,\,\,(12)$
%
% Since the problem in (12) is convex, the authors consider its Lagrangian given by
%
% $$\mathcal{L}(\tau,\lambda)=-\sum_{i=1}^K \lambda_i (R_i(\tau)-\overline{R}) \,\,\,\,\,\,\,\,\,\,\,\,\,\,(13)$
%
% where $$\lambda = [\lambda_1,\ldots, \lambda_K] \succeq 0$ consists of the Lagrange multipliers associated
% with the K user throughput constraints in problem (12). The dual function of problem (12) is then given by
% 
% $$\mathcal{G}(\lambda)=\min_{\tau \in \mathcal{D}} \mathcal{L}(\tau,\lambda) \,\,\,\,\,\,\,\,\,\,\,\,\,\,(14)$
%
% The dual function $$\mathcal{G}(\lambda)$ can be used to determine whether
% problem (12) is feasible, as provided in lemma 4.1: For a given $$\overline{R} > 0$, problem (12) is infeasible
% if and only if there exists an $$\lambda \geq 0$ such that $$\mathcal{G}(\lambda)> 0$.
%
% Thus, $$\mathcal{G}(\lambda)$ in (14) can be obtained for a given $$\lambda \geq 0$ by
% solving the following weighted sum-throughput maximization problem, which follows from (13).
%
% $$\max_{\tau} \sum_{i=1}^K \lambda_i R_i(\tau)$
%
% $$s.t. \tau \in \mathcal{D}.    \,\,\,\,\,\,\,\,\,\,\,\,\,\,(15)$
%
% The problem in (15) is convex and thus can be solved by convex optimization techniques. Therefore, for the sum-throughput 
% maximization case with $$\lambda_i = 1, \forall i$, the optimal time allocation solution for the weighted sum-throughput 
% maximization problem in (15) is given as following proposition.
% 
% Proposition 4.1: Given $$\lambda \geq 0$, the optimal time allocation solution for (15), denoted by 
% $$\tau^\ast=[\tau_0^\ast \tau_1^\ast \ldots \tau_K^\ast]$, is 
% 
% $$\tau_0^\ast = \frac{1}{1+ \sum_{j=1}^K (\gamma_j/z_j^\ast)}$,
%
% $$\tau_i^\ast = \frac{\gamma_i/z_i^\ast}{1+ \sum_{j=1}^K (\gamma_j/z_j^\ast)}, i = 1,\ldots, K$
% 
% where $$z_i^\ast$, $$i = 1,\ldots, K$, is the solution of the following
% equations: 
% 
% $$ln(1+z_i) - \frac{z_i}{1+z_i} = \frac{\mu^\ast}{\lambda_i}ln2$,
%
% $$\sum_{i=1}^K \frac{\lambda_i \gamma_i}{1+z_i} = \mu^\ast ln2$,
%
% with $$\mu^\ast > 0$ being a constant. 
% 
% Based on the above analysis, one *iterative optimization algorithm* to solve (P2) is given as
% follows,
% 
% 1) *Initialize* $$R_{min} = 0$, $$R_{max} > \overline{R}^\ast.$
%
% 2) *Repeat*
% 
% # $$\overline{R} = \frac{1}{2}(R_{min} + R_{max})$.
% # Initialize $$\lambda \geq 0$.
% # Given $$\lambda$, solve the problem in (15) by Proposition 4.1.
% # Compute $$\mathcal{G}(\lambda)$ using (13).
% # If $$\mathcal{G}(\lambda) > 0$, $$\overline{R}$ is infeasible, set $$R_{max} \longleftarrow \overline{R}$, go to step 1. 
% Otherwise, update $$\lambda$ using sub-gradient based algorithms (e.g., the ellipsoid method) and the subgradient
% of $$\mathcal{G}(\lambda)$ given by (20). If the stopping criteria of the sub-gradient based algorithm is not met, go to step 3.
% # Set $$R_{min} \longleftarrow \overline{R}$.
%
% 3) *Until* $$R_{max}-R_{min} < \delta$, where $$\delta > 0$ is a given error tolerance.
% 
% $$\nu_i=\tau_i^\ast log_2 (1 + \frac{\gamma_i \tau_0^\ast}{\tau_i^\ast})-\overline{R}, 1\geq i \leq K, \,\,\,\,\,\,\,\,\, (20)$
% 
% where $$\nu ={[\nu_1 \nu_2 \ldots \nu_K]}^T$ is the subgradient of $$\mathcal{G}(\lambda)$.  
%
% Fig. 7 shows the common-throughput in bps/Hz versus $$\tau_1$ and $$\tau_2$ for the same two-user channel setup as for Fig. 5.
[tau_1,tau_2] = meshgrid(0.001:0.001:1);
R_aver = zeros(size(tau_1));
for i = 1:length(tau_1)
    for j = 1:length(tau_1)
        if tau_1(i,j) + tau_2(i,j) < 1
            R1 = tau_1(i,j) * log2(1 + gamma(1)*(1-tau_1(i,j)-tau_2(i,j)) / tau_1(i,j));
            R2 = tau_2(i,j) * log2(1 + gamma(2)*(1-tau_1(i,j)-tau_2(i,j)) / tau_2(i,j));
            R_aver(i,j) = min(R1, R2);
        end
    end
end
figure
contour(tau_1,tau_2,R_aver,10);
hold on
grid on;
xlabel('\tau_1 (allocated to U_1)','fontsize',12);
ylabel('\tau_2 (allocated to U_2)','fontsize',12);
title('Fig. 7 in the Paper')

% find the optimal value of tau
[m,n] = find(R_aver == max(max(R_aver)));
tau_opt = [tau_1(m,n), tau_2(m,n)];
plot(tau_opt(1), tau_opt(2),'*')

% calculate the corresponding R(tau)
R1_opt_7 = tau_opt(1) * log2(1 + gamma(1)*(1-sum(tau_opt)) / tau_opt(1));
R2_opt_7 = tau_opt(2) * log2(1 + gamma(2)*(1-sum(tau_opt)) / tau_opt(2));
disp(['the optimal value of tau is: [tau_0, tau_1, tau_2] = ', num2str([1-sum(tau_opt), tau_opt])])
disp(['the optimal values of R1 and R2 are: [R1, R2] = ', num2str([R1_opt_7, R2_opt_7])])
disp(['the optimal common throughput is: min(R1, R2) = ', num2str(min(R1_opt_7, R2_opt_7))])
% calculate the throughput based on the optimal tau given in the paper
tau_given = [0.3683, 0.1386, 0.4932];
R1 = tau_given(2) * log2(1 + gamma(1)*tau_given(1) / tau_given(2));
R2 = tau_given(3) * log2(1 + gamma(2)*tau_given(1) / tau_given(3));
disp(['Based on the optimal tau obtained in the paper, [tau_0, tau_1, tau_2] = ', num2str([tau_given])])
disp(['the optimal values of R1 and R2 are: [R1, R2] = ', num2str([R1, R2])])
disp(['the optimal common throughput is: min(R1, R2) = ', num2str(min(R1, R2))])

%%
% From the comparing results, the value of $\tau^*$ obtained in the paper is not actually the optimal value.
% 
% We can also use the iterative algorithm given above to find the optimal
% value of $$\tau$. 

[tau_opt, R_tau, lambda] = optimal_tau_cal(K, gamma);

disp(['Based on the algorithm, the optimal tau is [tau_0, tau_1, tau_2] = ', num2str([tau_opt])])
disp(['the optimal values of R1 and R2 are: [R1, R2] = ', num2str([R_tau])])
disp(['the optimal common throughput is: min(R1, R2) = ', num2str(min(R_tau))])

%%
% The problem is, although using the iterative algorithm can obtain an
% acceptable solution, it takes a long time to converge. Actually, it takes
% longer time than using the brute force algorithm and we can see from the
% result that the optimal throughput achieved by using the iterative
% algorithm is worse than that obtained by using brute force.
% However, this obtained throughput is still larger than that obtained by
% using the $$\tau$ given in the paper. In other words, the algorithm
% implemented in this project converges better than the one used in the
% paper.
%
%
% Fig. 8 shows the ratio of the optimal time allocated to the far user
% $$U_2$ over that to the near user $$U_1$, i.e., $$\tau_2^\ast /\tau_1^\ast$ in (P1)
% versus (P2), with different values of the common pathloss exponent $$\alpha$ in both DL and UL, where the same two-user
% channel setup as for Fig. 5 is considered. It is assumed that
% $$\gamma_1$ for the near user $$U_1$ is fixed as $$\gamma_1 = 22dB$ and $$\gamma_2$ for the far user $$U_2$ is set the same as for Fig. 6.
% For P1, we can find the optimal $$\tau$ based on Proposition 3.1
R_p1 = zeros(1, N);
for i = 1 : N
    tau_opt = Optimal_tau_func_pro3(K, [gamma1_wpcn, gamma2_wpcn(i)]);
    R_p1(i) = tau_opt(3)/tau_opt(2); % the ratio of tau2 / tau1
end

% For P2, we can use the iterative optimization algorithm to optimize tau
R_p2 = zeros(1, N);
for i = 1 : N
    gamma = [gamma1_wpcn,gamma2_wpcn(i)];
    [tau_opt, R_tau, lambda] = optimal_tau_cal(K, gamma);
    R_p2(i) = tau_opt(3)/tau_opt(2);
end

figure
semilogy(alpha,R_p1, 'r-o ', alpha,R_p2, 'b-d ', 'linewidth',2)
xlabel('Pathloss Exponent(\alpha)','fontsize',12);
ylabel('\tau^*_2 / \tau^*_1','fontsize',12);
legend('(P1)','(P2)');
grid on
title('Fig. 8 in the Paper (Iterative Algorithm)')
%%
% However, the figure is not exactly the same with Fig. 8 in the paper. A
% possible reason is that the convergence performance of the algorithm here may
% be different from that implemented in the original paper. We tried  
% to find the optimal value of $$\tau$ by using brute force algorithm, and the
% obtained figure is almost identical with Fig. 8 in the paper.

% For P2, using the brute force algorithm to find the optimal tau
R_p2_bf = zeros(1, N);
for i = 1 : N
    gamma = [gamma1_wpcn,gamma2_wpcn(i)];
    [R_tau, tau_opt] =  Brute_force(K, gamma, 0.01, 2);
    R_p2_bf(i) = tau_opt(3)/tau_opt(2);
end

figure
semilogy(alpha,R_p1, 'r-o ', alpha,R_p2_bf, 'b-d ', 'linewidth',2)
xlabel('Pathloss Exponent(\alpha)','fontsize',12);
ylabel('\tau^*_2 / \tau^*_1','fontsize',12);
legend('(P1)','(P2)');
grid on
title('Fig. 8 in the Paper (Brute Force Algorithm)')

%%
% Fig. 9 shows the achievable throughput region of a two-user WPCN by solving
% the weighted sum-throughput maximization problem in (15) with different throughput weights for the near and far users,
% under the same channel setup as for Fig. 5. In the figure, "*" corresponds to the
% maximum sum-throughput and "o" corresponds to the maximum common-throughput. 
[tau_2,tau_1] = meshgrid(0.001:0.001:1);
N_9 = length(tau_1);
R_o1 = zeros(1, N_9);
R_o2 = zeros(1, N_9);
% plot the achievable throughput region of the two-user WPCN
for i = 1:N_9
    R_temp = R_sum_5(:, 1:i);
    [m,n] = find(R_temp == max(max(R_temp)));
    R_o1(i) = R1_5(m,n);
    R_o2(i) = R2_5(m,n);
end

figure
plot(R_o1, R_o2)
hold on
% plot the solution of P1 (obtained in Fig 5)
plot(R1_opt_5,R2_opt_5,'*');
hold on
% plot the solution of P2 (obtained in Fig 7)
plot(R1_opt_7,R2_opt_7,'o');
hold on
grid on;
xlabel('Throughput of U_1 (bps/Hz)','fontsize',12);
ylabel('Throughput of U_2 (bps/Hz)','fontsize',12);
title('Fig. 9 in the Paper');

%% Data Sources 
% In the simulation part, it is assumed that the channel reciprocity holds for the DL and
% UL and thus $$h_i = g_i, i = 1, \ldots, K$, with the same pathloss exponent $$\alpha_d= \alpha_u= \alpha$. Accordingly, 
% both the DL and UL channel power gains are modeled as $$h_i = g_i =
% 10^(-3)\rho_i^2 D_i^{-\alpha}, i = 1, \ldots, K$, where $$\rho_i$ represents the additional channel
% short-term fading which is assumed to be Rayleigh distributed, and $$\rho_i^2$ is an exponentially distributed random variable
% with unit mean. The values of parameters are shown in the following table,
%
% <html>
% <table border=1><tr><td>Parameter</td><td>Value</td></tr>
% <tr><td>Bandwidth</td><td>1 MHz</td></tr>
% <tr><td>Additional channel short-term fading </td><td>Rayleigh distribution</td></tr>
% <tr><td>Average signal power attenuation</td><td>30dB</td></tr>
% <tr><td>AWGN at the H-AP receiver</td><td>-160dBm/Hz</td></tr>
% <tr><td>Signal-to-noise ratio (SNR) gap</td><td>9.8dB</td></tr>
% <tr><td>Energy harvesting efficiency for WET</td><td>0.5</td></tr></table>
% </html>
% 
%
%% Visualization of Results 
% Fig. 10 shows the maximum sum-throughput versus the maximum common-throughput in the same WPCN with $$K = 2$, $$D_1 = 5m$, 
% and $$D_2 = 10m$ for different values of transmit power at H-AP, PA, in dBm. 
K = 2; % number of users
xi = 0.5; % energy harvesting efficiency
Gamma  = 9.8; % SNR gap, dB
D1 = 5; % distance from user i to the AP
D2 = 10;
alpha = 2; % channel pathloss exponent
sigma_1 = -160; % noise power spectral density, dBm/HZ;
sigma = -160 + 10*log10(10^6); % noise power, dBm;

%%
% Considering that it takes a long time for the iterative algorithm to
% converge, we do not average over 1000 randomly generated fading channel
% realizations as stated in the paper. Instead, we take the mean value of
% the fading, i.e., $\rho_i^2 = 1$, in this part.
h_1 = 10^(-3)*D1^(-alpha);
g_1 = 10^(-3)*D1^(-alpha);
h_2 = 10^(-3)*D2^(-alpha);
g_2 = 10^(-3)*D2^(-alpha);
PA = 0:5:30;
N_PA = length(PA);
gamma1 = xi*h_1*g_1*10.^(PA./10)/(10^((Gamma+sigma)/10));
gamma2 = xi*h_2*g_2*10.^(PA./10)/(10^((Gamma+sigma)/10));

R_o1_p1 = zeros(1, N_PA); R_o2_p1 = zeros(1, N_PA); R_sum_o2_p1 = zeros(1, N_PA);
R_sum_Norm_p1 = zeros(1, N_PA);
R_o1_p2 = zeros(1, N_PA); R_o2_p2 = zeros(1, N_PA); R_min = zeros(1, N_PA);
for i = 1 : N_PA
    % for P1, we can find the optimal tau based on Proposition 3.1
    gamma_p = [gamma1(i),gamma2(i)];
    [tau_p1] = Optimal_tau_func_pro3(K, gamma_p);
    R_o1_p1(i) =  tau_p1(2)*log2(1+gamma1(i)*tau_p1(1)/tau_p1(2));
    R_o2_p1(i) =  tau_p1(3)*log2(1+gamma2(i)*tau_p1(1)/tau_p1(3));
    R_sum_o2_p1(i) = R_o1_p1(i) + R_o2_p1(i);
    R_sum_Norm_p1(i) = R_sum_o2_p1(i)/K;
    
    % for P2, using the iterative algorithm to obtain the optimal solution
    [tau_opt, R_tau, lambda] = optimal_tau_cal(K, gamma_p);
    R_o1_p2(i) =  R_tau(1);
    R_o2_p2(i) =  R_tau(2);
    R_min(i) = min(R_tau);
end
figure
plot(PA,R_min,'-db',PA,R_sum_o2_p1,':or',PA,R_sum_Norm_p1,'-or',PA,R_o1_p1,'--or',PA,R_o2_p1,'-.or','linewidth',2);
xlabel('P_A(dBm)','fontsize',12);
ylabel('Average Thoughput (Mbps)','fontsize',12);
legend('Max.Common-Thoughput for(P2)','Max.Sum-Thoughput for(P1)','Normalized Max.Sum-Thoughput',...
    'Thoughput of U_1 in (P1)','Thoughput of U_2 in (P1)');
grid on;
title('Fig. 10 in the Paper (Iterative Algorithm)');

%%
% The achievable throughput is generally larger than that
% shown in the paper. A possible reason is that we ignored the impact
% of fast fading. To consider the impact of the fading, we can use the
% brute borce algorithm and average over 1000 randomly generated fading
% channel realizations, we can obtain almost identical figure as Fig 10 in
% the paper.

num_rand = 1000;
R_o1_p1 = zeros(num_rand, N); R_o2_p1 = zeros(num_rand, N); R_sum_o2_p1 = zeros(num_rand, N);
R_sum_Norm_p1 = zeros(num_rand, N);
R_min = zeros(num_rand, N);

for iter_rand = 1: num_rand
    rho = exprnd(1);
    
    h_1 = 10^(-3)*rho*D1^(-alpha);
    g_1 = 10^(-3)*rho*D1^(-alpha);
    h_2 = 10^(-3)*rho*D2^(-alpha);
    g_2 = 10^(-3)*rho*D2^(-alpha);
    gamma1 = xi*h_1*g_1*10.^(PA./10)/(10^((Gamma+sigma)/10));
    gamma2 = xi*h_2*g_2*10.^(PA./10)/(10^((Gamma+sigma)/10));
    
    for i = 1 : N_PA
        % for P1, we can find the optimal tau based on Proposition 3.1
        gamma_p1 = [gamma1(i),gamma2(i)];
        [tau_p1] = Optimal_tau_func_pro3(K, gamma_p1);
        R_o1_p1(iter_rand, i) =  tau_p1(2)*log2(1+gamma1(i)*tau_p1(1)/tau_p1(2));
        R_o2_p1(iter_rand, i) =  tau_p1(3)*log2(1+gamma2(i)*tau_p1(1)/tau_p1(3));
        R_sum_o2_p1(iter_rand, i) = R_o1_p1(iter_rand, i) + R_o2_p1(iter_rand, i);
        R_sum_Norm_p1(iter_rand, i) = R_sum_o2_p1(iter_rand, i)/2;
        
        % for P2, using the brute force algorithm to obtain the optimal solution
        [R_max, tau_opt] = Brute_force(K, gamma_p1, 0.02, 2);
        R_min(iter_rand, i) = min(R_max);
    end
end

figure
plot(PA,mean(R_min),'-db',PA,mean(R_sum_o2_p1),':or',PA,mean(R_sum_Norm_p1),'-or',PA,mean(R_o1_p1),'--or',PA,mean(R_o2_p1),'-.or','linewidth',2);
xlabel('P_A(dBm)','fontsize',12);
ylabel('Average Thoughput (Mbps)','fontsize',12);
legend('Max.Common-Thoughput for(P2)','Max.Sum-Thoughput for(P1)','Normalized Max.Sum-Thoughput',...
    'Thoughput of U_1 in (P1)','Thoughput of U_2 in (P1)');
grid on
title('Fig. 10 in the Paper (Brute Force Algorithm)');

%%
% By fixing $$PA = 20dBm$, Fig. 11 shows the throughput comparison for different values of the common pathloss exponent
% $$\alpha$ in both the DL and UL in the same WPCN as for Fig.10. 
PA = 20;
alpha = 2:0.5:4;
N_alpha = length(alpha);

% Channel power gains
h_1 = 10^(-3)*D1.^(-alpha);
g_1 = 10^(-3)*D1.^(-alpha);
h_2 = 10^(-3)*D2.^(-alpha);
g_2 = 10^(-3)*D2.^(-alpha);
gamma1 = xi.*h_1.*g_1*10^(PA/10)/(10^((Gamma+sigma)/10));
gamma2 = xi.*h_2.*g_2*10^(PA/10)/(10^((Gamma+sigma)/10));

R_o1_p1 = zeros(1, N_alpha); R_o2_p1 = zeros(1, N_alpha); R_sum_o2_p1 = zeros(1, N_alpha);
R_sum_Norm_p1 = zeros(1, N_alpha);
R_o1_p2 = zeros(1, N_alpha); R_o2_p2 = zeros(1, N_alpha); R_min = zeros(1, N_alpha);
for i = 1 : N_alpha
    gamma_p1 = [gamma1(i),gamma2(i)];
    % for P1, we can find the optimal tau based on Proposition 3.1
    [tau_p1] = Optimal_tau_func_pro3(K,gamma_p1);
    R_o1_p1(i) =  tau_p1(2)*log2(1+gamma1(i)*tau_p1(1)/tau_p1(2));
    R_o2_p1(i) =  tau_p1(3)*log2(1+gamma2(i)*tau_p1(1)/tau_p1(3));
    R_sum_o2_p1(i) = R_o1_p1(i) + R_o2_p1(i);
    R_sum_Norm_p1(i) = R_sum_o2_p1(i)/2;
    
    % for P2,  using the iterative algorithm to obtain the optimal solution
    [tau_opt, R_tau, lambda] = optimal_tau_cal(K, gamma_p1);
    R_o1_p2(i) =  R_tau(1);
    R_o2_p2(i) =  R_tau(2);
    R_min(i) = min(R_tau);
end

figure
plot(alpha,R_min,'-db',alpha,R_sum_o2_p1,':or',alpha,R_sum_Norm_p1,'-or',alpha,R_o1_p1,'--or',alpha,R_o2_p1,'-.or','linewidth',2);
xlabel('Pathloss Exponent(\alpha)','fontsize',12);
ylabel('Average Thoughput(Mbps)','fontsize',12);
legend('Max.Common-Thoughput for(P2)','Max.Sum-Thoughput for(P1)','Normalized Max.Sum-Thoughput',...
    'Thoughput of U_1 in (P1)','Thoughput of U_2 in (P1)');
grid on
title('Fig. 11 in the Paper (Iterative Algorithm)');
%%
% Similarly, the achievable performance is generally larger than the
% results given in the paper. To simulate the impact of channel fading, we
% can use the brute force algorithm and average over 1000 randomly generated fading
% channel realizations
num_rand = 1000;
R_o1_p1 = zeros(num_rand, N_alpha); R_o2_p1 = zeros(num_rand, N_alpha); R_sum_o2_p1 = zeros(num_rand, N_alpha);
R_sum_Norm_p1 = zeros(num_rand, N_alpha);
R_min = zeros(num_rand, N_alpha);

for iter_rand = 1: num_rand
    rho = exprnd(1);
    
    h_1 = 10^(-3)*rho*D1.^(-alpha);
    g_1 = 10^(-3)*rho*D1.^(-alpha);
    h_2 = 10^(-3)*rho*D2.^(-alpha);
    g_2 = 10^(-3)*rho*D2.^(-alpha);
    gamma1 = xi.*h_1.*g_1*10^(PA/10)/(10^((Gamma+sigma)/10));
    gamma2 = xi.*h_2.*g_2*10^(PA/10)/(10^((Gamma+sigma)/10));
    
    for i = 1 : N_alpha
        gamma_p1 = [gamma1(i),gamma2(i)];
        % for P1, we can find the optimal tau based on Proposition 3.1
        [tau_p1] = Optimal_tau_func_pro3(K,gamma_p1);
        R_o1_p1(iter_rand, i) =  tau_p1(2)*log2(1+gamma1(i)*tau_p1(1)/tau_p1(2));
        R_o2_p1(iter_rand, i) =  tau_p1(3)*log2(1+gamma2(i)*tau_p1(1)/tau_p1(3));
        R_sum_o2_p1(iter_rand, i) = R_o1_p1(iter_rand, i) + R_o2_p1(iter_rand, i);
        R_sum_Norm_p1(iter_rand, i) = R_sum_o2_p1(iter_rand, i)/2;
        
        % for P2, using the brute force algorithm to obtain the optimal solution
        [R_max, tau_opt] = Brute_force(K, gamma_p1, 0.02, 2);
        R_min(iter_rand, i) = min(R_max);
    end
end

figure
plot(alpha,mean(R_min),'-db',alpha,mean(R_sum_o2_p1),':or',alpha,mean(R_sum_Norm_p1),...
    '-or',alpha,mean(R_o1_p1),'--or',alpha,mean(R_o2_p1),'-.or','linewidth',2);
xlabel('Pathloss Exponent(\alpha)','fontsize',12);
ylabel('Average Thoughput(Mbps)','fontsize',12);
legend('Max.Common-Thoughput for(P2)','Max.Sum-Thoughput for(P1)','Normalized Max.Sum-Thoughput',...
    'Thoughput of U_1 in (P1)','Thoughput of U_2 in (P1)');
grid on
title('Fig. 11 in the Paper (Brute Force Algorithm)');
%%
% Fig. 12 shows the throughput over number of users, $$K$. It is assumed that $$K$ users in the network are
% equally separated from the H-AP according to $$D_i = \frac{D_K}{K} \times i$, $$i = 1, \ldots, K$, where $$D_K = 10m$. The transmit power at
% the H-AP and the pathloss exponent are set to be fixed as $$PA = 20dBm$ and $$\alpha = 2$, respectively. In addition,
% the authors also compare with the throughput achievable by equal time
% allocation (ETA), i.e., $$\tau_i = \frac{1}{K+1}, i = 0,1, \ldots, K$, as a lowcomplexity time allocation scheme.
% 
alpha = 2;
DK = 10;
PA = 20;
K_1 = 2:1:6;
N_K = length(K_1);

% calculate gamma with different numbers of users
gamma = cell(1, N_K);
for j = 1 : N_K
    gamma{j} = zeros(1, j);
    for z = 1 : j+1
        D = DK/(j+1)*z;
        h = 10^(-3)*D^(-alpha);
        g = 10^(-3)*D^(-alpha);
        gamma{j}(z) = xi*h*g*10^(PA/10)/(10^(Gamma/10)*10^(sigma/10));
    end
end
        
% Normalized Max sum-thoughput P1
R_Norm_max = zeros(1, N_K);
for j = 2:1:6
    % For max sum-throughput (P1), optimize based on Proposition 3.1
    [tau_p1] = Optimal_tau_func_pro3(j,gamma{j-1});
    R_o_p1 = zeros(1, j);
    for i=1:j
        R_o_p1(i) = tau_p1(i+1)*log2(1+gamma{j-1}(i)*tau_p1(1)/tau_p1(i+1));
    end
    R_Norm_max(j-1) = sum(R_o_p1)/j;
end

% Max. Common-throughput P2
R_max_common = zeros(1, N_K);
for j = 2:1:6
    % use the iterative algorithm to maximize the common throughput
    [tau_p2, R_tau, lambda] = optimal_tau_cal(j, gamma{j-1});
    R_o_p2 = zeros(1, j);
    for i=1:j
        R_o_p2(i) = tau_p2(i+1)*log2(1+gamma{j-1}(i)*tau_p2(1)/tau_p2(i+1));
    end
    R_max_common(j-1) = min(R_o_p2);
end

% Equal time allocation
R_min_eta = zeros(1, N_K);
R_norm_eta = zeros(1, N_K);
for i = 1 : N_K
    R = zeros(1, K_1(i));
    for z = 1 : K_1(i)
        tau = 1 / (K_1(i)+1);
        R(z) = tau * log2(1 + gamma{i}(z));
    end
    R_min_eta(i) = min(R);
    R_norm_eta(i) = sum(R)/K_1(i);
end

figure;
plot(K_1,R_Norm_max,'-or', K_1,R_norm_eta,'-sk', K_1,R_max_common,'-db', K_1,R_min_eta,':sk', 'linewidth',2);
xlabel('Number of users(K)','fontsize',12);
ylabel('Average Thoughput(Mbps)','fontsize',12);
legend('Normalized Max Sum-Thoughput','Normalized Sum-Thoughput by ETA','Max.Common-Thoughput','Min Thoughput by ETA');
grid on
title('Fig. 12 in the Paper (Iterative Algorithm)');
%%
% Similarly, we can have the resultes obtained by brute force algorithm.

num_rand = 1000;
R_Norm_max = zeros(num_rand, N_K); R_max_common = zeros(num_rand, N_K); 
R_min_eta = zeros(num_rand, N_K); R_norm_eta = zeros(num_rand, N_K);

for iter_rand = 1: num_rand
    rho = exprnd(1);
    
    % calculate gamma with different numbers of users
    gamma = cell(1, N_K);
    for j = 1 : N_K
        gamma{j} = zeros(1, j);
        for z = 1 : j+1
            D = DK/(j+1)*z;
            h = 10^(-3)*rho*D^(-alpha);
            g = 10^(-3)*rho*D^(-alpha);
            gamma{j}(z) = xi*h*g*10^(PA/10)/(10^(Gamma/10)*10^(sigma/10));
        end
    end

    % Normalized Max sum-thoughput P1
    for j = 2:1:6
        [tau_p1] = Optimal_tau_func_pro3(j,gamma{j-1});        
        R_o_p1 = zeros(1, j);
        for i=1:j
            R_o_p1(i) = tau_p1(i+1)*log2(1+gamma{j-1}(i)*tau_p1(1)/tau_p1(i+1));
        end
        R_Norm_max(iter_rand, j-1) = sum(R_o_p1)/j;
    end
    
    % Max. Common-throughput P2, using brute force algorithm 
    for j = 1 : N_K
        [R_max, tau_opt] = Brute_force(K_1(j), gamma{j}, 0.05, 2);
        R_max_common(iter_rand, j) = R_max;
    end
    
    % Equal time allocation
    for i = 1 : N_K
        R = zeros(1, K_1(i));
        for z = 1 : K_1(i)
            tau = 1 / (K_1(i)+1);
            R(z) = tau * log2(1 + gamma{i}(z));
        end
        R_min_eta(iter_rand, i) = min(R);
        R_norm_eta(iter_rand, i) = sum(R)/K_1(i);
    end
    
end

figure
plot(K_1,mean(R_Norm_max),'-or', K_1,mean(R_norm_eta),'-sk', K_1,mean(R_max_common),'-db', K_1,mean(R_min_eta),':sk', 'linewidth',2);
xlabel('Number of users(K)','fontsize',12);
ylabel('Average Thoughput(Mbps)','fontsize',12);
legend('Normalized Max Sum-Thoughput','Normalized Sum-Thoughput by ETA','Max.Common-Thoughput','Min Thoughput by ETA');
grid on
title('Fig. 12 in the Paper (Brute Force Algorithm)');

%% Analysis and Conclusions 
% 1. The two proposed maximization problems, i.e.,sum-throughput
% maximization problem (P1) and common-throughput maximization problem
% (P2), are reproduced. We were able to reprodece all of the results figures presented in the original paper [1]. 
%
% 2. For problem (P1), the authors have given out a closed form of the solution and shown the corresponding proofs. 
% In this report, we also use the brute force method to solve the original
% problem (P1) to further show the correctness of the provided closed 
% form expression. For problem (P2), the authors have mentioned to use sub-gradient
% based algorithms to update $$\lambda$ in their proposed algorithm, thus,
% we choose the gradient descent method to realize the provided algorithm in this report. Our 
% simulation results have shown the correctness of the original simulation
% results. 
%
% 3. However, for problem (P2), we noticed that the running time of the proposed iterative algorithm is $$10$ or more times longer than
% the brute force method. One-time simulation of the proposed iterative algorithm lasts for hours, while only 10 minutes are required for the
% brute force method. Thus, for Fig.10 - 12, instead of showing the results by averaging over 1000 randomly generated fading channel 
% realizations as the authors mentioned in the paper, we have shown the results of one-time simulation of the proposed iterative algorithm and also shown the results of using brute force algorithm and
% averaging over 1000 times. Moreover, we also noticed that, the iterative
% algorithm implemented in this project converges to a better result than
% that given in the paper, although it may take a longer time.
%
%
%% References
% [1] Hyungsik Ju and Rui Zhang, `Throughput Maximization in Wireless
% Powered Communication Networks', IEEE Trans. Wireless Commun., vol.13,
% no. 1, pp. 418-428, 2014.

%% Print all the functions used in the project
type('Optimal_tau_func_pro3')
%%
type('optimal_tau_cal')
%%
type('tau_opt_cal')
%%
type('gradient_descent')
%%
type('Brute_force')
