% DPC_wo_M_sync.m
% Distributed power control simulation without the knowledge of exact number of
% secondary user pairs.
% With the same initial value of lambda_i
% Assuming each SU has been able to synchronize for their updating
% Date: 2009/07/25
% Author: Senhua Huang
% Email: senhua@ece.ucdavis.edu
% Address: Dept. of ECE, UC Davis, Davis, CA 95616 USA
% All rights reserved

clear;
clc;

% Simulation Parameter
rand('state',0);
% Lenght of simulation
lenSim = 1000000; 

% Normalized received signal power at the PU-Rx 
a = 1;   
% SINR threshold for a successful reception of the PU packet
cinr_th = 6;  
% Outage probability constraint
eta = 0.1;

% Noise power level normalized to 1 (-100dBm)
N0 = 1; 
%path loss exponent
alpha = 4; 
%number of SU pairs
M = 4;

% transmit power of the PU-Tx: 2000mW, i.e., 33dBm
powerPUTX = 2000; 

% maximum transmit power of the SU-Tx
maxPowerSU = 10000; 
minPowerSU = 0;

% Topology
% for normalized noise power -100dBm -> 0dB
d0 = (10^(10/alpha)); 
x = 3000/d0;
y = 3000/d0;

rangePU = 500;
x_PUTX = x/2;
y_PUTX = y/2;
x_PURX = x_PUTX + rangePU/d0;
y_PURX = y_PUTX;
d_PU = sqrt((x_PUTX-x_PURX)^2 + (y_PUTX-y_PURX)^2);

rangeSU = 100;
x_SUTX(1:M)= x * rand(M,1);
y_SUTX(1:M)= y * rand(M,1);
x_SURX(1:M) = x_SUTX + rangeSU/d0;
y_SURX(1:M) = y_SUTX + rangeSU/d0;

d_S2P(1:M) = sqrt((x_SUTX-x_PURX).^2 + (y_SUTX-y_PURX).^2);
d_P2S(1:M) = sqrt((x_SURX-x_PUTX).^2 + (y_SURX-y_PUTX).^2);
d_SU(1:M) = sqrt((x_SUTX-x_SURX).^2 + (y_SUTX-y_SURX).^2);

% Trace of SU powers
log_powerSU = zeros(M,lenSim);

%channel parameter 
G0 = d_PU^(-alpha);
P0 = powerPUTX;
% outage probability of PU without any SU transmission
eta0 = 1-exp(-N0*cinr_th/P0/G0);
% relative interference margin
mu = (1-eta0)/(1-eta);
 
% SU channels
Hi = (d_SU(1:M)'.^(-alpha))./(N0+powerPUTX*(d_P2S(1:M)'.^(-alpha)));
Gi = d_S2P(1:M).^(-alpha);
bi = Gi*cinr_th./(P0*G0);
bi = bi';

% Solve problem using the closed-form solution assuming knowledge of the
% realization of M (M)
 
opt_lambda =1/(mu^(1/M)-1) + 1;
opt_p = (mu^(1/M)-1)./bi;
opt_util = sum(log(1+Hi.*opt_p));

% Distributed Power Control Algorithm
% Solve using dual decomposition with pricing update done by a genie
prim_obj = [];dual_obj = []; infeas = []; price = [];

% initial value of stepsize
a0 = 20;
lambda = 100;
powerSU = zeros(M,1);

% Length of learning period with fixed transmit power 
T = 100;

% # of learning periods
NN = floor(lenSim/T);

% outage prob. estimate at each secondary user
ro = zeros(1,M);

for i = 1:NN,  
    % transmit with the following power
    for j = 1:M,
        powerSU(j) = 1/((lambda-1)*bi(j));
        if (powerSU(j)<0),
            powerSU(j) = 0;
        end
        if(powerSU(j)>maxPowerSU),
            powerSU(j) = maxPowerSU;
        end
    end
    log_powerSU(1:M,i+1) = powerSU;

    prim_obj = [prim_obj sum(log(1+powerSU.*Hi))];

    % Calculate the value of outage probability
    prob_h(i) =  1 - (1-eta0)*1/(exp(sum(log(1+bi.*powerSU))));
    rh(i) = prob_h(i);
    rh(i) = sum((rand(1,T)<prob_h(i)))/T; %#ok<AGROW>

    if(rh(i)==0),
      rh(i) = 1/T;
    end
    
    % Dual subgradient
    lambda = lambda + a0/sqrt(i)*(log(1-eta)-log(1-rh(i)));
    price = [price lambda];
    
end

close all;
plot_ind = [1:1:NN];
figure(10), clf
set(gca,'FontSize',16);
plot(plot_ind,prim_obj,'k-',plot_ind,ones(1,NN)*opt_util,'b-','LineWidth',1.5);
legend('Simulation','Optimal');
xlabel('time k'), ylabel('objective')
print -depsc objective_evol_T10.eps

figure(20), clf
set(gca,'FontSize',16);
plot(plot_ind,prob_h,'k-',plot_ind,ones(1,NN)*eta,'b-','LineWidth',1.5);
xlabel('time k'), ylabel('Outage prob.')
legend('Simulation','Constraint');
print -depsc Pout_evol_T10.eps

figure(30), clf
set(gca,'FontSize',16);
plot(plot_ind,price,'k-',plot_ind,ones(1,NN)*opt_lambda,'b-','LineWidth',1.5);
xlabel('time k'), ylabel('price')
legend('Simulation','Optimal');
print -depsc price_evol_T10.eps