Minimize beamwidth of an array with arbitrary 2-D geometry¶
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:
- "Convex optimization examples" lecture notes (EE364) by S. Boyd
- "Antenna array pattern synthesis via convex optimization" by H. Lebret and S. Boyd
Introduction¶
This algorithm designs an antenna array such that:
- it has unit sensitivity at some target direction
- it obeys a constraint on a minimum sidelobe level outside the beam
- it minimizes the beamwidth of the pattern.
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
$$\Theta_b = \left(\theta_{\mbox{tar}} -\frac{1}{2}\Delta \theta_{\mbox{bw}},\; \theta_{\mbox{tar}} + \frac{1}{2}\Delta \theta_{\mbox{bw}}\right). $$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:
\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}$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.
Once the optimal beamwidth is found, the solution $w$ is refined with the following optimization:
\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 $\theta$.
Problem specification and data¶
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.
In [1]:
import cvxpy as cvx
import numpy as np
# Select array geometry:
ARRAY_GEOMETRY = '2D_RANDOM'
#ARRAY_GEOMETRY = '1D_UNIFORM_LINE'
#ARRAY_GEOMETRY = '2D_UNIFORM_LATTICE'
Data generation¶
In [2]:
#
# 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.matrix(range(0,n)).T, \
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.matrix(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.mat(range(1, 360+1)).T
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¶
In [3]:
# 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 seperately 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.bmat('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.bmat('As_R neg_As_I; As_I As_R')
As_RI_top = np.bmat('As_R neg_As_I')
As_RI_bot = np.bmat('As_I As_R')
# 1-vector as [1; 0] since no imaginary part
realones_ri = np.mat( np.vstack( \
(np.ones(Atar.shape[0]),
np.zeros(Atar.shape[0])) ))
# Create cvxpy variables and constraints
w_ri = cvx.Variable(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( \
cvx.norm(As_ri_row*w_ri) <= 10**(min_sidelobe/20) )
# Form and solve problem.
obj = cvx.Minimize(0)
prob = cvx.Problem(obj, constraints)
prob.solve(solver=cvx.CVXOPT)
# Bisection (or fail).
if prob.status == cvx.OPTIMAL:
print ('Problem is feasible for half beam-width = {}'
' degress').format(halfbeam_cur)
halfbeam_top = halfbeam_cur
elif prob.status == cvx.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.bmat('As_R neg_As_I; As_I As_R')
As_RI_top = np.bmat('As_R neg_As_I')
As_RI_bot = np.bmat('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( \
cvx.norm(As_ri_row*w_ri) <= 10**(min_sidelobe/20) )
# Form and solve problem.
# Note the new objective!
obj = cvx.Minimize(cvx.norm(w_ri))
prob = cvx.Problem(obj, constraints)
prob.solve(solver=cvx.SCS)
#if prob.status != cvx.OPTIMAL:
# raise Exception('CVXPY Error')
Out[3]:
Result plots¶
In [4]:
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(np.array(loc[:, 0]), np.array(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.bmat('A_R neg_A_I; A_I A_R');
y = A_RI*w_ri.value
y = y[0:y.shape[0]/2] + 1j*y[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(np.array(dBY)*np.array(np.cos(np.pi*theta/180)), \
np.array(dBY)*np.array(np.sin(np.pi*theta/180)))
plt.xlim(-zerodB, zerodB)
plt.ylim(-zerodB, zerodB)
plt.axis('off')
# 0 dB level.
plt.plot(zerodB*np.array(np.cos(np.pi*theta/180)), \
zerodB*np.array(np.sin(np.pi*theta/180)), 'k:')
plt.text(-zerodB,0,'0 dB', fontsize=16)
# Max sideband level.
m=min_sidelobe + zerodB
plt.plot(m*np.array(np.cos(np.pi*theta/180)), \
m*np.array(np.sin(np.pi*theta/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()
