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