Вторая часть синтеза диаграммы направленности антенной решетки: Оптимизация

В этом примере показано, как использовать методы оптимизации, чтобы выполнить синтез шаблона. Используя минимальное отклонение beamforming как рисунок, этот пример покрывает, как настроить синтез шаблона как задачу оптимизации и затем решить его с помощью решателей оптимизации.

Синтез шаблона может быть достигнут с помощью множества методов. Первая часть Синтеза Диаграммы направленности антенной решетки: Обнуление, Работа с окнами и Утончение примера ввели методы обнуления и работы с окнами.

Этот пример требует Optimization Toolbox™.

Минимальное отклонение Beamforming

По сравнению с одним элементом главное преимущество фазированной решетки является возможностью формирования лучей. В общем случае фазированная решетка N-элемента обеспечивает степени свободы N при синтезировании шаблона, означать, что можно настроить веса N, один для каждого элемента, чтобы управлять формой луча, чтобы удовлетворить некоторым предопределенным ограничениям.

Много популярных beamforming методов могут быть описаны как задача оптимизации. Например, популярный формирователь луча минимального отклонения ответа без искажений (MVDR) используется, чтобы минимизировать общий шумовой выход при сохранении сигнала от данного направления. Математически, MVDR beamforming веса может быть найден путем решения следующей задачи оптимизации

wo=minwwHRw            s.t.B(θ0,w)=1

где R ковариационная матрица, B(θ0,w) ответ луча к направлению θ0 с вектором весов w, и wo вектор, содержащий оптимальные веса. Поскольку диаграмма направленности может быть записана как

B(θ,w)=wHv(θ)

где v(θ) держащийся вектор, соответствующий углу θ, ограничения имеют линейное отношение в комплексной области.

Ограничение без искажений в θ0 только берет 1 степень свободы на пробел проекта. Поскольку массив N-элемента может обработать N-1 ограничения, формирователь луча MVDR может быть расширен, чтобы включать больше ограничений, который приводит к формирователю луча линейного ограничительного отклонения минимума (LCMV). Дополнительные ограничения в формирователе луча LCMV часто используются к пустому указателю интерференция от данных направлений. Поэтому такой формирователь луча LCMV является решением следующей задачи оптимизации

wo=minwwHRw            s.t.B(θ0,w)=1                    B(θk,w)=0      k=1,2,,K

где θk представляет kth обнуление направления и существуют K обнуление направлений всего. Поскольку все ограничения линейны в комплексной области, LCMV beamforing веса может быть выведен аналитически с помощью метода множителя Лагранжа.

Затем ULA с 32 элементами с половиной интервала длины волны используется, чтобы продемонстрировать, как выполнить синтез массивов с помощью подхода LCMV. Обратите внимание на то, что даже при том, что этот пример использует линейную матрицу, метод применяется к любой геометрии массивов. Примите, что сигнал интереса в 0 азимутах степеней, и интерференция в-70,-40, и-20 азимутов степеней. К шумовой власти приходят, чтобы быть в на 40 дБ ниже сигнала. Закрытое решение для формы для весов LCMV может быть вычислено с помощью lcmvweights функция.

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)),'b');
hold on;
plot([ang_i;ang_i],repmat([-100;0],1,size(ang_i,2)),'k--','LineWidth',1.5);
hold off;
ylim([-100 0]); grid on;
legend(plcmv,'LCMV - analytical');
xlabel('Azimuth Angle (deg)');
ylabel('Beam Pattern (dB)');

Figure contains an axes object. The axes object contains 4 objects of type line. This object represents LCMV - analytical.

Рисунок показывает, что вычисленные веса удовлетворяют всем ограничениям.

LCMV Beamforming Используя квадратичное программирование

С тех пор LCMV beamforming может быть описан как задача оптимизации, было бы полезно, чтобы видеть, как задача может быть решена с помощью методов оптимизации. Когда закрытое решение для формы не доступно, например, когда ограничения включают неравенства, такой рабочий процесс становится важным.

При использовании методов оптимизации, чтобы выполнить задачи синтеза шаблона, первый шаг должен выбрать решатель. Поскольку решатели часто разрабатываются, чтобы обратиться к конкретной формулировке задачи, она имеет наибольшее значение, чтобы сделать правильный выбор.

Решатель квадратичного программирования (QP) часто используется, чтобы решить квадратичную целевую функцию с линейными ограничениями. Формулировкой дают

wo=minw12wTHw+fTw            s.t.Awb                     Aeqw=beq                     wlbwwub

где A и Aeq матрицы, задающие ограничения неравенства и ограничения равенства, соответственно; b и beq векторы, задающие желаемые границы для ограничений неравенства и ограничений равенства, соответственно; и wlb и w укажите на нижние и верхние границы для элементов в w, соответственно.

На первый взгляд формулировка задачи квадратичного программирования, кажется, идеальная пара для LCM beamforming, поскольку мы знаем, что диаграмма направленности может быть записана как B(θ,w)=wHv(θ). К сожалению, это не верно, потому что несмотря на то, что диаграмма направленности линейна в комплексной области, это не находится так в действительной области. С другой стороны, большинство решателей оптимизации, включая квадратичное программирование, работает в действительной области, когда пробел решения становится действительно сложным в комплексной области. Кроме того, неравенство не четко определено в комплексной области также. Таким образом вопрос состоит в том, если LCMV beamforming проблема может быть сформулирован в действительную область и все еще получить форму квадратичного программирования.

Начните с целевой функции wHRw. Обратите внимание на то, что R ковариационная матрица в этом контексте, таким образом wHRw isreal. Поэтому цель состоит в том, чтобы действительно минимизировать действительную часть wHRw, который может быть написан в следующей форме:

{wHRw}=[{w}T{w}T][{R}-{R}{R}{R}][{w}{w}]=wTRw

где w=[{w}{w}] и R=[{R}-{R}{R}{R}].

Затем посмотрите на ограничения. То, что ограничения в LCMV beamforming являются равенствами, делает преобразование немного проще, поскольку одно комплексное линейное ограничение равенства может только быть разделено на два соответствующих действительных линейных ограничения равенства, один для действительной части и другого для мнимой части. Например, следующее ограничение

wHv=rH

может быть переписан как

vHw=r[{v}{v}{v}-{v}][{w}{w}]=[{v}{v}{v}-{v}]w=[{r}{r}]

С этими отношениями квадратичное программирование может затем использоваться, чтобы получить LCMV beamforming веса.

% 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;

% Plot
hold on;
plcmvq = plot(ang_plot,mag2db(abs(w_opt'*sv_plot)),'r--','LineWidth',2);
hold off;
legend([plcmv plcmvq],{'LCMV - Lagrange','LCMV - QP'})

Figure contains an axes object. The axes object contains 5 objects of type line. These objects represent LCMV - Lagrange, LCMV - QP.

Рисунок показывает, что метод QP может получить тот же LCMV beamforming веса.

Минимальное отклонение Beamforming Используя квадратичное программирование

Другое общее требование для синтеза шаблона должно гарантировать, что ответ ниже определенного порога для набора данных углов. Это требования неравенства, и полная задача оптимизации становится

wo=minwwHRw            s.t.B(θ0,w)=1                    B(θk,w)=0      k=1,2,,K                   |B(w,θl)|rl      l=1,2,,L

где θl lугол th, под которым уровень бокового лепестка должен быть меньше порога rl. Существуют L требования неравенства всего.

Как пример, добавьте следующие требования в нашу демонстрационную проблему:

  1. Боковой лепесток должен быть ниже-40 дБ между-30 и-10 градусами в области азимута.

  2. Везде за пределами основного лепестка, боковой лепесток должен быть ниже-20 дБ.

Следующий код строит маску желаемого шаблона. Это также возвращает ограничения неравенства.

[rc,angc,sv_cineq,r_ineq]=generatePatternMask(pos,r_c,[ang_d ang_i]);

cla;
pdp = plot(angc,rc,'k--','LineWidth',1.5);
ylim([-100 0]);
grid on;
xlabel('Azimuth Angle (deg)');
ylabel('Beam Pattern (dB)');

Figure contains an axes object. The axes object contains an object of type line.

С введением нового ограничения неравенства проблема больше не имеет решение закрытой формы. Кроме того, из-за присутствия операции нормы, вышеупомянутый прием, чтобы перевести комплексные ограничения в действительную область больше не работает. Поэтому новый способ перевести комплексное неравенство необходим для того, чтобы использовать решатель QP для этого набора ограничений,

Учитывая ограничение |B(θ,w)|<r, один возможный способ преобразовать это состоит в том, чтобы принять, что и действительная и мнимая часть способствует одинаково норме. В этом случае одно комплексное ограничение неравенства может быть разделено на два действительных неравенства как

-r2{B(θ,w)}r2-r2{B(θ,w)}r2

т.е. .

{B(θ,w)}r2{B(θ,w)}r2-{B(θ,w)}r2-{B(θ,w)}r2

Действительная и мнимая часть ограничений может быть извлечена как описано в предыдущем разделе. Следовательно, следующие кодовые наборы ограничения и решают задачу синтеза шаблона.


% 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;

% Plot
hold on;
pqp=plot(ang_plot,mag2db(abs(w_opt'*sv_plot)),'b','LineWidth',2);
hold off;
legend(pqp,'Pattern - QP')

Figure contains an axes object. The axes object contains 2 objects of type line. This object represents Pattern - QP.

Получившийся шаблон удовлетворяет требованиям.

Минимальное отклонение Beamforming Используя программирование конуса второго порядка

Даже при том, что решатель QP решает минимальное отклонение beamforming проблема успешно, преобразование ограничений неравенства от комплекса до действительного не идеально. Например, нет никакой причины, что действительные и мнимые части должны всегда одинаково способствовать норме. Однако нет действительно никакого простого способа обработать ограничение с нормой в настройке задач QP.

К счастью существует другой решатель оптимизации, который может использоваться в этом случае, названный Программированием конуса второго порядка (SOCP). SOCP разрабатывается, чтобы эффективно решить следующую задачу оптимизации:

wo=minwfTw            s.t.Ascw-bscdTw-γ                     Awb                     Aeqw=beq                     wlbwwub

где Asc матрица, bscи d векторы, и γ скаляр. Asc, bsc, d, и γвместе задайте ограничение конуса второго порядка, которое является естественным выбором задать ограничение на нормы. Решатель может обработать несколько ограничений конуса второго порядка одновременно.

По сравнению с формулировкой задачи для QP SOCP может не только обработать ограничения второго порядка, но также и обработать все ограничения, заданные в QP. Однако целевая функция очень отличается. Вместо того, чтобы минимизировать отклонение, SOCP минимизируют линейную функцию переменной оптимизации, w. К счастью, существует способ преобразовать минимальную целевую функцию отклонения, чтобы поместиться в эту среду SOCP. Прием должен задать вспомогательную переменную, t, и гарантируйте, что отклонение ограничено t. Таким образом, проблема становится, чтобы минимизировать t в то время как исходная квадратичная цель становится коническим ограничением, ограниченным t.

После этой мысли, исправления t к нижней части w. Теперь пространство поиска оптимизации задано w=[wt]. Поскольку цель состоит в том, чтобы минимизировать t, f задан как fT=[0,,01].

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

Затем исследуйте, как преобразовать квадратичное ограничение в ограничение конуса второго порядка. Поскольку идея является к ограничению отклонением как

wTRwt

это эквивалентно, чтобы сказать

R0.5w2t+14

Таким образом,

[R0.5001][wt]2R0.5w2+t2t2+t+14=(t+12)2

и второе коническое ограничение становится

[R0.5001]wt+12

% 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);

Остающиеся вторые конические ограничения от исходных ограничений неравенства и являются довольно прямыми, чтобы быть преобразованными. Каждое ограничение неравенства в форме

B(θ,w)r

может быть writtedn как ограничение конуса второго порядка как

[{B(θ,w)}0{B(θ,w)}000]wr

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

Ограничения равенства также должны быть увеличены, чтобы вместить вспомогательную переменную, t.

% patch Aeq and beq
Aeqc = [Aeq zeros(size(Aeq,1),1);zeros(1,2*N) 0]; 
beqc = [beq;0];

Веса массивов могут теперь быть вычислены с помощью решателя SOCP. Поскольку ограничения неравенства являются более естественными, чтобы быть представленными в конических ограничениях, наивно можно было бы ожидать, что в этом примере, получившийся шаблон от SOCP будет лучше по сравнению с шаблоном от 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;
pqp.LineStyle = '-.';
psocp = plot(ang_plot,mag2db(abs(w_op'*sv_plot)),'r-','LineWidth',2);
hold off;
legend([pqp psocp],{'Pattern - QP','Pattern - SOCP'})

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Pattern - QP, Pattern - SOCP.

Действительно, график показывает, что шаблон, полученный через SOCP, имеет лучшее подавление бокового лепестка по сравнению с шаблоном, полученным с помощью QP.

Сводные данные

В этом примере показано, как использовать методы оптимизации, чтобы выполнить синтез шаблона. Пример начат со знакомой проблемы LCMV и расширенный желаемый шаблон, чтобы включать ограничения неравенства. По сравнению с решением закрытой формы методы оптимизации более гибки и могут быть применены, чтобы обратиться к более комплексным ограничениям. Пример также показал, как обратиться к ограничениям неравенства с помощью или квадратичного программирования или программирования конуса второго порядка.

Ссылка

[1] Херв Лебрет и Стивен Бойд, Синтез Шаблона Антенной решетки через Выпуклую Оптимизацию, Сделку IEEE на Обработке сигналов, Издании 45, № 3, март 1997

function [rcplot,angc,svcineq,rineq]=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];

end