%% ECE602 Project: A Robust Maximin Approach for MIMO Communications With Imperfect Channel State Information Based on Convex Optimization [1]
%
% Student Name: ALONG JIN (20687178)
%
%% *1. Background of MIMO Systems and Related Work*
% 
% The multiple-input-multiple-output (MIMO) systems have been widely
% studied and used to achieve the spatial diversity and improve the quality
% of wireless communications. The performance of MIMO systems depend
% heavily on the quantity and quality of the channel state information
% (CSI) at both the transmitters and receivers, and such channel knowledge
% is generally imperfect in realistic scenarios.
%
% Let the $\mathrm{\hat{H}}$ be the estimation of the imperfect channel. In
% the non-robust design, only the maximum eigenmode of $\mathrm{\hat{H}}^H
% \mathrm{\hat{H}}$ is used. On the other hand, the pure orthogonal
% space-time block code (OSTBC) approach uses all the eigenmodes, over
% which the transmit power is uniformly allocated. However, without taking
% into account explicitly the errors in the CSI estimation, both the
% non-robust and pure OSTBC systems suffer from severe performance
% degradation with increasing error level.
%
% In this paper, the authors propose a robust design for MIMO communication
% system by explicitly considering the errors in the channel estimate.
% Particularly, the robust design is based on the notion of maximin
% philosophy. The authors first give the system model and formulate an
% maximin problem. Then, by assuming the convex uncertainty regions of
% the errors, the optimization problem is transformed to a simple convex
% problem. Besides, different convex uncertainty regions are investigated,
% in terms of the Gaussian noise and quantization error.
%
%
%
%% *2. System Model and Problem Formulation*
%
% Consider a MIMO system with $n_T$ transmit and $n_R$ receive antennas,
% where the actual channel matrix is $\mathrm{H}$. Let $\Delta$ be the
% error in the estimate, there is 
%
% $$ \mathrm{H} = \hat{\mathrm{H}} + \Delta \quad (1) $$
%
% To design robust transmitter, the authors adopt OSTBC block and
% distribute transmit power over the all the eigenmodes of
% $\mathrm{\hat{H}}^H \mathrm{\hat{H}}$. Let $\mathrm{\hat{U}} =
% [\mathrm{\hat{u}_1} \cdots \mathrm{\hat{u}_{n_T}}]$ be the normalized
% eigenvectors of $\hat{\mathrm{H}}^H \hat{\mathrm{H}}$, and the
% corresponding power allocation be $\mathbf{p}^T = [p_1 \cdots p_{n_T}]$.
% The performance of proposed system can be measured by the SNR, given by
%
% $$ \mathrm{SNR} = \frac{1}{\sigma_n^2} \mathrm{Tr}(\hat{\mathrm{U}}^H
% \mathrm{H}^H \mathrm{H \hat{U} diag(\mathbf{p})}) \quad (2) $$
%
% where $\sigma_n^2$ is the power of additive white Gaussian noise (AWGN)
% at the receiver. Based on this, the performance function $f$ in this
% system can be defined as
%
% $$ f(\mathbf{p}, \Delta) = \mathrm{Tr}(\hat{\mathrm{U}}^H \mathrm{H}^H
% \mathrm{H \hat{U} diag(\mathbf{p})}) = \mathrm{Tr}(\hat{\mathrm{U}}^H
% \mathrm{(\hat{H} + \Delta)}^H \mathrm{(\hat{H} + \Delta) \hat{U}
% diag(\mathbf{p})}) \quad (3) $$
%
% The objective is to optimize the power allocation $\mathbf{p}$ with the
% worst case error $\Delta$ in the uncertainty region $\mathcal{R}$. Then,
% the robust approach is formulated as
%
% $$ \max_{\mathbf{p}} \inf_{\Delta \in \mathcal{R}} f(\mathbf{p}, \Delta)
% \quad \mathrm{s.t.} \quad \mathbf{1}^T \mathbf{p} \le P_0, \quad p_i \ge
% 0 ~~\forall i \quad (4) $$
%
%
%
%% *3. Solving the Proposed Maximin Problem*
%
% Assuming that the constraint set $\mathcal{R}$ is convex, the maximin
% problem in (4) can be transformed into the following simplified convex
% optimization problem
%
% $$ \min_{t, \Delta \in \mathcal{R}} t \quad \mathrm{s.t.} \quad t \ge P_0
% \mathrm{\hat{u}_i}^H \mathrm{(\hat{H} + \Delta)}^H \mathrm{(\hat{H} +
% \Delta)} \mathrm{u_i} ~~\forall i  \quad (5) $$
%
% where $t$ is a dummy variable. The optimal robust power allocation
% $\mathbf{p}^* = [p_1^* \cdots p_{n_T}^*]^T$ is closely related to the
% optimal dual variables $\{ \gamma_i^* \}$ associated with the convex
% optimization problem in (5), i.e., $p_i^* = P_0 \gamma_i^*$.
%
% It should be noted that the optimal dual variables $\{ \gamma_i^* \}$ and
% worst case error $\Delta^*$ do not depend on the power budget $P_0$.
% Thus, to achieve a target $\mathrm{SNR}_0$, the minimum transmit power is
% given by
%
% $$ P_0^* = \mathrm{SNR}_0
% \frac{\sigma_n^2}{\mathrm{Tr}(\mathrm{\hat{U}}^H \mathrm{(\hat{H} +
% \Delta^*)}^H \mathrm{(\hat{H} + \Delta^*) \hat{U} diag(\{\gamma_i^*\})})}
% \quad (6) $$
%
%
%
%% *4. Convex Uncertainty Regions*
%
% In the following, two sources of errors are identified and described, as
% well as the joint uncertainty region.
%
%% $\quad$ 4-1. Estimation Gaussian Noise
%
% Let the unbiased channel estimation be $\mathrm{\tilde{H}} =
% \mathrm{H}+\mathrm{E}$, where $\mathrm{E}$ is the zero-mean estimation
% noise, independent from $\mathrm{H}$. If the channel matrix $\mathrm{H}$
% and error matrix $\mathrm{E}$ have i.i.d. components with zero-mean and
% variance $\sigma_h^2$ and $\sigma_e^2$, respectively, the uncertainty
% region for the channel estimation reduces to a sphere of radius
% $\sqrt{\epsilon}$, given by
%
% $$ \mathcal{R} = \{ \Delta: ||\Delta||_F^2 \le \epsilon, ~\epsilon = r^2
% \frac{\sigma_h^2}{1+\mathrm{SNR_{est}}} \} \quad (7) $$
%
% where $\mathrm{SNR_{est}} = \sigma_h^2 / \sigma_e^2$. Denote $P_{in} =
% \mathrm{Pro}(\Delta \in \mathcal{R})$, which is the probability of
% providing required QoS, i.e., probability of having a SNR higher than the
% target $\mathrm{SNR}_0$. The relationship between $P_{in}$ and $r^2$ is
% given by
%
% $$ P_{in} = \phi (r^2) \quad (8) $$
%
% where $\phi (\cdot)$ is the cumulative density function of the chi-square
% distribution with $2 n_R n_T$ degrees of freedom.
%
% Besides, for spherical uncertainty regions, a closed-form solution is
% provided, as shown in the subsection IV-E of the paper [1].
%
%
%% $\quad$ 4-2. Quantization Error
%
% Due to finite capacity of feedback channel, the channel response has to
% be quantized. Assuming that uniform quantization is used with the step
% $\Delta_q$, the quantization SNR is defined as $\mathrm{SNR}_q = 6
% \sigma_h^2 / \Delta_q^2$. The uncertainty region for the channel
% estimation is a hypercube centered at $\mathrm{\hat{H}}$, given by
%
% $$ \mathcal{R} = \{ \Delta: |\mathrm{Re}\{[\Delta]_{ij}\}| \le
% \frac{\Delta_q}{2}, ~|\mathrm{Im}\{[\Delta]_{ij}\}| \le
% \frac{\Delta_q}{2} \} \quad (9) $$
%
%
%% $\quad$ 4-3. Combined Estimation and Quantization Errors
%
% In realistic scenarios, both Gaussian noise and quantization error effect
% the channel estimation. Thus, an appropriate uncertainty region for the
% errors can be defined as
%
% $$ \mathcal{R} = \{ \Delta = \Delta_1 + \Delta_2: ||\Delta_1||_F^2 \le
% \epsilon, ~|\mathrm{Re}\{[\Delta_2]_{ij}\}| \le \frac{\Delta_q}{2},
% ~|\mathrm{Im}\{[\Delta_2]_{ij}\}| \le \frac{\Delta_q}{2} \} \quad (10) $$
%
% According to the region, the optimization problem in (5) can be rewritten
% as the following convex quadratic problem:
%
% $$ \min_{t, \Delta_1, \Delta_2} t \quad \mathrm{s.t.} \quad t \ge P_0
% \mathrm{\hat{u}_i}^H \mathrm{(\hat{H} + \Delta_1 + \Delta_2)}^H
% \mathrm{(\hat{H} + \Delta_1 + \Delta_2)} \mathrm{u_i} ~~\forall i, \quad
% \mathrm{Tr}(\Delta_1^H \Delta_1) \le \epsilon, \quad
% |\mathrm{Re}\{[\Delta_2]_{ij}\}| \le \frac{\Delta_q}{2}, \quad
% |\mathrm{Im}\{[\Delta_2]_{ij}\}| \le \frac{\Delta_q}{2}  \quad (11) $$
%
%
%
%% *5. Repreducing Simulation Results*
%
% In the following, the main results in the paper [1] will be reproduced
% with CVX in MATLAB.
%
%% $\quad$ 5-1. Mean Value of Normalized Power Allocation (Gaussian Noise with Spherical Uncertainty Regions)
%
%   close all;
% clear all;
% %
% % Parameters
% varH = 2; % variance of channel matrix H
% varE = 1; % variance of error matrix E
% nT = 4; % number of transmit antennas
% nR = 6; % number of receiver antennas
% g = 0:0.1:0.9;
% g = [g 0.95 0.99]; % parameter g
% %
% % Optimization problem with CVX (Eq.(23) in paper)
% N = 300; % number of iterations
% pDist = zeros(nT,12); % normalized power allocation for different radius
% for a = 1:N
%     H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%     E = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varE/2); % estimation error matrix
%     estH = (H+E)*varH/(varH+varE); % estimated CSI matrix
%     [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
%     r = norm(estH(:))*g; % radius of spherical uncertainty region
%     for b = 1:12
%         cvx_begin
%             variable t;
%             variable delta(nR,nT) complex;     
%             dual variables gamma{nT};
%             minimize (t);
%             norm(delta(:)) <= r(b);
%             for c = 1:nT
%                 gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
%             end
%         cvx_end
%         pDist(:,b) = pDist(:,b)+cell2mat(gamma);
%     end
% end
% pDist = pDist/N; % mean value of normalized power allocation
% save('data/powerDist', 'g', 'pDist'); % save data

% Fig.4 in Paper [1]
load('data/powerDist');
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 16);
hold on;
grid on;
plot(g, pDist(4,:), '-o', g, pDist(3,:), '--x', g, pDist(2,:), '-.*', g, pDist(1,:), ':^', 'linewidth', 2);
legend('Mean value of \gamma_1', 'Mean value of \gamma_2', 'Mean value of \gamma_3', 'Mean value of \gamma_4');
xlabel('Parameter g');
ylabel('Mean Value of \gamma_i');
title('Fig. 4 in Paper [1]');

% Code: powerDist.m


%% $\quad$ 5-2. Mean Value of Minimum Transmit Power (dB) - TDD System - Spherical Uncertainty Regions
%
%   close all;
% clear all;
% %
% % Parameters
% SNR0 = 10; % target SNR
% varG = 1; % variance of AWGN at receiver
% varH = 2; % variance of channel matrix
% estSNRdB = 8:14; % estimation SNR (dB)
% estSNR = 10.^(estSNRdB/10); % estimation SNR
% nT = 2; % number of transmit antennas
% nR = 2; % number of receive antennas
% %
% % Algorithms
% P0 = zeros(4,7); % minimum required transmit power
% N = 2; % number of iterations
% for a = 1:N
%     H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%     E = randn(nR,nT) + i*randn(nR,nT); % estimation error matrix without specifying variance
%     %    
%     qosP = 0.85; % probability of providing required QoS
%     r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR)); % radius of spherical uncertainty region
%     for b = 1:7
%         varE = varH/estSNR(b); % variance of error matrix
%         estH = (H+E*sqrt(varE/2))*varH/(varH+varE); % estimated CSI matrix
%         [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
%         %        
%         % Robust with all eigenmodes
%         cvx_begin
%             variable t;
%             variable delta(nR,nT) complex;     
%             dual variables gamma{nT}; % dual variables
%             minimize (t);
%             norm(delta(:)) <= r(b);
%             for c = 1:nT
%                 gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
%             end
%         cvx_end
%         P0(1,b) = P0(1,b)+SNR0*varG/real(cell2mat(gamma)'*diag(U'*(estH+delta)'*(estH+delta)*U));
%         %    
%         % Non-robust with maximum eigenmode
%         cvx_begin
%             variable delta(nR,nT) complex;
%             minimize (norm((estH+delta)*U(:,nT), 'fro'));
%             norm(delta(:)) <= r(b);
%         cvx_end
%         P0(2,b) = P0(2,b)+SNR0*varG/real(cvx_optval^2);
%         %    
%         % OSTBC with equal power over all eigenmodes
%         cvx_begin
%             variable delta(nR,nT) complex;
%             f = 0;
%             for c = 1:nT
%                 f = f+quad_form((estH+delta)*U(:,c), diag(ones(nR,1)))/nT;
%             end
%             minimize (f);
%             norm(delta(:)) <= r(b);
%         cvx_end
%         P0(3,b) = P0(3,b)+SNR0*varG/real(cvx_optval);
%     end
%     %    
%     % Robust for qosP = 0.6
%     qosP = 0.6; % probability of providing required QoS
%     r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR)); % radius of spherical uncertainty region
%     for b = 1:7
%         varE = varH/estSNR(b); % variance of error matrix
%         estH = (H+E*sqrt(varE/2))*varH/(varH+varE); % estimated CSI matrix
%         [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
%         cvx_begin
%             variable t;
%             variable delta(nR,nT) complex;     
%             dual variables gamma{nT}; % dual variables
%             minimize (t);
%             norm(delta(:)) <= r(b);
%             for c = 1:nT
%                 gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
%             end
%         cvx_end
%         P0(4,b) = P0(4,b)+SNR0*varG/real(cell2mat(gamma)'*diag(U'*(estH+delta)'*(estH+delta)*U));
%     end
% end
% P0 = 10*log10(P0/N); % minimum required transmit power (dB)
% save('data/minPowerTDD.mat', 'estSNRdB', 'P0'); % save data

% Fig. 6 in Paper [1]
load('data/minPowerTDD');
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 16);
hold on;
grid on;
plot(estSNRdB, P0(1,:), '-o', estSNRdB, P0(2,:), '--o', estSNRdB, P0(3,:), '-.o', estSNRdB, P0(4,:), '-*', 'linewidth', 2);
legend('Robust: P_{in}=0.85, n_T=2, n_R=2', 'Non-Robust: P_{in}=0.85, n_T=2, n_R=2', ...
    'OSTBC: P_{in}=0.85, n_T=2, n_R=2', 'Robust: P_{in}=0.6, n_T=2, n_R=2');
xlabel('Estimation SNR (dB)');
ylabel('Mean Value of Minimum Transmit Power (dB)');
title('Fig. 6 in Paper [1]');

% Code: minPowerTDD.m


%% $\quad$ 5-3. Mean Value of Minimum Transmit Power (dB) - FDD System - Cubic Uncertainty Regions
%
%   close all;
% clear all;
% %
% % Parameters
% SNR0 = 10; % target SNR
% varG = 1; % variance of AWGN at receiver
% varH = 2; % variance of channel matrix
% qSNRdB = 5:14; % quantization SNR (dB)
% qSNR = 10.^(qSNRdB/10); % quantization SNR
% qSTEP = sqrt(6*varH./qSNR); % quantization step
% %
% % Algorithms
% P0 = zeros(4,10); % minimum required transmit power
% N = 300; % number of iterations
% for a = 1:N
%     nT = 4; % number of transmit antennas
%     nR = 4; % number of receive antennas
%     H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%     for b = 1:10
%         estH = qSTEP(b).*round(H./qSTEP(b)); % estimated CSI matrix
%         [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)
%         %    
%         % Robust with all eigenmodes
%         cvx_begin
%             variable t;
%             variable delta(nR,nT) complex;
%             dual variables gamma{nT}; % dual variables
%             minimize (t);
%             max(abs(real(delta(:)))) <= qSTEP(b)/2;
%             max(abs(imag(delta(:)))) <= qSTEP(b)/2;
%             for c = 1:nT
%                 gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
%             end
%         cvx_end
%         P0(1,b) = P0(1,b)+SNR0*varG/real(cell2mat(gamma)'*diag(U'*(estH+delta)'*(estH+delta)*U));
%         %    
%         % Non-robust with maximum eigenmode
%         cvx_begin
%             variable delta(nR,nT) complex;
%             minimize (norm((estH+delta)*U(:,nT), 'fro'));
%             max(abs(real(delta(:)))) <= qSTEP(b)/2;
%             max(abs(imag(delta(:)))) <= qSTEP(b)/2;
%         cvx_end
%         P0(2,b) = P0(2,b)+SNR0*varG/real(cvx_optval^2);
%         %    
%         % OSTBC with equal power over all eigenmodes
%         cvx_begin
%             variable delta(nR,nT) complex;
%             f = 0;
%             for c = 1:nT
%                 f = f+quad_form((estH+delta)*U(:,c), diag(ones(nR,1)))/nT;
%             end
%             minimize (f);
%             max(abs(real(delta(:)))) <= qSTEP(b)/2;
%             max(abs(imag(delta(:)))) <= qSTEP(b)/2;
%         cvx_end
%         P0(3,b) = P0(3,b)+SNR0*varG/real(cvx_optval);
%     end
%     %    
%     % Robust for nT = 6 and nR = 6
%     nT = 6; % number of transmit antennas
%     nR = 6; % number of receive antennas
%     H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%     for b = 1:10
%         estH = qSTEP(b).*round(H./qSTEP(b)); % estimated CSI matrix
%         [U, ~] = eig(estH'*estH); % normalized eigenvectors and eigenvalues (ascending order)   
%         cvx_begin
%             variable t;
%             variable delta(nR,nT) complex;
%             dual variables gamma{nT}; % dual variables
%             minimize (t);
%             max(abs(real(delta(:)))) <= qSTEP(b)/2;
%             max(abs(imag(delta(:)))) <= qSTEP(b)/2;
%             for c = 1:nT
%                 gamma{c} : t >= quad_form((estH+delta)*U(:,c), diag(ones(nR,1)));
%             end
%         cvx_end
%         P0(4,b) = P0(4,b)+SNR0*varG/real(cell2mat(gamma)'*diag(U'*(estH+delta)'*(estH+delta)*U));
%     end
% end
% P0 = 10*log10(P0/N); % minimum required transmit power (dB)
% save('data/minPowerFDD.mat', 'qSNRdB', 'P0'); % save data

% Fig.7 in Paper [1]
load('data/minPowerFDD');
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 16);
hold on;
grid on;
plot(qSNRdB, P0(1,:), '-o', qSNRdB, P0(2,:), '--o', qSNRdB, P0(3,:), '-.o', qSNRdB, P0(4,:), '-*', 'linewidth', 2);
legend('Robust: n_T=4, n_R=4', 'Non-Robust: n_T=4, n_R=4', 'OSTBC: n_T=4, n_R=4', 'Robust: n_T=6, n_R=6');
xlabel('Quantization SNR (dB)');
ylabel('Mean Value of Minimum Transmit Power (dB)');
title('Fig. 7 in Paper [1]');

% Code: minPowerFDD.m


%% $\quad$ 5-4. Probability of Service Provision (Gaussian Noise with Spherical Uncertainty Regions)
%
%   close all;
% clear all;
% %
% % Parameters
% varH = 2; % variance of channel matrix
% noQoS = [0.00001 0.0001 0.001 0.01 0.05 0.1 0.2 0.4 0.6 0.8 1.0]; % probability of no QoS
% qosP = 1 - noQoS; % probability of providing required QoS
% estSNRdB = [9 7 6 5 3]; % estimation SNR (dB)
% estSNR = 10.^(estSNRdB/10); % estimation SNR
% varE = varH./estSNR; % variance of error matrix
% %
% % Algorithms
% N = 1000; % number of iterations
% pSev = N*ones(6,11); % number of service provision
% for a = 1:N
%     nT = 4; % number of transmit antennas
%     nR = 4; % number of receive antennas
%     H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%     E = randn(nR,nT) + i*randn(nR,nT); % estimation error matrix without specifying variance
%     for b = 1:5
%         estH = (H+E*sqrt(varE(b)/2))*varH/(varH+varE(b)); % estimated CSI matrix
%         r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR(b))); % radius of spherical uncertainty region
%         c = (norm(estH(:)) <= r);
%         pSev(b,c) = pSev(b,c)-1; % actual CSI matrix H=0 belongs to the uncertainty region
%     end
%     %    
%     nT = 2; % number of transmit antennas
%     nR = 2; % number of receive antennas
%     H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%     E = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varE(1)/2); % estimation error matrix
%     estH = (H+E)*varH/(varH+varE(1)); % estimated CSI matrix
%     r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR(1))); % radius of spherical uncertainty region
%     b = (norm(estH(:)) <= r);
%     pSev(6,b) = pSev(6,b)-1; % actual CSI matrix H=0 belongs to the uncertainty region
% end
% pSev = pSev/N; % probability of service provision
% save('data/probServiceTDD.mat', 'noQoS', 'pSev'); % save data

% Fig.8 in Paper [1]
load('data/probServiceTDD');
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
hold on;
grid on;
plot(noQoS, pSev(1,:), '-o', noQoS, pSev(2,:), '-+', noQoS, pSev(3,:), '-*', 'linewidth', 2);
plot(noQoS, pSev(4,:), '-s', noQoS, pSev(5,:), '-v', noQoS, pSev(6,:), '--', 'linewidth', 2);
set(gca, 'XScale', 'log', 'Fontsize', 16);
legend('n_T=4, n_R=4, SNR_{est}=9dB', 'n_T=4, n_R=4, SNR_{est}=7dB', ...
    'n_T=4, n_R=4, SNR_{est}=6dB', 'n_T=4, n_R=4, SNR_{est}=5dB', ...
    'n_T=4, n_R=4, SNR_{est}=3dB', 'n_T=2, n_R=2, SNR_{est}=9dB', 'location', 'southeast');
xlabel('1 - P_{in} (Probability of no QoS)');
ylabel('Probability of Service Provision');
title('Fig. 8 in Paper [1]');

% Code: probServiceTDD.m


%% $\quad$ 5-5. Probability of Service Provision (Quantization Error with Cubic Uncertainty Regions)
%
%   close all;
% clear all;
% %
% % Parameters
% varH = 2; % variance of channel matrix
% qSNRdB = -5:8; % quantization SNR (dB)
% qSNR = 10.^(qSNRdB/10); % quantization SNR
% qSTEP = sqrt(6*varH./qSNR); % quantization step
% nTR = [8 6 4 3 2 1]; % number of transmit and receive antennas
% %
% % Algorithms
% N = 1000; % number of iterations
% pSev = N*ones(6,14); % number of service provision
% for a = 1:N
%     for b = 1:6
%         nT = nTR(b); % number of transmit antennas
%         nR = nTR(b); % number of receive antennas
%         H = (randn(nR,nT) + i*randn(nR,nT))*sqrt(varH/2); % actual CSI matrix
%         for c = 1:14
%             estH = qSTEP(c)*round(H./qSTEP(c)); % estimated CSI matrix
%             if max(abs(real(estH(:)))) <= qSTEP(c)/2 && max(abs(imag(estH(:)))) <= qSTEP(c)/2
%                 pSev(b,c) = pSev(b,c)-1; % actual CSI matrix H=0 belongs to the uncertainty region
%             end
%         end
%     end
% end
% pSev = pSev/N; % probability of service provision
% save('data/probServiceFDD.mat', 'qSNRdB', 'pSev'); % save data

% Fig.9 in Paper [1]
load('data/probServiceFDD');
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 16);
hold on;
grid on;
plot(qSNRdB, pSev(1,:), '-o', qSNRdB, pSev(2,:), '-*', qSNRdB, pSev(3,:), '--o', 'linewidth', 2);
plot(qSNRdB, pSev(4,:), '--*', qSNRdB, pSev(5,:), '-.o', qSNRdB, pSev(6,:), '-.*', 'linewidth', 2);
xlim([-5 8]);
legend('n_T=8, n_R=8', 'n_T=6, n_R=6', 'n_T=4, n_R=4', ...
    'n_T=3, n_R=3', 'n_T=2, n_R=2', 'n_T=1, n_R=1', 'location', 'southeast');
xlabel('Quantization SNR (dB)');
ylabel('Probability of Service Provision');
title('Fig. 9 in Paper [1]');

% Code: probServiceFDD.m


%%
close all;
clear all;



%% *6. Concluding Remarks*
%
% In this project, I have reproduced the main results of the assigned
% paper with the CVX package, while the authors implement the proposed
% algorithm with the function $\mathrm{fmincon}$ of the optimization
% toolbox in MATLAB. For the assigned paper, the investigated problem is
% interesting and important, and the solution to the optimization problem
% of the proposed algorithm is very impressive. 
%
% According to the reproduced results, we can make the same key
% observations as in paper [1]. Nevertheless, here are several aspects to
% clarify:
%
% $\quad$ (a). There are two typos which effect the simulation results: Eq.
% (17) in paper [1] should be $\mathbf{m}_{h|\tilde{h}} = \mathbf{m}_h +
% \mathbf{C}_h (\mathbf{C}_h + \mathbf{C}_e)^{-1} (\tilde{h} -
% \mathbf{m}_h)$. Beside, the service provision probability in Fig. 8 and
% Fig. 9 should be $\mathrm{Pr} (\mathrm{\hat{H}} + \Delta \neq 0, ~\forall
% \Delta \in \mathcal{R})$.
%
% $\quad$ (b). For the simulation of Fig. 6 in paper [1], the result is
% averaged over only two simulation rounds. The main reason is that it is
% of high probability that there is no service provision according to Fig.
% 8, which means the minimum transmit power is infinite. One way to deal
% with this problem is to add the condition: $\tilde{\mathcal{R}} = \{
% \mathrm{H}: ||\mathrm{H}||_F^2 \ge \rho \}$ , where $\rho$ is a positive
% real value. However, this will make uncertainty region nonconvex.
% Although the result is obtained with two iterations, we can make the same
% key observations as the authors.
%
% $\quad$ (c). As the authors do not specify the simulation settings in
% this paper, the parameters used in this project are very likely to be
% different from the original paper. Thus, the results or settings of this
% project are slightly different from the original paper (e.g., Fig. 7 and
% Fig. 8).
%
% $\quad$ (d). There are seven figures in the paper, and I reproduce five
% of them in this project. Actually, the code of other two figures can be
% easily implemented by adapting existing implementations, however, it would
% take a long time to run the code. Besides, all the core algorithms and
% reference methods have been implemented through the five experiments.
%
%
%
%% *7. Reference*
% [1] A. Pascual-Iserte, et al., "A Robust Maximin Approach for MIMO
% Communications With Imperfect Channel State Information Based on Convex
% Optimization," IEEE Trans. Signal Processing, vol. 54, no. 1, pp.
% 346-360, 2006.
%
%