{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Power Assignment in a Wireless Communication System\n",
    "\n",
    "by Robert Gowers, Roger Hill, Sami Al-Izzi, Timothy Pollington and Keith Briggs\n",
    "\n",
    "from Boyd and Vandenberghe, Convex Optimization, exercise 4.20 page 196\n",
    "\n",
    "Convex optimization can be used to maximise the minimum signal to inteference plus noise ratio (SINR) of a wireless communication system.  Consider a system with $n$ transmitters, each with power $p_j \\geq 0$, transmitting to $n$ receivers.  Let $G_{ij} \\geq 0$ denote the path gain from transmitter $j$ to receiver $i$.  These path gains form the matrix $G \\in \\mathbb{R}^{n \\times n}$.\n",
    "\n",
    "Each receiver is assigned to a transmitter such that the signal power at receiver $i$, $S_i = G_{ii}p_i$ and the interefence power at receiver $i$ is $I_i = \\sum_{k\\neq i} G_{ik}p_k$.  Given a noise power $\\sigma_i$ at each receiver, the SINR at receiver $i$, $\\gamma_i = \\frac{S_i}{I_i + \\sigma_i}$.\n",
    "\n",
    "The objective is to maximise the minimum SINR of the system under certain power constraints.  These constraints are:\n",
    "\n",
    "i - Each transmitter power $p_j \\leq P_j^{\\text{max}}$\n",
    "\n",
    "ii - If the transmitters are partitioned into $m$ nonoverlapping groups, $K_1, ..., K_m$, which share a common power supply with total power $P_l^{\\text{gp}}$: $\\sum_{k\\in K_l}p_k \\leq P_l^{\\text{gp}}$.\n",
    "\n",
    "iii - There is a maximum power that each receiver can receive $P_i^{\\text{rc}}$, $\\sum_{k=1}^{n}G_{ik}p_k \\leq P_i^{\\text{rc}}$.\n",
    "\n",
    "The objective function can be rewritten as:\n",
    "\n",
    "minimise $\\max_{i=1,...,n}\\frac{I_i + \\sigma_i}{S_i}$\n",
    "\n",
    "However, since this is a quasiconvex objective function we cannot solve it directly using CVXPY.  Instead we must use a bisection method.  First we take the step of rewriting the objective, $\\alpha = \\gamma^{-1} \\geq 0$, as a constraint:\n",
    "\n",
    "$I_i+\\sigma_i \\leq S_i\\alpha$\n",
    "\n",
    "Then we choose initial lower and upper bounds $L_0$ and $U_0$ for $\\alpha$, which should be chosen such that $L < \\alpha^* < U$, where $\\alpha^*$ is the optimal value of $\\alpha$.  Starting with an initial value $\\alpha_0 = \\frac{1}{2}(L_0+U_0)$, feasibility is checked for $\\alpha_0$ by using an arbitrary objective function.  The new upper and lower bounds are determined from the feasibility:\n",
    "\n",
    "If $\\alpha_0$ is feasible then $L_1 = L_0$, $U_1 = \\alpha_0$ and $\\alpha_1 = \\frac{1}{2}(L_1+U_1)$.\n",
    "\n",
    "If $\\alpha_0$ is infeasible then $L_1 = \\alpha_1$, $U_1 = U_0$ and $\\alpha_1 = \\frac{1}{2}(L_1+U_1)$.\n",
    "\n",
    "This bisection process is repeated until $U_N - L_N < \\epsilon$, where $\\epsilon$ is the desired tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "# @author: R. Gowers, S. Al-Izzi, T. Pollington, R. Hill & K. Briggs\n",
    "\n",
    "import cvxpy as cp\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def maxmin_sinr(G, P_max, P_received, sigma, Group, Group_max, epsilon = 0.001):\n",
    "    # find n and m from the size of the path gain matrix\n",
    "    n, m = np.shape(G)\n",
    "    \n",
    "    # Checks sizes of inputs\n",
    "    if m != np.size(P_max):\n",
    "        print('Error: P_max dimensions do not match gain matrix dimensions\\n')\n",
    "        return 'Error: P_max dimensions do not match gain matrix dimensions\\n', np.nan, np.nan, np.nan\n",
    "    \n",
    "    if n != np.size(P_received):\n",
    "        print('Error: P_received dimensions do not match gain matrix dimensions\\n')\n",
    "        return 'Error: P_received dimensions do not match gain matrix dimensions', np.nan, np.nan, np.nan\n",
    "    \n",
    "    if n != np.size(sigma):\n",
    "        print('Error: σ dimensions do not match gain matrix dimensions\\n')\n",
    "        return 'Error: σ dimensions do not match gain matrix dimensions', np.nan, np.nan, np.nan\n",
    "\n",
    "    #I = np.zeros((n,m))\n",
    "    #S = np.zeros((n,m))\n",
    "\n",
    "    delta = np.identity(n)\n",
    "    S = G*delta # signal power matrix\n",
    "    I = G-S # interference power matrix\n",
    "\n",
    "    # group matrix: number of groups by number of transmitters\n",
    "    num_groups = int(np.size(Group,0))\n",
    "\n",
    "    if num_groups != np.size(Group_max):\n",
    "        print('Error: Number of groups from Group matrix does not match dimensions of Group_max\\n')\n",
    "        return ('Error: Number of groups from Group matrix does not match dimensions of Group_max',\n",
    "                np.nan, np.nan, np.nan, np.nan)\n",
    "\n",
    "    # normalising the max power of a group so it is in the range [0,1]\n",
    "    Group_norm = Group/np.sum(Group,axis=1).reshape((num_groups,1))\n",
    "    \n",
    "    # create scalar optimisation variable p: the power of the n transmitters\n",
    "    p = cp.Variable(shape=n)\n",
    "    best = np.zeros(n)\n",
    "\n",
    "    # set upper and lower bounds for sub-level set\n",
    "    u = 1e4\n",
    "    l = 0\n",
    "\n",
    "    # alpha defines the sub-level sets of the generalised linear fractional problem\n",
    "    # in this case α is the reciprocal of the minimum SINR\n",
    "    alpha = cp.Parameter(shape=1)\n",
    "    \n",
    "    # set up the constraints for the bisection feasibility test\n",
    "    constraints = [I*p + sigma <= alpha*S*p, p <= P_max, p >= 0, G*p <= P_received, Group_norm*p <= Group_max]\n",
    "\n",
    "    # define objective function, in our case it's constant as only want to test the solution's feasibility\n",
    "    obj = cp.Minimize(alpha)\n",
    "    \n",
    "    # now check whether the solution lies between u and l\n",
    "    alpha.value = [u]\n",
    "    prob = cp.Problem(obj, constraints)\n",
    "    prob.solve()\n",
    "    \n",
    "    if prob.status != 'optimal':\n",
    "        # in this case the level set u is below the solution\n",
    "        print('No optimal solution within bounds\\n')\n",
    "        return 'Error: no optimal solution within bounds', np.nan, np.nan, np.nan\n",
    "    \n",
    "    alpha.value = [l]\n",
    "    prob = cp.Problem(obj, constraints)\n",
    "    prob.solve()\n",
    "\n",
    "    if prob.status == 'optimal':\n",
    "        # in this case the level set l is below the solution\n",
    "        print('No optimal solution within bounds\\n')\n",
    "        return 'Error: no optimal solution within bounds', np.nan, np.nan, np.nan\n",
    "    \n",
    "    # Bisection algortithm starts\n",
    "    maxLoop = int(1e7)\n",
    "    for i in range(1,maxLoop):\n",
    "        # First check that u is in the feasible domain and l is not, loop finishes here if this is not the case\n",
    "        # set α as the midpoint of the interval\n",
    "        alpha.value = np.atleast_1d((u + l)/2.0)\n",
    "\n",
    "        # test the size of the interval against the specified tolerance\n",
    "        if u-l <= epsilon:\n",
    "            break\n",
    "        \n",
    "        # form and solve problem\n",
    "        prob = cp.Problem(obj, constraints)\n",
    "        prob.solve()\n",
    "\n",
    "        # If the problem is feasible u -> α, if not l -> α, best takes the last feasible value as the optimal one as\n",
    "        # when the tolerance is reached the new α may be out of bounds\n",
    "        if prob.status == 'optimal':\n",
    "            u = alpha.value\n",
    "            best = p.value\n",
    "        else:\n",
    "            l = alpha.value\n",
    "            \n",
    "        # final condition to check that the interval has converged to order ε, i.e. the range of the optimal sublevel set is <=ε\n",
    "        if u - l > epsilon and i == (maxLoop-1):\n",
    "            print(\"Solution not converged to order epsilon\")\n",
    "    \n",
    "    return l, u, float(alpha.value), best\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "As a simple example, we will consider a case with $n=5$, where $G_{ij} = 0.6$ if $i=j$ and $0.1$ otherwise. \n",
    "\n",
    "$P_j^{\\text{max}} = 1$ for all transmitters and the transmitters are split into two groups, each with $P_l^{\\text{gp}} = 1.8$.  The first group contains transmitters 1 & 2, while the second group contains 3,4 & 5.\n",
    "\n",
    "For all receivers $P_i^{\\text{rc}} = 4$ and $\\sigma_i = 0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum SINR=1.148\n",
      "Power=[0.8 0.8 0.8 0.8 0.8]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(precision=3)\n",
    "\n",
    "# in this case we will use a gain matrix with a signal weight of 0.6 and interference weight of 0.1\n",
    "G = np.array([[0.6,0.1,0.1,0.1,0.1],\n",
    "              [0.1,0.6,0.1,0.1,0.1],\n",
    "              [0.1,0.1,0.6,0.1,0.1],\n",
    "              [0.1,0.1,0.1,0.6,0.1],\n",
    "              [0.1,0.1,0.1,0.1,0.6]])\n",
    "\n",
    "# in this case m=n, but this generalises if we want n receivers and m transmitters\n",
    "n, m = np.shape(G)\n",
    "\n",
    "# set maximum power of each transmitter and receiver saturation level\n",
    "P_max = np.array([1.]*n)\n",
    "\n",
    "# normalised received power, total possible would be all power from all transmitters so 1/n\n",
    "P_received = np.array([4.,4.,4.,4.,4.])/n\n",
    "\n",
    "# set noise level\n",
    "sigma = np.array([0.1,0.1,0.1,0.1,0.1])\n",
    "\n",
    "# group matrix: number of groups by number of transmitters\n",
    "Group = np.array([[1.,1.,0,0,0],[0,0,1.,1.,1.]])\n",
    "\n",
    "# max normalised power for groups, number of groups by 1\n",
    "Group_max = np.array([1.8,1.8])\n",
    "\n",
    "# now run the optimisation problem\n",
    "l, u, alpha, best = maxmin_sinr(G, P_max, P_received, sigma, Group, Group_max)\n",
    "\n",
    "print('Minimum SINR={:.4g}'.format(1/alpha))\n",
    "print('Power={}'.format(best))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
