% File: V_imperfect_sense_wi_ACK.m
% Objective: 
% Find the optimal sensing-transmission control for opportunistic spectrum 
% access of cognitive radio networks using dynamic programming
% Author: Senhua Huang
% Department of Electrical and Computer Engineering
% UC Davis
% Email: senhua AT ece.ucdavis.edu
% Date: 07/14/2008
% All rights reserved.
% If you have any questions or comments, please contact me. Thanks.

% More descriptions:
% Get the optimal performance and the optimal threshold for the access
% control when the idle time follows UNIFORM distribution
% Other distributions can be tested by trivial modifications
% Imperfect Sensing modeled as (Pf, Pd) pairs

clear;
clc;

% Simulation Parameters setting
% Cost of collision
C = 5;

% Captured effect + imperfect collision detection
% a0: Pr[NACK| No collision at PU]
% a1: Pr[NACK| collision at PU]
a0 = 0.1;
a1 = 0.9;

% Imperfect Sensing
% Pf: Prob. of false-alarm, recognizing idle PU as busy
% Pd: Prob. of missed detection, recognizing busy PU as idle
Pf = 0.01;
Pd = 0.99;

% average idle time of PU
vp = 500;
% PU idle time uniformly distributed over [0, Imax]
Imax = 2*vp;
% average busy time of PU
lp = 500;

% Length of SU transmission packet
T = 30;
% Reward for successful transmission
R = 1;
% Length of SU sensing time
Ts = 10; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the optimal policy: threshold-based policy
% belief value of SU: Prob[PU idle @ time t]
% quantified
p_set = [0:0.0001:1];
len_p = length(p_set);
p_ind = [1:1:len_p];

% Time index
k_set = [1:1:Imax];
len_k = length(k_set);

% Transition probability at time k with regard to sensing
% Specific to different idle time distribution
g_ks = (Imax - (k_set+Ts-1))./(Imax - (k_set-1));
g_ks(find(g_ks<0)) = 0;

%transition probability at time k with regard to transmission
% Specific to different idle time distribution
g_kt = (Imax - (k_set+T-1))./(Imax - (k_set-1));
g_kt(find(g_kt<0)) = 0;

% Technical constraint: critical value
cc = (C+a1-1)/(a1-a0+C);

% vector that stores the value function for each p, and each k
Pro = 30;%protection for overflow
V = zeros(len_k+Pro,len_p+1);
% vector that stores the optimal action for each p, and each k
F = zeros(1,len_p);

%T_star: minimum time that the optimal action is to sense even when the
%probability of primary being idle is 1
if(cc>0)
    T_star = max(find(g_kt>=cc));
else
    T_star = Imax;
end

T_star = max(find(g_kt>=cc));

% vector that stores the threshold for each time index
pp_star = ones(1,len_k);


for ii = 1:T_star-1,
    t = T_star - ii;
    for pp = 1:len_p,
        p = p_set(pp);
        
        % probability of sensing idle
        p_idle = (p*g_ks(t)*(1-Pf)+(1-p*g_ks(t))*(1-Pd));
        % update belief value if sensing idle
        if(p_idle),
           pl_idle = p*g_ks(t)*(1-Pf)/p_idle; 
        else % p_idle = 0
            % does not have impact on value function
            pl_idle = 0.9; 
        end
        % probability of sensing busy
        p_busy = 1-p_idle;
        % update belief value if sensing busy
        if(p_busy),
           pl_busy =p*g_ks(t)*Pf/p_busy;
        else
            pl_busy = 0.9;
        end
        % quantization of the updated belief
        % note: other ways (e.g., interpolation is also feasible here)
        pl_busy_ind = floor(pl_busy*len_p) + 1;
        pl_idle_ind = floor(pl_idle*len_p) + 1;
        
        % calculate the expected reward of sensing action
        r_s = p_idle*V(t+Ts,pl_idle_ind) + p_busy*V(t+Ts,pl_busy_ind);
        
        % update belief value if transmit a packet
        pl_ind = floor(p*g_kt(t)*len_p)+1;
        
        % probability of receiving an ACK at SU-Rx
        p_ack = (p*g_kt(t)*(1-a0)+(1-p*g_kt(t))*(1-a1));
        % update belief value if an ACK is received 
        if(p_ack ==0),
            pl_ack = 0.9;
        else
            pl_ack = p*g_kt(t)*(1-a0)/p_ack;
        end

        % probability of receiving a NACK at SU-Rx
        p_nack = 1-p_ack;
        % update belief value if a NACK is received
        if(p_nack == 0),
            pl_nack = 0.9;
        else
            pl_nack =p*g_kt(t)*a0/p_nack;
        end

        % quantization of the updated belief
        pl_ack_ind = floor(pl_ack*len_p)+1;
        pl_nack_ind = floor(pl_nack*len_p)+1;

        % calculate the immediate reward of transmission
        im_r_t = p*g_kt(t)*R*(1-a0)+(1-p*g_kt(t))*(1-a1-C);
        % calculate the expected reward of transmission
        r_t = im_r_t*T + p_ack*V(t+T,pl_ack_ind) + p_nack*V(t+T,pl_nack_ind);
        % Bellman equation
        V(t,pp) =max(r_s,r_t);
        % store the optimal action
        if(r_s>=r_t)
            F(pp) = 0;
        else
            F(pp) = 1;
        end
    end
    
    % get the optimal threshold
    % for uniform distribution, we probably have at most one threshold
    if(sum((F(:)>0))), % For some value of p, optimal action at t is to transmit
        star_ind = min(find(F(:)>0));
        pp_star(t) = p_set(star_ind);
    else
        pp_star(t) = 1.0;
    end

end
% Maximum average reward per time unit for the SU 
reward_V = V(1,len_p)/(vp+lp);
% plot the threshold
figure(1);clf;
set(gca,'FontSize',16);
plot(pp_star);
xlabel('Time index');
ylabel('p^*');