RIS Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RIS Phase Hopping for Ultra-Reliable Communications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Author:_ Karl-Ludwig Besser, Technische Universität Braunschweig\n",
    "\n",
    "This notebook is part of the papers \"Reconfigurable Intelligent Surface Phase Hopping for Ultra-Reliable Communications\" (K.-L. Besser and E. Jorswieck, IEEE Transactions on Wireless Communications, vol. 21, no. 11, pp. 9082–9095, Nov 2022. [doi:10.1109/TWC.2022.3172760](https://doi.org/10.1109/TWC.2022.3172760), [arXiv:2107.11852](https://arxiv.org/abs/2107.11852)) and \"Artificial Fast Fading from Reconfigurable Surfaces Enables Ultra-Reliable Communications\" (SPAWC 2021 [doi:10.1109/SPAWC51858.2021.9593269](https://doi.org/10.1109/SPAWC51858.2021.9593269)).  \n",
    "If you use any of this work, please cite the above paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> If you are not familiar with Jupyter notebooks: The easiest way to use this notebook interactively, is to hit `Kernel --> Restart & Run All` in the menu. This will execute all cells and enable the interactive elements of the plots.  \n",
    "> Alternatively, you can execute the cells one by one using Shift+Return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "%matplotlib widget\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, interact_manual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from phases import gains_constant_phase, rvs_channel_phases, rvs_ris_phases, rvs_ris_phases_quant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constant RIS Phases\n",
    "\n",
    "First, we consider constant RIS phases _without_ phase hopping.\n",
    "\n",
    "It can be seen that the $\\varepsilon$-outage capacity is close to zero for $\\varepsilon$ close to zero (for all $N$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def constant_ris_phases(num_samples=50000, num_elements=[5, 10, 20, 50]):\n",
    "    fig, axs = plt.subplots()\n",
    "    def update(show_approx=True):\n",
    "        axs.clear()\n",
    "        axs.set_xlabel(\"Rate $R$\")\n",
    "        axs.set_ylabel(\"Outage Probability $\\\\varepsilon$\")\n",
    "        for _num_elements in num_elements:\n",
    "            channel_realizations = rvs_channel_phases(_num_elements, num_samples)\n",
    "            const_phase = gains_constant_phase(channel_realizations)\n",
    "            capac_const_phase = np.log2(1 + const_phase)\n",
    "            _hist = np.histogram(capac_const_phase, bins=100)\n",
    "            _r_ax = np.linspace(min(capac_const_phase)*.9, max(capac_const_phase)*1.1, 1000)\n",
    "            cdf_hist = stats.rv_histogram(_hist).cdf(_r_ax)\n",
    "            axs.plot(_r_ax, cdf_hist, label=\"Empirical CDF N={:d}\".format(_num_elements))\n",
    "            if show_approx:\n",
    "                cdf_appr = 1. - np.exp(-(2**_r_ax-1)/_num_elements)  # N --> oo, for sum\n",
    "                axs.plot(_r_ax, cdf_appr, '--', label=\"Approximate N={:d}\".format(_num_elements))\n",
    "        axs.legend()\n",
    "    interact(update, show_approx=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "constant_ris_phases()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RIS Phase Hopping\n",
    "\n",
    "Next, we consider RIS phase hopping with randomly varied phases $\\theta_i$, $i=1, \\dots{}, N$.\n",
    "All $\\theta_i$ are iid with a uniform distribution over $[0, 2\\pi]$.\n",
    "\n",
    "It can be seen that the outage probability is a step function when phase hopping is employed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def n_element_fast_fading(num_samples_slow=1000, num_samples_fast=5000):\n",
    "    fig, axs = plt.subplots()\n",
    "    def update(num_elements=3):\n",
    "        axs.clear()\n",
    "        channel_realizations = rvs_channel_phases(num_elements, num_samples_slow)\n",
    "        channel_realizations = np.tile(channel_realizations, (num_samples_fast, 1, 1))\n",
    "        ris_phases = rvs_ris_phases(num_elements, num_samples_slow, num_samples_fast, copula=\"indep\")\n",
    "        channel_realizations = channel_realizations + ris_phases\n",
    "        const_phase = gains_constant_phase(channel_realizations)\n",
    "        capac_const_phase = np.log2(1 + const_phase)\n",
    "        expect_capac = np.mean(capac_const_phase, axis=0)\n",
    "        _hist = np.histogram(expect_capac, bins=100)\n",
    "        _r_ax = np.linspace(0, 3, 1000)\n",
    "        cdf_hist = stats.rv_histogram(_hist).cdf(_r_ax)\n",
    "        axs.plot(_r_ax, cdf_hist)\n",
    "        axs.set_title(\"RIS Phase Hopping with N={:d} RIS Elements\".format(num_elements))\n",
    "        axs.set_xlabel(\"Rate $R$\")\n",
    "        axs.set_ylabel(\"Outage Probability $\\\\varepsilon$\")\n",
    "    interact_manual(update, num_elements=(2, 10, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_element_fast_fading()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantized Phases\n",
    "\n",
    "Previously, continuous RIS phases over the set $[0, 2\\pi]$ have been considered.\n",
    "In the following, we assume that only a discrete set of possible phase values is available, i.e.,\n",
    "$$\\theta_i\\in\\mathcal{Q}=\\left\\{k\\frac{2\\pi}{K}\\;\\middle|\\; k=0,\\dots{},K-1\\right\\}.$$\n",
    "\n",
    "The next plot shows an example with $N=2$ and $K=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantized_two_phases_two_elements(num_samples_slow=1000, num_samples_fast=5000):\n",
    "    fig, axs = plt.subplots()\n",
    "    num_elements = 2\n",
    "    dependency = \"indep\"\n",
    "    channel_realizations = rvs_channel_phases(num_elements, num_samples_slow)\n",
    "    channel_realizations = np.tile(channel_realizations, (num_samples_fast, 1, 1))\n",
    "    ris_phases = rvs_ris_phases_quant(num_elements, num_samples_slow, num_samples_fast,\n",
    "                                      copula=dependency, K=2)\n",
    "    total_phases = channel_realizations + ris_phases\n",
    "    const_phase = gains_constant_phase(total_phases)\n",
    "    capac_const_phase = np.log2(1 + const_phase)\n",
    "    expect_capac = np.mean(capac_const_phase, axis=0)\n",
    "    _hist = np.histogram(expect_capac, bins=100)\n",
    "    _r_ax = np.linspace(0, 3, 1000)\n",
    "    cdf_hist = stats.rv_histogram(_hist).cdf(_r_ax)\n",
    "    zoc = 0.5*np.log2(5)\n",
    "    axs.plot(_r_ax, cdf_hist, label=\"ECDF\")\n",
    "    axs.vlines(zoc, 0, 1, 'r', label=\"ZOC -- Exact\")\n",
    "    axs.legend()\n",
    "    axs.set_title(\"Artificial Fast Fading with N={:d} RIS Elements\\nQuantized Phases with 2 Quantization Steps\".format(num_elements))\n",
    "    axs.set_xlabel(\"Rate $R$\")\n",
    "    axs.set_ylabel(\"Outage Probability $\\\\varepsilon$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quantized_two_phases_two_elements(num_samples_fast=100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we visualize a scheme that allows for a positive ZOC using an RIS with $N=2$ elements and $K=2$ possible phase values.\n",
    "Each symbol is repeated $L=2$ times.\n",
    "The overall capacity is then given as\n",
    "\\begin{equation}\n",
    "\t\t\\frac{1}{L} \\log_2\\left(1 + \\sum_{l=1}^{L}\\left|{\\sum_{i=1}^{N}\\exp\\left(\\mathrm{j} \\left(\\theta_{i,l} + \\varphi_i\\right)\\right)}\\right|^2\\right)\\,.\n",
    "\t\\end{equation}\n",
    "\n",
    "The phases of the RIS elements are set to $\\theta_{1,1}=\\theta_{1,2}=0$, and $\\theta_{2,1}=0$ and $\\theta_{2,2}=\\pi$.\n",
    "By this, we can achieve a constant value of the resulting SNR as\n",
    "$$\\sum_{l=1}^{2}\\left|{\\sum_{i=1}^{2}\\exp\\left(\\mathrm{j} \\left(\\theta_{i,l} + \\varphi_i\\right)\\right)}\\right|^2 = 4,$$\n",
    "independent of the channel realizations $\\varphi_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def positive_zoc_repetition():\n",
    "    fig, (axs, axs2) = plt.subplots(1,2)\n",
    "    axs.set_xlim([-2, 2])\n",
    "    axs.set_ylim([-2, 2])\n",
    "    axs2.set_xlim([0, 1])\n",
    "    axs2.set_ylim([0, 5])\n",
    "    axs2.set_ylabel(\"Absolute Value Squared\")\n",
    "    plot1 = axs.plot([0, 1], [0, 1], '--')[0]\n",
    "    plot21 = axs.plot([0, 1], [0, 1], '--')[0]\n",
    "    plot22 = axs.plot([0, 1], [0, 1], '--')[0]\n",
    "    plot_abs1 = axs.plot([0, 1], [0, 1], 'o-')[0]\n",
    "    plot_abs2 = axs.plot([0, 1], [0, 1], 'o-')[0]\n",
    "    plot2_abs1 = axs2.plot([0, 1], [0, 1], '-', color=plot_abs1.get_color(), lw=4)[0]\n",
    "    plot2_abs2 = axs2.plot([0, 1], [0, 1], '-', color=plot_abs2.get_color(), lw=4)[0]\n",
    "    def update(phi1=np.pi/4, phi2=0):\n",
    "        z1 = np.exp(1j*phi1)\n",
    "        z21 = z1 + np.exp(1j*phi2)\n",
    "        z22 = z1 + np.exp(1j*(phi2+np.pi))\n",
    "        plot1.set_data([0, np.real(z1)], [0, np.imag(z1)])\n",
    "        plot21.set_data([np.real(z1), np.real(z21)], [np.imag(z1), np.imag(z21)])\n",
    "        plot22.set_data([np.real(z1), np.real(z22)], [np.imag(z1), np.imag(z22)])\n",
    "        \n",
    "        plot_abs1.set_data([0, np.real(z21)], [0, np.imag(z21)])\n",
    "        plot_abs2.set_data([0, np.real(z22)], [0, np.imag(z22)])\n",
    "        \n",
    "        plot2_abs1.set_data([0.5, 0.5], [0, np.abs(z21)**2])\n",
    "        plot2_abs2.set_data([0.5, 0.5], [np.abs(z21)**2, np.abs(z21)**2 + np.abs(z22)**2])\n",
    "    interact(update, phi1=(0, 2*np.pi, .01), phi2=(0, 2*np.pi, .01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "positive_zoc_repetition()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}