{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimal Power and Bandwidth Allocation in a Gaussian Channel\n",
    "by Robert Gowers, Roger Hill, Sami Al-Izzi, Timothy Pollington and Keith Briggs\n",
    "\n",
    "from Boyd and Vandenberghe, Convex Optimization, exercise 4.62 page 210\n",
    "\n",
    "Consider a system in which a central node transmits messages to $n$ receivers.  Each receiver channel $i \\in \\{1,...,n\\}$ has a transmit power $P_i$ and bandwidth $W_i$.  A fraction of the total power and bandwidth is allocated to each channel, such that $\\sum_{i=1}^{n}P_i = P_{tot}$ and $\\sum_{i=1}^{n}W_i = W_{tot}$.  Given some utility function of the bit rate of each channel, $u_i(R_i)$, the objective is to maximise the total utility $U = \\sum_{i=1}^{n}u_i(R_i)$.\n",
    "\n",
    "Assuming that each channel is corrupted by Gaussian white noise, the signal to noise ratio is given by $\\beta_i P_i/W_i$.  This means that the bit rate is given by:\n",
    "\n",
    "$R_i = \\alpha_i W_i \\log_2(1+\\beta_iP_i/W_i)$\n",
    "\n",
    "where $\\alpha_i$ and $\\beta_i$ are known positive constants.\n",
    "\n",
    "One of the simplest utility functions is the data rate itself, which also gives a convex objective function.\n",
    "\n",
    "The optimisation problem can be thus be formulated as:\n",
    "\n",
    "minimise $\\sum_{i=1}^{n}-\\alpha_i W_i \\log_2(1+\\beta_iP_i/W_i)$\n",
    "\n",
    "subject to $\\sum_{i=1}^{n}P_i = P_{tot} \\quad \\sum_{i=1}^{n}W_i = W_{tot} \\quad P \\succeq 0 \\quad W \\succeq 0$\n",
    "\n",
    "Although this is a convex optimisation problem, it must be rewritten in DCP form since $P_i$ and $W_i$ are variables and DCP prohibits dividing one variable by another directly.  In order to rewrite the problem in DCP format, we utilise the $\\texttt{kl_div}$ function in CVXPY, which calculates the Kullback-Leibler divergence.\n",
    "\n",
    "$\\text{kl_div}(x,y) = x\\log(x/y)-x+y$\n",
    "\n",
    "$-R_i = \\text{kl_div}(\\alpha_i W_i, \\alpha_i(W_i+\\beta_iP_i)) - \\alpha_i\\beta_iP_i$\n",
    "\n",
    "Now that the objective function is in DCP form, the problem can be solved using CVXPY."
   ]
  },
  {
   "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 numpy as np\n",
    "import cvxpy as cp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimal_power(n, a_val, b_val, P_tot=1.0, W_tot=1.0):\n",
    "    # Input parameters: α and β are constants from R_i equation\n",
    "    n = len(a_val)\n",
    "    if n != len(b_val):\n",
    "        print('alpha and beta vectors must have same length!')\n",
    "        return 'failed', np.nan, np.nan, np.nan\n",
    "    \n",
    "    P = cp.Variable(shape=n)\n",
    "    W = cp.Variable(shape=n)\n",
    "    alpha = cp.Parameter(shape=n)\n",
    "    beta = cp.Parameter(shape=n)\n",
    "    alpha.value = np.array(a_val)\n",
    "    beta.value = np.array(b_val)\n",
    "\n",
    "    # This function will be used as the objective so must be DCP; \n",
    "    # i.e. elementwise multiplication must occur inside kl_div, \n",
    "    # not outside otherwise the solver does not know if it is DCP...\n",
    "    R = cp.kl_div(cp.multiply(alpha, W),\n",
    "                  cp.multiply(alpha, W + cp.multiply(beta, P))) - \\\n",
    "                  cp.multiply(alpha, cp.multiply(beta, P))\n",
    "\n",
    "    objective = cp.Minimize(cp.sum(R))\n",
    "    constraints = [P>=0.0,\n",
    "                   W>=0.0,\n",
    "                   cp.sum(P)-P_tot==0.0,\n",
    "                   cp.sum(W)-W_tot==0.0]\n",
    "    \n",
    "    prob = cp.Problem(objective, constraints)\n",
    "    prob.solve()\n",
    "      \n",
    "    return prob.status, -prob.value, P.value, W.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Example\n",
    "\n",
    "Consider the case where there are 5 channels, $n=5$, $\\alpha = \\beta = (2.0,2.2,2.4,2.6,2.8)$, $P_{\\text{tot}} = 0.5$ and $W_{\\text{tot}}=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: optimal\n",
      "Optimal utility value = 2.451\n",
      "Optimal power level:\n",
      "[1.151e-09 1.708e-09 2.756e-09 5.788e-09 5.000e-01]\n",
      "Optimal bandwidth:\n",
      "[3.091e-09 3.955e-09 5.908e-09 1.193e-08 1.000e+00]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(precision=3)\n",
    "n = 5               # number of receivers in the system\n",
    "\n",
    "a_val = np.arange(10,n+10)/(1.0*n)  # α\n",
    "b_val = np.arange(10,n+10)/(1.0*n)  # β\n",
    "P_tot = 0.5\n",
    "W_tot = 1.0\n",
    "status, utility, power, bandwidth = optimal_power(n, a_val, b_val, P_tot, W_tot)\n",
    "\n",
    "print('Status: {}'.format(status))\n",
    "print('Optimal utility value = {:.4g}'.format(utility))\n",
    "print('Optimal power level:\\n{}'.format(power))\n",
    "print('Optimal bandwidth:\\n{}'.format(bandwidth))"
   ]
  }
 ],
 "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
}