%% A Robust Maximin Approach for MIMO Communications With Imperfect Channel State Information Based on Convex Optimization
%
% By
%
% Shahed Ali Khan 
%
% Srivatsava Reddy Dubbudu 
%
% Vishnu Sai Adusumalli 
% 
%% Introduction to MIMO systems
% Multiple-input multiple-output, or MIMO, is a radio communications technology 
% or RF technology that is being mentioned and used in many new technologies these days.
% MIMO channels have become a popular means to increase the spectral
% efficiency and the quality of wireless communications by means of spatial diversity at both sides of
% the link. The design of the communication system depends on
% the quantity and the quality of the channel state information
% (CSI) available at both the transmitter and the receiver. In case
% that the transmitter does not have any information about the
% channel, then space–time coding techniques should be applied. The optimum exploitation of the benefits provided by
% MIMO systems can be achieved when a perfect CSI is available
% at both the transmitter and the receiver. In such a case, the
% optimum solution consists in the joint design of the transmitter
% and the receiver taking into account the channel response and
% the noise plus interference profile.
%% System Model and Problem Statement
% Consider the transmission through a multi-antenna flat 
% fading wireless MIMO channel, the $\mathrm{\hat{H}}$ represents channel gain.
% In case that the channel knowledge is imperfect then the maximum eigenmode of 
% $\mathrm{\hat{H}}^H \mathrm{\hat{H}}$ is used.
% A transmitter is designed
% according to different degrees of channel knowledge, showing
% that in some situations, the capacity achieving solution may
% need to use multiple beamforming. 
% The design of the transmitter will be based on a linear processing
% or beamforming scheme, whereas at the receiver, a maximum
% likelihood (ML) detector will be employed assuming a
% perfect channel knowledge. 
% The architecture for the
% robust transmitter is composed of an OSTBC block, whose outputs
% are transmitted through all the eigenmodes of $\mathrm{\hat{H}}^H
% \mathrm{\hat{H}}$ using
% an adequate power distribution among them, as opposed to the
% nonrobust design, which uses only the maximum eigenmode.
% 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) $$
%% The MAXMIN problem
% The direct way to solve the maximin problem is to obtain the
% minimization of  $f$ analytically and then solve the outer maximization,
% either numerically or analytically. Such an approach,
% however, is difficult because it is not clear what is the minimizing  $\Delta$
% for a given $p$ in closed-form. 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) $$
%% Convex Uncertainity Regions
% The definition of the uncertainty region $\mathcal{R}$ may impact importantly
% on the system performance. The size and the shape of this
% region should take into account the quality of the channel estimate
% and the imperfections that generate the error, linking the
% mathematical optimization problem and the physical phenomenon
% producing the error.
% In the following, two sources of errors are identified and three
% different uncertainty regions, jointly with their sizes, are described.
%
% 1. Estimation Gaussian Noise - 
% 
% A usual error in the channel estimate comes from the
% Gaussian noise, especially in TDD systems, where the transmitter
% can estimate the channel using the signals received in the
% reverse link, and use it as an estimate in the forward link, due
% to the channel reciprocity principle.
% 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.
%
% 2. Quantization Errors - 
%
% In frequency division duplexing (FDD) systems, the estimate
% of the channel at the transmitter has to be obtained through a
% feedback channel from the receiver to the transmitter. Since this
% feedback channel has a finite capacity, the channel response has
% to be quantized introducing an error in the CSI available at the
% transmitter. 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) $$
%
% 3. Combined Estimation and Quantization Errors - 
%
% In a realistic scenario, the two effects considered previously,
% i.e., the Gaussian noise from the estimation process and the
% quantization errors, are expected to be combined. This can be
% modeled mathematically by defining an appropriate uncertainty
% region for the error, which can be expressed 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) $$
%
% 4. Other Uncertainity Regions - 
%
% In addition to the previous uncertainty regions, there are
% many other possibilities and feedback strategies, whose error
% and imperfection models also lead to convex uncertainty regions.
% The regions described previously make sense when
% the whole channel matrix $\mathrm{H}$ is estimated or each element of  $\mathrm{H}$
% is independently quantized.
%
% 5. A Closed-Form Solution for Spherical Uncertainty Regions - 
%
% The general convex problem presented in (8) can be extremely
% simplified for the case of spherical uncertainty regions.
% A maximin approach has to be taken to find
% a solution to the problem of the maximization of the SNR assuming spherical uncertainty regions with a prefixed size,
% where no discussion is given on how to calculate the value
% of the radius of the region.
%% Making a facsimile of the simulation using CVX instead of fmincon function of the optimization toolbox in MATLAB
%%
% 1. Mean value of the power distribution for different sizes of the
% uncertainty region
%
close all;
clear all;
varE = 1; % Var of the error matrix 
varH = 2; % Var of the channel matrix 
nR = 6; % No. of receiving antennas
nT = 4; % No. of transmiting antennas
g = 0:0.1:0.9;
g = [g 0.95 0.99];
N = 3; % No. of iterations
pDist = zeros(nT,12); % Normalized power allocation
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); % Est error matrix
  estH = (H+E)*varH/(varH+varE); % Est 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
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 18);
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('Figure 1')
%%
%
% 2. Mean value of the minimum required transmit power in a TDD system
% (Gaussian noise, spherical uncertainty regions) according to a target SNR equal
% to $\mathrm{SNR}_0$ = 10 dB
%
close all;
clear all;
SNR0 = 10; % Target Signal to Noise ratio
varH = 2; % Var of channel matrix
varG = 1; % Var of AWGN at Rx
estSNRdB = 8:14; % Estimated SNR (dB)
estSNR = 10.^(estSNRdB/10); % Estimated SNR
nT = 2; % No. of Tx antennas
nR = 2; % No. of Rx antennas
P0 = zeros(4,7); % Min power
N = 2; % No. 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 with n variance
  qosP = 0.85; % Prob 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); % Var of error matrix
      estH = (H+E*sqrt(varE/2))*varH/(varH+varE); % Est CSI matrix
      [U, ~] = eig(estH'*estH); % Normalized eigenvectors and eigenvalues
      %
      % 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; % Prob of providing QoS
  r = sqrt(chi2inv(qosP, 2*nR*nT)*varH./(1+estSNR)); % Radius of spherical uncertainty region
  for b = 1:7
      varE = varH/estSNR(b); % Var of error matrix
      estH = (H+E*sqrt(varE/2))*varH/(varH+varE); % Est 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); % Min required transmit power (dB)
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 18);
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('Figure 2');
%%
%
% 3. Mean value of the minimum required transmit power in a FDD system
% (quantization errors, cubic uncertainty regions) according to a target SNR equal
% to $\mathrm{SNR}_0$ = 10 dB
%
close all;
clear all;
SNR0 = 10; % Target Signal to Noise ratio
varH = 2; % Var of channel matrix
varG = 1; % Var of AWGN at Rx
qSNRdB = 5:14; % Quantized SNR (dB)
qSNR = 10.^(qSNRdB/10); % Quantized SNR
qSTEP = sqrt(6*varH./qSNR); % Quantization step
P0 = zeros(4,10); % Min required transmit power
N = 15; % No. of iterations
for a = 1:N
  nR = 4; % No. of Rx antennas
  nT = 4; % No. of Tx 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)); % Est CSI matrix
      [U, ~] = eig(estH'*estH); % Normalized eigenvectors and eigenvalues
      %
      % 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 max 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
  nR = 6; % No. of Rx antennas
  nT = 6; % No. of Tx 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)); % Est CSI matrix
      [U, ~] = eig(estH'*estH); % Normalized eigenvectors and eigenvalues
      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); % Min required transmit power (dB)
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('Figure 3');

%%
%
% 4. Probability of service provision versus quantization SNR, assuming
% quantization errors and cubic uncertainty regions
%
close all;
clear all;
varH = 2; % Var of channel matrix
qSNRdB = -5:8; % Quantized SNR (dB)
qSNR = 10.^(qSNRdB/10); % Quantized SNR
qSTEP = sqrt(6*varH./qSNR); % Quantized step
nTR = [8 6 4 3 2 1]; % No. of Tx and Rx antennas
N = 1000; % No. of iterations
pSev = N*ones(6,14); % No. of service provision
for a = 1:N
  for b = 1:6
      nT = nTR(b); % No. of transmit antennas
      nR = nTR(b); % No. 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)); % Est CSI matrix
          if max(abs(real(estH(:)))) <= qSTEP(c)/2 && max(abs(imag(estH(:)))) <= qSTEP(c)/2
              pSev(b,c) = pSev(b,c)-1;
          end
      end
  end
end
pSev = pSev/N; % Prob of service provision
fig = figure;
set(fig, 'Units', 'Pixels', 'Position', [100 600 800 500]);
set(gca, 'Fontsize', 18);
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('Figure 4')
%
%% Conclusion
% A robust design strategy has been proposed for a wireless multi-antenna 
% MIMO channel. The transmitter is de- signed based on a channel estimate that 
% may have errors from different origins, such as the Gaussian noise from the
% estimation process, or the errors from the quantization of the channel estimate,
% among others. These errors have been taken into account explicitly in the design,
% obtaining a robust solution less sensitive to these errors. The transmitter 
% architecture that has been proposed is based on the 
% concatenation of an OSTBC block, a power allocation,and a set of beamformers, 
% each one connected to one of the outputs of 
% the OSTBC and corresponding to an eigenmode of the MIMO channel estimate.
%
% The mathematical optimization problem corresponding to the maximin 
% robust power allocation has been transformed into a simple convex optimization 
% problem that can be solved efficiently by existing software packages.
% The authors have used fmincon of the optimization toolbox in MATLAB for 
% simulation but we have used CVX
% package for simulating 4 out of 7 figures. From the simulations results, it has been observed that 
% the gains in terms of saving in transmit power are quite important 
% when compared to the nonrobust and the pure OSTBC techniques,
% especially when the estimation and the quantization SNR are low.
%
%% References 
% [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.

