% Matlab script for the paper "Beamforming with Uncertain Weights", % by A. Mutapcic, S.-J. Kim, and S. Boyd. % % We consider beamformers with complex multiplicative uncertainty % in the element weights, and magnitude bound on each uncertainty. % In this case, the robust optimal beamformers can be designed using % the weighted complex l_1-regularization of the nominal problem. % % Nominal problem: minimize the rejection level over the rejection % band, subject to unit array gain in the desired % direction theta_des. % % minimize max |y(theta)| for theta in rejection band % s.t. Re y(theta_des) >= 1 % % This script requires CVX software, which is available at % http://www.stanford.edu/~boyd/cvx/ %******************************************************************** % antenna specs for our problem instance %******************************************************************** lambda = 1; % wavelength theta_des = 60; % desired direction (should be an integer -- discretization) half_beamwidth = 20; % half beamwidth around the target direction % rectangular 6 x 6 grid of 36 sensor (antenna) elements m = 6; n = m^2; d = 0.45*lambda; % sensor locations loc = zeros(n,2); for x = 0:m-1 for y = 0:m-1 loc(m*y+x+1,:) = [x y]; end end loc = loc*d; % angle sampling theta = [1:1:360]'; %******************************************************************** % construct optimization data %******************************************************************** % build matrix A that relates weights w and composite array output y = A*w A = kron(cos(pi*theta/180), loc(:,1)') + kron(sin(pi*theta/180), loc(:,2)'); A = exp(2*pi*i/lambda*A); % desired direction constraint matrix [diff_closest, ind_closest] = min( abs(theta - theta_des) ); Ades = A(ind_closest,:); thdes = theta(ind_closest); % rejection band constraint matrix ind = find(theta <= (theta_des-half_beamwidth) | ... theta >= (theta_des+half_beamwidth) ); Astop = A(ind,:); thstop = theta(ind); %******************************************************************** % nominal beamforming problem %******************************************************************** cvx_begin variable w_nom(n) complex minimize( max( abs( Astop*w_nom ) ) ) subject to real( Ades*w_nom ) >= 1; cvx_end % check if problem was successfully solved disp(['Problem is ' cvx_status]) if ~strcmp(cvx_status,'Solved') return end % compute the optimal rejection level mrl_nom = max(abs(Astop*w_nom)); fprintf(1,'The minimum rejection level is %3.2f dB.\n\n', 20*log10(mrl_nom) ); %******************************************************************** % robust beamforming problem %******************************************************************** % regularization constants (they are unity in this example) mu = max( abs( Astop ) ); nu = abs( Ades ); % loop over different rhos and recompute robust optimal solution rho = linspace(0,0.15,40); mrl_rob = []; mrl_tik_wc = []; mrl_nom_wc = []; cvx_quiet(true); for k = 1:length(rho) % optimal robust solution cvx_begin variable w_rob(n) complex minimize( max( abs( Astop*w_rob ) ) + rho(k)*mu*abs(w_rob) ) subject to real( Ades*w_rob ) >= 1 + rho(k)*nu*abs(w_rob); cvx_end % check if problem was successfully solved disp(['rho = ' num2str(rho(k))]) disp([' Optimal worst-case robust problem is ' cvx_status '.']) if ~strcmp(cvx_status,'Solved') return end % worst-case robust performance mrl_rob(k) = 20*log10( max(abs(Astop*w_rob)) + rho(k)*mu*abs(w_rob) ); % Tikhonov l_2-regularized robust solution % we have found that MU = 2 gives excellent performance (experimental) MU = 2; cvx_begin variable w_tik(n) complex minimize( max( abs( Astop*w_tik ) ) + MU*rho(k)*sum_square_abs(w_tik) ) subject to % note that the solution will not satisfy this constraint robustly real( Ades*w_tik ) >= 1; cvx_end % check if problem was successfully solved disp([' Tikhonov l_2-regularized problem is ' cvx_status '.']) if ~strcmp(cvx_status,'Solved') return end % scale Tikhonov solution so it satisfies real( Ades*w_tik ) == 1 w_tik = w_tik / ( real(Ades*w_tik) - rho(k)*nu*abs(w_tik) ); % worst-case Tikhonov performance mrl_tik_wc(k) = 20*log10( max(abs(Astop*w_tik)) + rho(k)*mu*abs(w_tik) ); % also compute the worst-case performance of the nominal design % w_nom = w_nom / ( real(Ades*w_nom) - rho(k)*nu*abs(w_nom) ); mrl_nom_wc(k) = 20*log10( max(abs(Astop*w_nom)) + rho(k)*mu*abs(w_nom) ); end cvx_quiet(false); %******************************************************************** % plots %******************************************************************** figure(1), clf plot(rho,mrl_nom_wc,'r--','LineWidth',1.15), hold on, plot(rho,mrl_tik_wc,'b--','LineWidth',1.15), hold on, plot(rho,20*log10(mrl_nom + rho./(1-rho)),'r:','LineWidth',1.15), hold on, plot(rho,mrl_rob,'-','Color',[0 .5 0],'LineWidth',1.15), hold off xlabel('rho'), ylabel('mrl') axis([0.01 .15 -25 -10]) % build matrix A for plotting A = kron(cos(pi*theta/180), loc(:,1)') + kron(sin(pi*theta/180), loc(:,2)'); A = exp(2*pi*i/lambda*A); % beamformer plots figure(2), clf Rho = rho(end); deg = 20; delta_rand = (2*Rho*rand(n,1)-1 + i*(2*Rho*rand(n,1)-1)); delta_nom = Rho*exp(i*angle(Astop(deg,:)*w_nom)/n)*... (Astop(deg,:)'.*w_nom)./abs(Astop(deg,:)'.*w_nom); delta_rob = Rho*exp(i*angle(Astop(deg,:)*w_rob)/n)*... (Astop(deg,:)'.*w_rob)./abs(Astop(deg,:)'.*w_rob); delta_tik = Rho*exp(i*angle(Astop(deg,:)*w_tik)/n)*... (Astop(deg,:)'.*w_tik)./abs(Astop(deg,:)'.*w_tik); y_nom = abs( A*w_nom ); y_nom_wc = abs( A*(w_nom.*(1 + delta_nom)) ); y_rob_wc = abs( A*(w_rob.*(1 + delta_rob)) ); y_tik_wc = abs( A*(w_tik.*(1 + delta_tik)) ); dBYnom = 20*log10(y_nom/y_nom(thdes)); dBYnom_wc = 20*log10(y_nom_wc/y_nom_wc(thdes)); dBYrob_wc = 20*log10(y_rob_wc/y_rob_wc(thdes)); dBYtik_wc = 20*log10(y_tik_wc/y_tik_wc(thdes)); plot(theta,dBYnom_wc,'r--','LineWidth',1.15), hold on % plot(theta,dBYtik_wc,'b:','LineWidth',1.15), hold on plot(theta,dBYrob_wc,'-','Color',[0 .5 0],'LineWidth',1.15), hold on ymin = -60; ymax = 20; plot([thdes+half_beamwidth thdes+half_beamwidth],[ymin ymax],'k:',... [thdes thdes],[ymin ymax],'k:',... [thdes-half_beamwidth thdes-half_beamwidth],[ymin ymax],'k:') hold off axis([0 360 ymin ymax]) xlabel('angle'), ylabel('abs pattern')