generateFigure49.m

%This m-file generates Figure 4.9 in "Optimal Resource Allocation in
%Coordinated Multi-Cell Systems" by Emil Björnson and Eduard Jorswieck.
%
%This is version 1.2. (Last edited: 2014-07-23)
%
%Figure 4.9 illustrates the performance regions for random channel
%realizations in a scenario with two users and global joint transmission.
%The main difference from Figure 3.1 is the transceiver impairments, which
%are modeled with linear distortion functions. The Pareto boundary is
%generated by solving FPOs over a grid of weighting vectors, by combining
%Theorem 3.4 and Corollary 4.3. The EVM at the transmitters and receivers
%is 0, 5, 10, or 15 percent.
%
%Please note that the channels are generated randomly, thus the results
%will not be exactly the same as in the book.

close all;
clear all;


%%Simulation parameters

Kt = 2; %Number of base stations (BSs)
Kr = 2; %Number of users (in total)
Nt = 2; %Number of antennas per BS

%Define the 4 EVM combinations that will be simulated. In this case, the
%transmitter and receiver always have the same EVMs.
EVMrangeTxAndRx = [0 5 10 15]/100;

rng('shuffle'); %Initiate the random number generators with a random seed

%%If rng('shuffle'); is not supported by your Matlab version, you can use
%%the following commands instead:
%randn('state',sum(100*clock));
%rand('state',sum(100*clock));

PdB = 10; %SNR (in dB)
P = 10.^(PdB/10); %SNR

%D-matrices for global joint transmission
D = repmat(eye(Kt*Nt),[1 1 Kr]);

%Definition of per base station power constraints
L = Kt;
Q = zeros(Kt*Nt,Kt*Nt,L); %The L weighting matrices
Qsqrt = zeros(Kt*Nt,Kt*Nt,L); %The matrix-square root of the L weighting matrices
for l = 1:L
    Q((l-1)*Nt+1:l*Nt,(l-1)*Nt+1:l*Nt,l) = (1/P)*eye(Nt);
    Qsqrt((l-1)*Nt+1:l*Nt,(l-1)*Nt+1:l*Nt,l) = sqrt(1/P)*eye(Nt);
end
q = ones(L,1); %Limits of the L power constraints. Note that these have been normalized.

%Generation of normalized Rayleigh fading channel realizations with unit
%variance to closest base station and 1/3 as variance from other base station.
pathlossOverNoise = 2/3*eye(Kr)+1/3*ones(Kr,Kr);
H = sqrt(kron(pathlossOverNoise,ones(1,Nt))).*(randn(Kr,Nt*Kt)+1i*randn(Kr,Nt*Kt))/sqrt(2);




%%Calculate sample points on the Pareto boundary by applying
%%Theorem 3.4 on fine grid of weighting vectors.

delta = 0.01; %Accuracy of the line searches in the FPO algorithm

nbrOfSamplesFPO = 1001; %Number of sample points
ParetoBoundary = zeros(Kr,nbrOfSamplesFPO,length(EVMrangeTxAndRx)); %Pre-allocation of sample matrix


%Computation of the utopia point using MRT with full transmit power, 
%which is optimal when there is only one active user.

wMRT = functionMRT(H,D); %Compute normalized beamforming vectors for MRT

utopia = zeros(Kr,1);
for k = 1:Kr
    W = sqrt(P)*[wMRT(1:2,k)/norm(wMRT(1:2,k)); wMRT(3:4,k)/norm(wMRT(3:4,k))]; %Scale to use all available power
    utopia(k) = log2(1+abs(H(k,:)*W)^2);
end


%Generate a grid of equally spaced search directions
interval = linspace(0,1,nbrOfSamplesFPO);
profileDirections = [interval; 1-interval];

for m = 1:nbrOfSamplesFPO

    %Output of the simulation progress - since this part takes a long time
    if mod(m,10) == 1
        disp(['Progress: Generating boundary, ' num2str(m) ' out of ' num2str(nbrOfSamplesFPO) ' search directions']); 
    end
    
    %Search on a line from the origin in each of the equally spaced search
    %directions. The upper point is either on Pareto boundary (unlikely) or
    %outside of the rate region.
    lowerPoint = zeros(Kr,1);
    upperPoint = sum(utopia) * profileDirections(:,m);
    
    %Find the intersection between the line and the Pareto boundary by
    %solving an FPO problem. Different input is given if the hardware is
    %ideal or has transceiver impairments.
    for impairIndex = 1:length(EVMrangeTxAndRx)
        
        if EVMrangeTxAndRx(impairIndex) == 0
            finalInterval = functionFairnessProfile_cvx(H,D,Qsqrt,q,delta,lowerPoint,upperPoint);
        else
            finalInterval = functionFairnessProfile_cvx(H,D,Qsqrt,q,delta,lowerPoint,upperPoint,1,EVMrangeTxAndRx(impairIndex)*ones(2,1));
        end
        
        ParetoBoundary(:,m,impairIndex) = finalInterval(:,1); %Save the best feasible point found by the algorithm
    end
    
    
end


%Plot the rate regions for different levels of impairments
figure; hold on; box on;

for impairIndex = 1:length(EVMrangeTxAndRx)
    plot(ParetoBoundary(1,:,impairIndex),ParetoBoundary(2,:,impairIndex),'k-','LineWidth',1); %Plot rate region
    
    %Search for the sample points that maximize the sum rate, geometric mean,
    %and max-min fairness, respectively. Note that the accuracy
    %of these points improves when the number of sample points is increased.
    [~,sumrateIndex] = max(ParetoBoundary(1,:,impairIndex)+ParetoBoundary(2,:,impairIndex));
    [~,geomeanIndex] = max(ParetoBoundary(1,:,impairIndex).*ParetoBoundary(2,:,impairIndex));
    [~,maxminIndex] = max(min(ParetoBoundary(:,:,impairIndex),[],1));

    %Plot sample points that maximize sum rate, geometric mean,
    %harmonic mean, and max-min fairness, respectively.
    plot(ParetoBoundary(1,sumrateIndex,impairIndex),ParetoBoundary(2,sumrateIndex,impairIndex),'ro');
    plot(ParetoBoundary(1,geomeanIndex,impairIndex),ParetoBoundary(2,geomeanIndex,impairIndex),'sg');
    plot(ParetoBoundary(1,maxminIndex,impairIndex),ParetoBoundary(2,maxminIndex,impairIndex),'k*');
end

legend('Pareto Boundary based on FPO','Max Arithmetic Mean','Max Geometric Mean','Max-Min Fairness','Location','SouthWest')

xlabel('log_2(1+SINR_1) [bits/channel use]')
ylabel('log_2(1+SINR_2) [bits/channel use]')