{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimize beamwidth of an array with arbitrary 2-D geometry\n",
    "\n",
    "A derivative work by Judson Wilson, 5/14/2014.<br>\n",
    "Adapted (with significant changes) from the CVX example of the same name, by Almir Mutapcic, 2/2/2006.\n",
    "\n",
    "\n",
    "Topic References:\n",
    "\n",
    "* \"Convex optimization examples\" lecture notes (EE364) by S. Boyd\n",
    "* \"Antenna array pattern synthesis via convex optimization\" by H. Lebret and S. Boyd\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This algorithm designs an antenna array such that:\n",
    "\n",
    "* it has unit sensitivity at some target direction  \n",
    "* it obeys a constraint on a minimum sidelobe level outside the beam\n",
    "* it minimizes the beamwidth of the pattern.\n",
    "\n",
    "This is a quasiconvex problem. Define the target direction as $\\theta_{\\mbox{tar}}$, and a beamwidth of $\\Delta \\theta_{\\mbox{bw}}$. The beam occupies the angular interval\n",
    "\n",
    "$$\\Theta_b = \\left(\\theta_{\\mbox{tar}}\n",
    "-\\frac{1}{2}\\Delta \\theta_{\\mbox{bw}},\\; \\theta_{\\mbox{tar}} \n",
    "+ \\frac{1}{2}\\Delta \\theta_{\\mbox{bw}}\\right).\n",
    "$$\n",
    "\n",
    "Solving for the minimum beamwidth $\\Delta \\theta_{\\mbox{bw}}$ is performed by bisection, where the interval which contains the optimal value is bisected  according to the result of the following feasibility problem:\n",
    "\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize}   &  0 \\\\\n",
    "\\mbox{subject to} & y(\\theta_{\\mbox{tar}}) = 1 \\\\\n",
    "                  &  \\left|y(\\theta)\\right| \\leq t_{\\mbox{sb}} \n",
    "                         \\quad \\forall \\theta \\notin \\Theta_b.\n",
    "\\end{array}\n",
    "    \n",
    "$y$ is the antenna array gain pattern (a complex-valued function), $t_{\\mbox{sb}}$ is the maximum allowed sideband gain threshold, and the variables are $w$ (antenna array weights or shading coefficients). The gain pattern is a linear function of $w$: $y(\\theta) = w^T a(\\theta)$ for some $a(\\theta)$ describing the antenna array configuration and specs.\n",
    "\n",
    "Once the optimal beamwidth is found, the solution $w$ is refined with the following optimization:\n",
    "\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize}   &  \\|w\\| \\\\\n",
    "\\mbox{subject to} & y(\\theta_{\\mbox{tar}}) = 1 \\\\\n",
    "                  & \\left|y(\\theta)\\right|  \\leq t_{\\mbox{sb}}\n",
    "                          \\quad \\forall \\theta \\notin \\Theta_b.\n",
    "\\end{array}\n",
    "\n",
    "The implementation below discretizes the angular quantities and their counterparts, such as $\\theta$.\n",
    "\n",
    "## Problem specification and data\n",
    "\n",
    "### Antenna array selection\n",
    "\n",
    "Choose either:\n",
    "\n",
    "* A random 2D positioning of antennas.\n",
    "* A uniform 1D positioning of antennas along a line.\n",
    "* A uniform 2D positioning of antennas along a grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "# Select array geometry:\n",
    "ARRAY_GEOMETRY = '2D_RANDOM'\n",
    "#ARRAY_GEOMETRY = '1D_UNIFORM_LINE'\n",
    "#ARRAY_GEOMETRY = '2D_UNIFORM_LATTICE'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Problem specs.\n",
    "#\n",
    "lambda_wl = 1         # wavelength\n",
    "theta_tar = 0           # target direction\n",
    "min_sidelobe = -20    # maximum sidelobe level in dB\n",
    "\n",
    "max_half_beam = 50    # starting half beamwidth (must be feasible)\n",
    "\n",
    "#\n",
    "# 2D_RANDOM: \n",
    "#     n randomly located elements in 2D.\n",
    "#\n",
    "if ARRAY_GEOMETRY == '2D_RANDOM':\n",
    "    # Set random seed for repeatable experiments.\n",
    "    np.random.seed(1)\n",
    "    # Uniformly distributed on [0,L]-by-[0,L] square.\n",
    "    n = 36\n",
    "    L = 5\n",
    "    loc = L*np.random.random((n,2))\n",
    "\n",
    "#\n",
    "# 1D_UNIFORM_LINE:\n",
    "#     Uniform 1D array with n elements with inter-element spacing d.\n",
    "#\n",
    "elif ARRAY_GEOMETRY == '1D_UNIFORM_LINE':\n",
    "    n = 30\n",
    "    d = 0.45*lambda_wl\n",
    "    loc = np.hstack(( d * np.array(range(0,n)).reshape(-1, 1), \\\n",
    "                          np.zeros((n,1)) ))\n",
    "\n",
    "#\n",
    "# 2D_UNIFORM_LATTICE:\n",
    "#     Uniform 2D array with m-by-m element with d spacing.\n",
    "#\n",
    "elif ARRAY_GEOMETRY == '2D_UNIFORM_LATTICE':\n",
    "    m = 6\n",
    "    n = m**2\n",
    "    d = 0.45*lambda_wl\n",
    "\n",
    "    loc = np.zeros((n, 2))\n",
    "    for x in range(m):\n",
    "        for y in range(m):\n",
    "            loc[m*y+x,:] = [x,y]\n",
    "    loc = loc*d\n",
    "\n",
    "else:\n",
    "    raise Exception('Undefined array geometry')\n",
    "\n",
    "\n",
    "#\n",
    "# Construct optimization data.\n",
    "#\n",
    "\n",
    "# Build matrix A that relates w and y(theta), ie, y = A*w.\n",
    "theta = np.array(range(1, 360+1)).reshape(-1, 1)\n",
    "A = np.kron(np.cos(np.pi*theta/180), loc[:, 0].T) \\\n",
    "  + np.kron(np.sin(np.pi*theta/180), loc[:, 1].T)\n",
    "A = np.exp(2*np.pi*1j/lambda_wl*A)\n",
    "\n",
    "# Target constraint matrix.\n",
    "ind_closest = np.argmin(np.abs(theta - theta_tar))\n",
    "Atar = A[ind_closest,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve using bisection algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We are only considering integer values of the half beam-width\n",
      "(since we are sampling the angle with 1 degree resolution).\n",
      "\n",
      "Problem is feasible for half beam-width = 26.0 degress\n",
      "Problem is feasible for half beam-width = 14.0 degress\n",
      "Problem is not feasible for half beam-width = 8.0 degress\n",
      "Problem is feasible for half beam-width = 11.0 degress\n",
      "Problem is feasible for half beam-width = 10.0 degress\n",
      "Problem is not feasible for half beam-width = 9.0 degress\n",
      "Optimum half beam-width for given specs is 10.0\n",
      "final objective value: 1.33693666264\n"
     ]
    }
   ],
   "source": [
    "# Bisection range limits. Reduce by half each step.\n",
    "halfbeam_bot = 1\n",
    "halfbeam_top = max_half_beam\n",
    "\n",
    "print('We are only considering integer values of the half beam-width')\n",
    "print('(since we are sampling the angle with 1 degree resolution).')\n",
    "print('')\n",
    "\n",
    "# Iterate bisection until 1 angular degree of uncertainty.\n",
    "while halfbeam_top - halfbeam_bot > 1:\n",
    "    # Width in degrees of the current half-beam.\n",
    "    halfbeam_cur = np.ceil( (halfbeam_top + halfbeam_bot)/2.0 )\n",
    "\n",
    "    # Create optimization matrices for the stopband,\n",
    "    # i.e. only A values for the stopband angles.\n",
    "    upper_beam = (theta_tar + halfbeam_cur) % 360\n",
    "    lower_beam = (theta_tar - halfbeam_cur) % 360\n",
    "    if upper_beam > lower_beam:\n",
    "        ind = np.nonzero(np.squeeze(np.array(np.logical_or( \\\n",
    "                        theta <= lower_beam, \\\n",
    "                        theta >= upper_beam ))))\n",
    "    else:\n",
    "        ind = np.nonzero(np.squeeze(np.array(np.logical_and( \\\n",
    "                        theta <= lower_beam, \\\n",
    "                        theta >= upper_beam ))))\n",
    "    \n",
    "    As = A[ind[0],:]\n",
    "    \n",
    "    #\n",
    "    # Formulate and solve the feasibility antenna array problem.\n",
    "    #\n",
    "\n",
    "    # As of this writing (2014/05/14) cvxpy does not do complex valued math,\n",
    "    # so the real and complex values must be stored separately as reals\n",
    "    # and operated on as follows:\n",
    "    #     Let any vector or matrix be represented as a+bj, or A+Bj.\n",
    "    #     Vectors are stored [a; b] and matrices as [A -B; B A]:\n",
    "    \n",
    "    # Atar as [A -B; B A]\n",
    "    Atar_R = Atar.real\n",
    "    Atar_I = Atar.imag\n",
    "    neg_Atar_I = -Atar_I\n",
    "    Atar_RI = np.block([[Atar_R, neg_Atar_I], [Atar_I, Atar_R]])\n",
    "\n",
    "    # As as [A -B; B A]\n",
    "    As_R = As.real\n",
    "    As_I = As.imag\n",
    "    neg_As_I = -As_I\n",
    "    As_RI = np.block([[As_R, neg_As_I], [As_I, As_R]])\n",
    "    As_RI_top = np.block([As_R, neg_As_I])\n",
    "    As_RI_bot = np.block([As_I, As_R])\n",
    "\n",
    "    # 1-vector as [1, 0] since no imaginary part\n",
    "    realones_ri = np.array([1.0, 0.0])\n",
    "\n",
    "    # Create cvxpy variables and constraints\n",
    "    w_ri = cp.Variable(shape=(2*n))\n",
    "    constraints = [ Atar_RI*w_ri == realones_ri]\n",
    "    # Must add complex valued constraint \n",
    "    # abs(As*w <= 10**(min_sidelobe/20)) row by row by hand.\n",
    "    # TODO: Future version use norms() or complex math\n",
    "    # when these features become available in cvxpy.\n",
    "    for i in range(As.shape[0]):\n",
    "        #Make a matrix whos product with w_ri is a 2-vector\n",
    "        #which is the real and imag component of a row of As*w\n",
    "        As_ri_row = np.vstack((As_RI_top[i, :], As_RI_bot[i, :]))\n",
    "        constraints.append( \\\n",
    "                cp.norm(As_ri_row*w_ri) <= 10**(min_sidelobe/20) )\n",
    "\n",
    "    # Form and solve problem.\n",
    "    obj = cp.Minimize(0)\n",
    "    prob = cp.Problem(obj, constraints)\n",
    "    prob.solve(solver=cp.CVXOPT)\n",
    "\n",
    "    # Bisection (or fail).\n",
    "    if prob.status == cp.OPTIMAL:\n",
    "        print('Problem is feasible for half beam-width = {}'\n",
    "              ' degress'.format(halfbeam_cur))\n",
    "        halfbeam_top = halfbeam_cur\n",
    "    elif prob.status == cp.INFEASIBLE:\n",
    "        print('Problem is not feasible for half beam-width = {}'\n",
    "              ' degress'.format(halfbeam_cur))\n",
    "        halfbeam_bot = halfbeam_cur\n",
    "    else:\n",
    "        raise Exception('CVXPY Error')\n",
    "\n",
    "# Optimal beamwidth.\n",
    "halfbeam = halfbeam_top\n",
    "upper_beam = (theta_tar + halfbeam) % 360\n",
    "lower_beam = (theta_tar - halfbeam) % 360\n",
    "print('Optimum half beam-width for given specs is {}'.format(halfbeam))\n",
    "\n",
    "# Compute the minimum noise design for the optimal beamwidth\n",
    "if upper_beam > lower_beam:\n",
    "    ind = np.nonzero(np.squeeze(np.array(np.logical_or( \\\n",
    "                    theta <= lower_beam, \\\n",
    "                    theta >= upper_beam ))))\n",
    "else:\n",
    "    ind = np.nonzero(np.squeeze(np.array(np.logical_and( \\\n",
    "                    theta <= lower_beam, \\\n",
    "                    theta >= upper_beam ))))\n",
    "        \n",
    "As = A[ind[0],:]\n",
    "\n",
    "# As as [A -B; B A]\n",
    "# See earlier calculations for real/imaginary representation\n",
    "As_R = As.real\n",
    "As_I = As.imag\n",
    "neg_As_I = -As_I\n",
    "As_RI = np.block([[As_R, neg_As_I], [As_I, As_R]])\n",
    "As_RI_top = np.block([As_R, neg_As_I])\n",
    "As_RI_bot = np.block([As_I, As_R])\n",
    "\n",
    "constraints = [ Atar_RI*w_ri == realones_ri]\n",
    "# Same constraint as a above, on new As (hense different\n",
    "# actual number of constraints). See comments above.\n",
    "for i in range(As.shape[0]):\n",
    "    As_ri_row = np.vstack((As_RI_top[i, :], As_RI_bot[i, :]))\n",
    "    constraints.append( \\\n",
    "        cp.norm(As_ri_row*w_ri) <= 10**(min_sidelobe/20) )\n",
    "\n",
    "# Form and solve problem.\n",
    "# Note the new objective!\n",
    "obj = cp.Minimize(cp.norm(w_ri))\n",
    "prob = cp.Problem(obj, constraints)\n",
    "prob.solve(solver=cp.SCS)\n",
    "\n",
    "#if prob.status != cp.OPTIMAL:\n",
    "#    raise Exception('CVXPY Error')\n",
    "print(\"final objective value: {}\".format(obj.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13352 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 45400 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 47544 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plot inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "\n",
    "#\n",
    "# First Figure: Antenna Locations\n",
    "#\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.scatter(loc[:, 0], loc[:, 1],\n",
    "            s=30, facecolors='none', edgecolors='b')\n",
    "plt.title('Antenna Locations', fontsize=16)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "#\n",
    "# Second Plot: Array Pattern\n",
    "#\n",
    "\n",
    "# Complex valued math to calculate y = A*w_im;\n",
    "# See comments in code above regarding complex representation as reals.\n",
    "A_R = A.real\n",
    "A_I = A.imag\n",
    "neg_A_I = -A_I\n",
    "A_RI = np.block([[A_R, neg_A_I], [A_I, A_R]]);\n",
    "\n",
    "y = A_RI.dot(w_ri.value)\n",
    "y = y[0:int(y.shape[0]/2)] + 1j*y[int(y.shape[0]/2):] #now native complex\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "ymin, ymax = -40, 0\n",
    "plt.plot(np.arange(360)+1, np.array(20*np.log10(np.abs(y))))\n",
    "plt.plot([theta_tar, theta_tar], [ymin, ymax], 'g--')\n",
    "plt.plot([upper_beam, upper_beam], [ymin, ymax], 'r--')\n",
    "plt.plot([lower_beam, lower_beam], [ymin, ymax], 'r--')\n",
    "plt.xlabel('look angle', fontsize=16)\n",
    "plt.ylabel(r'mag $y(\\theta)$ in dB', fontsize=16)\n",
    "plt.ylim(ymin, ymax)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "#\n",
    "# Third Plot: Polar Pattern\n",
    "#\n",
    "plt.figure(figsize=(6,6))\n",
    "zerodB = 50\n",
    "dBY = 20*np.log10(np.abs(y)) + zerodB\n",
    "plt.plot(dBY * np.cos(np.pi*theta.flatten()/180),\n",
    "         dBY * np.sin(np.pi*theta.flatten()/180))\n",
    "plt.xlim(-zerodB, zerodB)\n",
    "plt.ylim(-zerodB, zerodB)\n",
    "plt.axis('off')\n",
    "\n",
    "# 0 dB level.\n",
    "plt.plot(zerodB*np.cos(np.pi*theta.flatten()/180),\n",
    "         zerodB*np.sin(np.pi*theta.flatten()/180), 'k:')\n",
    "plt.text(-zerodB,0,'0 dB', fontsize=16)\n",
    "# Max sideband level.\n",
    "m=min_sidelobe + zerodB\n",
    "plt.plot(m*np.cos(np.pi*theta.flatten()/180),\n",
    "         m*np.sin(np.pi*theta.flatten()/180), 'k:')\n",
    "plt.text(-m,0,'{:.1f} dB'.format(min_sidelobe), fontsize=16)\n",
    "#Lobe center and boundaries angles.\n",
    "theta_1 = theta_tar+halfbeam\n",
    "theta_2 = theta_tar-halfbeam\n",
    "plt.plot([0, 55*np.cos(theta_tar*np.pi/180)], \\\n",
    "         [0, 55*np.sin(theta_tar*np.pi/180)], 'k:')\n",
    "plt.plot([0, 55*np.cos(theta_1*np.pi/180)], \\\n",
    "         [0, 55*np.sin(theta_1*np.pi/180)], 'k:')\n",
    "plt.plot([0, 55*np.cos(theta_2*np.pi/180)], \\\n",
    "         [0, 55*np.sin(theta_2*np.pi/180)], 'k:')\n",
    "\n",
    "#Show plot.\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
