array pattern synthesis part ii: optimization -凯发k8网页登录

this example uses:

this example shows how to use optimization techniques to perform pattern synthesis. using minimum variance beamforming as an illustration, the example covers how to set up the pattern synthesis as an optimization problem and then solve it using optimization solvers.

pattern synthesis can be achieved using a variety of techniques. the example introduced nulling and windowing techniques.

this example requires optimization toolbox™.

minimum variance beamforming

compared to a single element, a major benefit of a phased array is the capability of forming beams. in general, an n-element phased array provides n degrees of freedom when synthesizing a pattern, meaning that one can adjust n weights, one for each element, to control the shape of the beam to satisfy some pre-defined constraints.

many popular beamforming techniques can be expressed as optimization problems. for example, the popular minimum variance distortionless response (mvdr) beamformer is used to minimize the total noise output while preserving the signal from a given direction. mathematically, the mvdr beamforming weights can be found by solving the optimization problem

$\begin{array}{l} w_{o} = \min_{w} w^{h} r w \\ s . t . b (θ_{0}, w) = 1 \end{array}$

where $r$ is the covariance matrix, $b (θ_{0}, w)$ is the beam response towards the direction $θ_{0}$ with the weights vector $w$ , and $w_{o}$ is the vector containing the optimal weights. because the beam pattern can be written as

$b (θ, w) = w^{h} v (θ)$

where $v (θ)$ is the steering vector corresponding to the angle $θ$ , the constraints have a linear relation in the complex domain.

the distortionless constraint at $θ_{0}$ only takes 1 degree of freedom in the design space. since an n-element array can handle n – 1 constraints, the mvdr beamformer can be extended to include more constraints, which results in a linear constraints minimum variance (lcmv) beamformer. the extra constraints in an lcmv beamformer are often used to null out interferences from given directions. therefore, such an lcmv beamformer is a solution to the following optimization problem

$\begin{array}{l} w_{o} = \min_{w} w^{h} r w \\ s . t . b (θ_{0}, w) = 1 \\ b (θ_{k}, w) = 0 k = 1, 2, \dots, k \end{array}$

where $θ_{k}$ represents the $k$ th nulling direction and there are $k$ nulling directions in total. because all constraints are linear in the complex domain, the lcmv beamforing weights can be derived analytically using the lagrange multiplier method.

next, a 32-element ula with half wavelength spacing is used to demonstrate how to perform array synthesis using the lcmv approach. note that even though this example uses a linear array, the technique applies to any array geometry. assume that the signal of interest is at 0 degrees azimuth and the interferences are at –70, -–40, and –20 degrees azimuth. the noise power is assumed to be at 40 db below signal. the closed form solution for lcmv weights can be computed using the lcmvweights function.

n = 32;
pos = (0:n-1)*0.5;         % position of elements
ang_i = [-70 -40 -20];     % interference angles
ang_d = 0;                 % desired angle
rn = sensorcov(pos,ang_i,db2pow(-40));  % noise covariance matrix
sv_c = steervec(pos,[ang_d ang_i]);     % linear constraints
r_c = [1 zeros(size(ang_i))]';          % desired response
w_lcmv = lcmvweights(sv_c,r_c,rn);      % lcmv weights
ang_plot = -90:0.1:90;
sv_plot = steervec(pos,ang_plot);
plcmv = plot(ang_plot,mag2db(abs(w_lcmv'*sv_plot)));
xline(ang_i,"--")
ylim([-100 0])
grid on
legend(plcmv,"lcmv - analytic")
xlabel("azimuth angle (deg)")
ylabel("beam pattern (db)")

figure contains an axes object. the axes object with xlabel azimuth angle (deg), ylabel beam pattern (db) contains 4 objects of type line, constantline. this object represents lcmv - analytic.

the figure shows that the computed weights satisfy all constraints.

lcmv beamforming using quadratic programming

since lcmv beamforming can be described as an optimization problem, it would be useful to see how the problem can be solved using optimization techniques. when a closed form solution is not available, for example when constraints include inequalities, such workflow becomes important.

when using optimization techniques to perform pattern synthesis tasks, the first step is to pick a solver. since solvers are often developed to address a particular problem formulation, it is of utmost importance to make the right choice.

the quadratic programming (qp) solver is often used to solve a quadratic objective function with linear constraints. the formulation is given by

$\begin{array}{l} w_{o} = \min_{w} \frac{1}{2} w^{t} h w f^{t} w \\ s . t . a \cdot w \leq b \\ a_{eq} \cdot w = b_{eq} \\ w_{lb} \leq w \leq w_{ub} \end{array}$

where $a$ and $a_{eq}$ are matrices specifying the inequality and equality constraints, respectively; $b$ and $b_{eq}$ are vectors specifying the desired bounds for inequality and equality constraints, respectively; and $w_{lb}$ and $w$ indicate lower and upper bounds for elements in $w$ , respectively.

at first sight, the problem formulation of quadratic programming seems to be a perfect match for lcm beamforming because we know that the beam pattern can be written as $b (θ, w) = w^{h} v (θ)$ . unfortunately, this is not true because although the beam pattern is linear in complex domain, it is not so in real domain. on the other hand, most optimization solvers, including quadratic programming, work in real domain as the solution space becomes really complicated in the complex domain. in addition, the inequality is not well defined in the complex domain either. thus, the question is if the lcmv beamforming problem can be formulated into real domain and still obtain the quadratic programming form.

start with the objective function $w^{h} r w$ . note that $r$ is a covariance matrix in this context, thus $w^{h} r w$ is real. therefore, the goal is really to minimize the real part of $w^{h} r w$ , which can be written in the following form:

$ℜ {w^{h} r w} = [ℜ {w}^{t} ℑ {w}^{t}] [\begin{array}{c} ℜ {r} & - ℑ {r} \\ ℑ {r} & ℜ {r} \end{array}] [\begin{array}{c} ℜ {w} \\ ℑ {w} \end{array}] = {\overline{w}}^{t} \overline{r} \overline{w}$

where $\overline{w} = [\begin{array}{c} ℜ {w} \\ ℑ {w} \end{array}]$ and $\overline{r} = [\begin{array}{c} ℜ {r} & - ℑ {r} \\ ℑ {r} & ℜ {r} \end{array}]$ .

next, look at the constraints. the fact that the constraints in lcmv beamforming are equalities makes the conversion a bit easier since one complex linear equality constraint can just be divided into two corresponding real linear equality constraints, one for the real part and the other for the imaginary part. for example, the following constraint

$w^{h} v = r^{h}$

can be rewritten as

$\begin{array}{l} v^{h} w = r \\ [\begin{array}{c} ℜ {v} & ℑ {v} \\ ℑ {v} & - ℜ {v} \end{array}] [\begin{array}{c} ℜ {w} \\ ℑ {w} \end{array}] = [\begin{array}{c} ℜ {v} & ℑ {v} \\ ℑ {v} & - ℜ {v} \end{array}] \overline{w} = [\begin{array}{c} ℜ {r} \\ ℑ {r} \end{array}] \end{array}$

with these relations, quadratic programming can then be used to obtain lcmv beamforming weights.

% objective function
h = [real(rn) -imag(rn);imag(rn) real(rn)];
f = [];
% equality constraints
aeq = [real(sv_c).' imag(sv_c).';imag(sv_c).' -real(sv_c).'];
beq = [real(r_c);imag(r_c)];
% compute weights
x_opt = quadprog(h,f,[],[],aeq,beq);

minimum found that satisfies the constraints.
optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

w_opt = x_opt(1:numel(x_opt)/2) x_opt(numel(x_opt)/2 1:end)*1i;
hold on
plcmvq = plot(ang_plot,mag2db(abs(w_opt'*sv_plot)),"--");
hold off
legend([plcmv plcmvq],["lcmv - lagrange" "lcmv - qp"])

figure contains an axes object. the axes object with xlabel azimuth angle (deg), ylabel beam pattern (db) contains 5 objects of type line, constantline. these objects represent lcmv - lagrange, lcmv - qp.

the figure shows that the qp technique can obtain the same lcmv beamforming weights.

minimum variance beamforming using quadratic programming

another common requirement for pattern synthesis is to ensure the response is below a certain threshold for a set of given angles. these are inequality requirements and the overall optimization problem becomes

$\begin{array}{l} w_{o} = \min_{w} w^{h} r w \\ s . t . b (θ_{0}, w) = 1 \\ b (θ_{k}, w) = 0 k = 1, 2, \dots, k \\ | b (w, θ_{l}) | \leq r_{l} l = 1, 2, \dots, l \end{array}$

where $θ_{l}$ is the $l$ th angle at which the sidelobe level should be less than the threshold $r_{l}$ . there are $l$ inequality requirements in total.

as an example, add the following requirements to our sample problem:

the sidelobe should be below –40 db between –30 and –10 degrees in azimuth.
everywhere outside of the mainlobe, the sidelobe should be below −20 db.

this code plots the mask of the desired pattern and returns the inequality constraints:

[rc,angc,sv_cineq,r_ineq,ang_ineq]=generatepatternmask(pos,r_c,[ang_d ang_i]);
pdp = plot(angc,rc,"k--");
ylim([-100 0])
grid on
xlabel("azimuth angle (deg)")
ylabel("beam pattern (db)")

figure contains an axes object. the axes object with xlabel azimuth angle (deg), ylabel beam pattern (db) contains an object of type line.

with the introduction of the new inequality constraint, the problem no longer has a closed-form solution. in addition, due to the presence of the norm operation, the aforementioned trick to translate complex constraints to real domain no longer works. therefore, a new way to translate the complex inequality is needed to use the qp solver for this set of constraints,

given a constraint of $| b (θ, w) | < r$ , one possible way to convert this is to assume that both real and imaginary parts contribute equally to the norm. in that case, the single complex inequality constraint can be divided into two real inequalities as

$\begin{array}{l} - \frac{r}{\sqrt{2}} \leq ℜ {b (θ, w)} \leq \frac{r}{\sqrt{2}} \\ - \frac{r}{\sqrt{2}} \leq ℑ {b (θ, w)} \leq \frac{r}{\sqrt{2}} \end{array}$

$\begin{array}{l} ℜ {b (θ, w)} \leq \frac{r}{\sqrt{2}} \\ ℑ {b (θ, w)} \leq \frac{r}{\sqrt{2}} \\ - ℜ {b (θ, w)} \leq \frac{r}{\sqrt{2}} \\ - ℑ {b (θ, w)} \leq \frac{r}{\sqrt{2}} \end{array}$

the real and imaginary part of the constraints can be extracted as described in the previous section. hence, the following code sets up the constraints and solves the pattern synthesis task.

% map complex inequality constraints to real constraints
sv_c_real = [real(sv_cineq).' imag(sv_cineq).'];
sv_c_imag = [imag(sv_cineq).' -real(sv_cineq).'];
a_ineq = [sv_c_real;sv_c_imag];
a_ineq = [a_ineq;-a_ineq];
b_ineq = [r_ineq;r_ineq]/sqrt(2);  
b_ineq = [b_ineq;b_ineq];
% compute weights
x_opt = quadprog(h,f,a_ineq,b_ineq,aeq,beq);

minimum found that satisfies the constraints.
optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

% reassemble complex weights from real and imaginary parts
w_opt = x_opt(1:numel(x_opt)/2) x_opt(numel(x_opt)/2 1:end)*1i;
hold on
pqp = plot(ang_plot,mag2db(abs(w_opt'*sv_plot)));
hold off
legend(pqp,"pattern - qp")

figure contains an axes object. the axes object with xlabel azimuth angle (deg), ylabel beam pattern (db) contains 2 objects of type line. this object represents pattern - qp.

the resulting pattern satisfies the requirements.

minimum variance beamforming using second order cone programming

even though the qp solver solves the minimum variance beamforming problem successfully, the conversion of the inequality constraints from complex to real is not ideal. for example, there is no reason that real and imaginary parts should always equally contribute to the norm. however, there is really no easy way to handle a constraint with norm in a qp problem setup.

luckily, there is another optimization solver that can be used in this case, called second order cone programming (socp). socp is developed to efficiently solve the following optimization problem:

$\begin{array}{l} w_{o} = \min_{w} f^{t} w \\ s . t . ‖ a_{sc} \cdot w - b_{sc} ‖ \leq d^{t} w - γ \\ a \cdot w \leq b \\ a_{eq} \cdot w = b_{eq} \\ w_{lb} \leq w \leq w_{ub} \end{array}$

where $a_{sc}$ is a matrix, $b_{sc}$ and $d$ are vectors, and $γ$ is a scalar. $a_{sc}$ , $b_{sc}$ , $d$ , and $γ$ together define a second order cone constraint, which is a natural choice to specify a constraint on norms. the solver can handle multiple second order cone constraints simultaneously.

comparing to the problem formulation for qp, socp can not only handle second order constraints, but also handle all constraints specified in qp. however, the objective function is quite different. instead of minimizing a variance, socp minimize a linear function of the optimization variable, $w$ . fortunately, there is a way to convert the minimum variance objective function to fit into this socp framework. the trick is to define an auxiliary variable, $t$ , and ensure that the variance is bounded by $t$ . this way, the problem becomes to minimize $t$ while the original quadratic objective becomes a cone constraint bounded by $t$ .

following this thought, amend $t$ to the bottom of $\overline{w}$ . now the optimization search space is defined by $\tilde{w} = [\begin{array}{c} \overline{w} \\ t \end{array}]$ . since the goal is to minimize $t$ , $f$ is defined as $f^{t} = [0, \dots, 0 1]$ .

f = [zeros(2*n,1);1];

next, explore how to convert the quadratic constraint to a second order cone constraint. since the idea is to constraint the variance as

${\overline{w}}^{t} \overline{r} \overline{w} \leq t$

it is equivalent to say

${‖ {\overline{r}}^{0.5} \overline{w} ‖}^{2} \leq t \frac{1}{4}$

thus,

${‖ [\begin{array}{c} {\overline{r}}^{0.5} & 0 \\ 0 & 1 \end{array}] [\begin{array}{c} \overline{w} \\ t \end{array}] ‖}^{2} \approx {‖ {\overline{r}}^{0.5} \overline{w} ‖}^{2} t^{2} \leq t^{2} t \frac{1}{4} = {(t \frac{1}{2})}^{2}$

and the second cone constraint becomes

$‖ [\begin{array}{c} {\overline{r}}^{0.5} & 0 \\ 0 & 1 \end{array}] \tilde{w} ‖ \leq t \frac{1}{2}$

% the cone constraints are defined as ||a_c*u-b_c||<=d_c^t*u-gamma_c
ac = [sqrtm(h) zeros(2*n,1);zeros(1,2*n) 1];
bc = zeros(2*n 1,1);
dc = [zeros(2*n,1);1];
gammac = -1/2;
m = numel(r_ineq);
socconstraints(m 1) = secondordercone(ac,bc,dc,gammac);

the remaining second cone constraints are from the original inequality constraints and are fairly straightforward to be converted. each inequality constraint in the form of

$‖ b (θ, w) ‖ \leq r$

can be written as a second order cone constraint as

$‖ [\begin{array}{c} ℜ {b (θ, w)} & 0 \\ ℑ {b (θ, w)} & 0 \\ 0 & 0 \end{array}] \tilde{w} ‖ \leq r$

for m = 1:m
    ac = [a_ineq([m m m],:) zeros(2,1);zeros(1,2*n) 0];
    bc = zeros(3,1);
    dc = zeros(2*n 1,1);
    gammac = -r_ineq(m);
    socconstraints(m) = secondordercone(ac,bc,dc,gammac);
end

the equality constraints also need to be augmented to accommodate the auxiliary variable, $t$ .

% patch aeq and beq
aeqc = [aeq zeros(size(aeq,1),1);zeros(1,2*n) 0]; 
beqc = [beq;0];

the array weights can now be computed using the socp solver. since the inequality constraints are more natural to be represented in cone constraints, naively one would expect that in this example, the resulting pattern from socp would be better compared to the pattern from qp.

[u_op,feval,exitflag] = coneprog(f,socconstraints,[],[],aeqc,beqc);

optimal solution found.

% reassemble complex weights from real and imaginary parts 
% the last entry in the optimized vector corresponds to the variable t
w_op = u_op(1:n) u_op(n 1:2*n)*1i;
hold on
psocp = plot(ang_plot,mag2db(abs(w_op'*sv_plot)));
hold off
legend([pqp psocp],["pattern - qp" "pattern - socp"])

figure contains an axes object. the axes object with xlabel azimuth angle (deg), ylabel beam pattern (db) contains 3 objects of type line. these objects represent pattern - qp, pattern - socp.

indeed, the plot shows that the pattern obtained via socp has better sidelobe suppression compared to the pattern obtained using qp. as you can see, the process is quite involved when an optimization solver is used. we can simplify this process by calling the minvarweights function. the code below shows how to define the mask constraints as inputs for minvarweights function. the mask sidelobe level and mask angle represent the inequality constraints for the optimization problem, and the null angle represents the equality constraint for socp. for more details, refer to the documentation.

w_opt = minvarweights(pos,ang_d,rn, ...
    masksidelobelevel=mag2db(r_ineq)', ...
    maskangle=ang_ineq,nullangle=ang_i);
plot(ang_plot,mag2db(abs(w_opt'*sv_plot)))
hold on
plot(angc,rc,'k--')
hold off
legend(["pattern - minvarweights" "sidelobe mask"], ...
    location="best")
ylim([-100 0])
grid on
xlabel("azimuth angle (deg)")
ylabel("beam pattern (db)")

figure contains an axes object. the axes object with xlabel azimuth angle (deg), ylabel beam pattern (db) contains 2 objects of type line. these objects represent pattern - minvarweights, sidelobe mask.

minimum variance beamforming using `minvarweights` with other mask constraints

in this section, we use minvarweights function to solve pattern synthesis problems with different optimization constraints. consider a linear array with 55 antenna elements and 0.5λ spacing ula in the direction of 25-degree in azimuth. the goal is to optimize the beamforming weights with the following constraints:

sidelobe levels less than a tapered sidelobe mask decrease linearly from −18db to −55db in the range of 1st sll to ±90 degrees in azimuth.
steer nulls at [−45 −35 40 60] degrees in azimuth.

n = 55;                     % number of elements
pos = (0:n-1)*0.5;          % position of elements
angmainbeam = 25;           % main beam direction
% define null directions 
angn = [-35 -45 40 60];     % null angles, degree
% define the sll mask  
angml = -90:.2:22.5;                        % mask angles for left slls
sllml = linspace(-55,-18,numel(angml));     % mask slls for angml ( decrease linearly from −18db to −55db)
angmr = 27.5:0.2:90;                        % mask angles for right slls
sllmr = linspace(-18,-55,numel(angmr));     % mask slls for angmr ( decrease linearly from −18db to −55db)
angm = [angml angmr];                       % mask angles
sllm = [sllml sllmr];                       % mask sll
% compute optimized weights 
optimizedweights = minvarweights(pos, angmainbeam, ...
    maskangle=angm,masksidelobelevel=sllm,nullangle=angn);
% calculate array pattern
sv_plot = steervec(pos,-90:.25:90);
pattern = optimizedweights'*sv_plot;
% plot the pattern
plot(-90:.25:90,mag2db(abs(pattern)))
xlabel("azimuth angle (deg)")
ylabel("beam pattern (db)")
axis([-90 90 -80 0])
grid on
% plot the mask with nulls
hold on
[angmask,ia] = sort([angm angn]);
sllmask = [sllm -1000*ones(1,size(angn,2))];
plot(angmask,sllmask(ia),'k')
legend(["pattern - minvarweights" "sidelobe mask"], ...
    location="best")
hold off

this animation shows the pattern synthesis results for the beam sweeping from −35 degree to 35 degree and slls decrease linearly from −18 db to −55 db in the range of 1st sll to ±90 degrees in azimuth.

summary

this example shows how to use optimization techniques to perform pattern synthesis. the example started with the familiar lcmv problem and extended the desired pattern to include inequality constraints. compared to the closed-form solution, optimization techniques are more flexible and can be applied to address more complex constraints. the example also showed how to address inequality constraints using either quadratic programming, second order cone programming, or minvarweights.

reference

[1] lebret, hervé, and stephen boyd. "antenna array pattern synthesis via convex optimization." ieee transactions on signal processing, vol. 45, no. 3 (march 1997), pp. 526–532.

function [rcplot,angc,svcineq,rineq,angineq]=generatepatternmask(pos,rcexist,angexist)
% generatepatternmask  generate mask for pattern synthesis
ang_c1 = [-30:-21 -19:-10];
ang_c2 = [-90:-71 -69:-41 -39:-31 -9:-5 5:90];  
% note the main beam is between -10 and 10 degrees
sv_c1 = steervec(pos,ang_c1);
sv_c2 = steervec(pos,ang_c2);
r1 = db2mag(-40)*ones(size(sv_c1,2),1);
r2 = db2mag(-20)*ones(size(sv_c2,2),1);
[angc,cidx] = sort([ang_c1 ang_c2 angexist]);
rc = [r1;r2;rcexist];
rcplot = mag2db(rc(cidx));
rcplot(rcplot<-100) = -100;
% inequality constraints
sv_c1 = steervec(pos,ang_c1);
sv_c2 = steervec(pos,ang_c2);
svcineq = [sv_c1 sv_c2];
rineq = [r1;r2];
angineq = [ang_c1 ang_c2];
end