Antenna Design Genetic Algorithm

Background

For certain classes of antennas, e.g. Yagi-Uda antennas, the design characteristics have no known “best case” numeric values.  That is, if you want to design a Yagi-Uda antenna for a particular frequency, there is no known numeric solution for the width of the dipoles, number of dipoles, and distance between each dipole in order to achieve the highest gain.  Instead, people rely on experimental evidence: the designs of a number of common frequencies have been tested in the field to produce certain amount of gain, so if you know what frequency you’re looking to design for, you go look up the tables based upon what others have tested for.  If you’re looking for an entirely unique frequency, you have to go experiment yourself.

New Approach

I wrote a MatLab program to use a genetic algorithm to modify the parameters of a antenna and eventually “give birth” to a “best known case” antenna based upon forward gain, etc. Currently, it uses NEC2 (Numerical Electromagnetics Code 2) as the processing engine. It writes out a text file to disk, then calls NEC2 to process that file.  This allows us to try a number of unknown antenna designs and permute possible solutions.  It can be run in a distributed fashion, with each machine “phoning home” to a central database which then redistributes the best-known designs to the worker machines.

The output is something like the following, which shows the forward gain of a given design through a set of frequencies

Output from a forward gain analysis

Code

Genetic Algorithm:

function [beta,stopcode]=ga(funstr,parspace,options,p1,p2,p3,p4,p5,p6,p7,p8,p9)
%[beta,stopcode]=ga(funstr,parspace,options,p1,p2,p3,p4,p5,p6,p7,p8,p9)
% Genetic Algorithm for function maximization.
%
% OUTPUTS:
%  beta       = (1 x K) parameter vector maximizing funstr
%  stopcode   = code for terminating condition
%                == 1 if terminated normally
%                == 2 if maximum number of iterations exceeded
%
% INPUTS:
%  funstr     = name of function to be maximized (string).
%  parspace   = (2 x K) matrix is [min; max] of parameter space dimensions
%               or, if (3 x K), then bottom row is a good starting value
%  options    = vector of option settings
%  p1,p2,...,p9 are optional parameters to be passed to funstr
%
% where:
% options(1) = m (size of generation, must be even integer)
% options(2) = eta (crossover rate in (0,1); use Booker's VCO if < 0)
% options(3) = gamma (mutation rate in (0,1))
% options(4) = printcnt (print status once every printcnt iterations)
%                Set printcnt to zero to suppress printout.
% options(5) = maxiter (maximum number of iterations)
% options(6) = stopiter (minimum number of gains < epsln before stop)
% options(7) = epsln (smallest gain worth recognizing)
% options(8) = rplcbest (every rplcbest iterations, insert best-so-far)
% options(9) = 1 if function is vectorized (i.e., if the function
%                can simultaneously evaluate many parameter vectors).
%    Default option settings: [20,-1,0.12,10,20000,2000,1e-4,50,0]
%
% Note: 
%    The function is maximized with respect to its first parameter,
%    which is expressed as a row vector.
%    Example: 
%      Say we want to maximize function f with respect to vector p,
%      and need also to pass to f data matrices x,y,z.  Then,
%      write the function f so it is called as f(p,x,y,z).  GA will
%      assume that p is a row vector.

defopt=[200,-1,0.12,10,20000,2000,1e-4,50,0];
months = ['Jan';'Feb';'Mar';'Apr';'May';'Jun';...
          'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];

if nargin>2
   if isempty(options)
        options=defopt;
   end
else
   options=defopt;
end
m=options(1); eta=options(2); gam=options(3);
printcnt=options(4);
maxiter=options(5);
stopiter=options(6); epsln=options(7);
rplcbest=options(8);
vecfun=options(9);

% Use Booker's VCO if eta==-1
vco=(eta<0);  eta=abs(eta);

% Cancel rplcbest if <=0
if rplcbest<=0, rplcbest=maxiter+1; end

K=size(parspace,2);

% Draw initial Generation
G=rand(m,K).*(parspace(2*ones(m,1),:)-parspace(ones(m,1),:))...
       +parspace(ones(m,1),:);
b0rows=size(parspace,1)-2;
if b0rows>0
  G(1:b0rows,:)=parspace(3:b0rows+2,:);
  parspace=parspace([1 2],:);
end

% Initial 'best' holders
inarow=0;
bestfun=-Inf; beta=zeros(1,K);

% Score for each of m vectors
f=zeros(m,1);

% Setup function string for evaluations
paramstr=',p1,p2,p3,p4,p5,p6,p7,p8,p9';
evalstr = [funstr,'(G'];
if ~vecfun
        evalstr=[evalstr, '(i,:)'];
end
if nargin>3, evalstr=[evalstr,paramstr(1:3*(nargin-3))]; end
evalstr = [evalstr, ')'];

% Print header
if printcnt>0
   disp(['Maximization of function ',funstr])
   disp('i      = Current generation')
   disp('best_i = Best function value in generation i')
   disp('best   = Best function value so far')
   disp('miss   = Number of generations since last hit')
   disp('psi    = Proportion of unique genomes in generation')
   disp(sprintf(['\n',blanks(20),'i     best_i        best     miss   psi']))
end

iter=0;  stopcode=0;
oldpsi=1;  % for VCO option
while stopcode==0
   iter=iter+1;
   % Call function for each vector in G
   if vecfun
        f=eval(evalstr);
   else
     for i=1:m
        f(i)=eval(evalstr);
     end
   end
   f0=f;
   [bf0,bx]=max(f);
   bf=max([bf0 bestfun]);
   fgain=(bf-bestfun);
   if fgain>epsln
        inarow=0;
   else
        inarow=inarow+1;
   end
   if fgain>0
        bestfun=bf;
        beta=G(bx(1),:);
   end
   if printcnt>0 & rem(iter,printcnt)==1
        psi=length(unique(G))/m;
        ck=clock;
        ckhr=int2str(ck(4)+100);  ckday=int2str(ck(3)+100);
        ckmin=int2str(ck(5)+100); cksec=int2str(ck(6)+100);
        timestamp=[ckday(2:3),months(ck(2),:),' ',...
           ckhr(2:3),':',ckmin(2:3),':',cksec(2:3),' '];
        disp([timestamp,sprintf('%6.0f %8.5e %8.5e %5.0f %5.3f',...
                [iter bf0 bestfun inarow psi])])
        disp(beta)
        save gabest beta timestamp iter funstr
   end
   % Reproduction
   f=(f-min(f)).^(1+log(iter)/100);
   pcum=cumsum(f)/sum(f);
   r=rand(1,m); r=sum(r(ones(m,1),:)>pcum(:,ones(1,m)))+1;
   G=G(r,:);
   % Crossover
   if vco
        psi=length(unique(G))/m;
        eta=max([0.2 min([1,eta-psi+oldpsi])]);
        oldpsi=psi;
   end   
   y=sum(rand(m/2,1)<eta);
   if y>0
     % choose crossover point
     x=floor(rand(y,1)*(K-1))+1;
     for i=1:y
        tmp=G(i,x(i)+1:K);
        G(i,x(i)+1:K)=G(i+m/2,x(i)+1:K);
        G(i+m/2,x(i)+1:K)=tmp;
     end
   end
   % Mutation
   M=rand(m,K).*(parspace(2*ones(m,1),:)-parspace(ones(m,1),:))...
       +parspace(ones(m,1),:);
   domuta=find(rand(m,K)<gam);
   G(domuta)=M(domuta);
   % Once every rplcbest iterations, re-insert best beta
   if rem(iter,rplcbest)==0
        G(m,:)=beta;
   end
   stopcode=(inarow>stopiter)+2*(iter>maxiter);
end

if printcnt>0
   if stopcode==1
        disp(sprintf('GA: No improvement in %5.0f generations.\n',stopiter))
   else
        disp(sprintf('GA: Maximum number of iterations exceeded.\n'))
   end
end
% end of GA.M
function [gain_t]=Yagi(p)
% This program is used with a getic optimization code.
% It creates an NEC input file given the parameters for a 3 element YAGI
% antenna.  Then, it runs NEC and reads in the parameters(Gain, Impedance, etc) generated by NEC.
%=========Create NEC input file========================================
Fname_nec='Yagi.nec';
FID_nec=fopen(Fname_nec,'wt');
D=0.0085; % diamter of elements in wavelengths
R=D/2;
Lr=p(1); %perimeter of reflector
Ls=p(2); %perimeter of driven element
Ld=p(3); %perimeter of director
Sr=-p(4); %location of reflector
Sd=p(5);  %location of director
%    Geometry input for NEC

rsl = Lr / 4; %reflector single side length: square loop, one side = permiter / 4
ssl = Ls / 4; %driven single side length: square loop, one side = permiter / 4
dsl = Lr / 4; %director single side length: square loop, one side = permiter / 4
num_segments_per_side = 7;

fprintf(FID_nec,strcat('CM UDA-YAGI SQUARE LOOP ANTENNA','\n'));
fprintf(FID_nec,strcat('CE FIle Generated by MatLab','\n'));
fprintf(FID_nec,'GW %3i %3i %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n',1,num_segments_per_side,Sr,0,-Lr/2,Sr,0,Lr/2,R); %Reflector
fprintf(FID_nec,'GW %3i %3i %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n',2,num_segments_per_side,0 ,0,-Ls/2, 0,0,Ls/2,R); %Driven Element
fprintf(FID_nec,'GW %3i %3i %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n',3,num_segments_per_side,Sd,0,-Ld/2,Sd,0,Ld/2,R); %Director
%    Program Control Commands for NEC
fprintf(FID_nec,strcat('GE','\n'));
fprintf(FID_nec,'EX %3i %3i %3i %3i %8.4f %8.4f\n',0,2,4,0,1,0);%Exitation Command wire 2 segment 4
fprintf(FID_nec,'FR %3i %3i %3i %3i %8.4f %8.4f\n',0,1,0,0,2400,0);%set freq to 299.8 MHz so wavelength will be 1m
fprintf(FID_nec,'RP %3i %3i %3i %3i %8.4f %8.4f %8.4f %8.4f\n',0,1,1,1000,90,0,0,0);%calculate gain at boresite
fprintf(FID_nec,strcat('EN','\n'));
fclose(FID_nec);
%=======Create file to pipe to NEC ===================================
FID_input=fopen('input_CMD','wt');
fprintf(FID_input,strcat(Fname_nec,'\n'));
fprintf(FID_input,strcat('NEC.out','\n'));
fclose(FID_input);
%=======Run NEC======================================================
!NEC2Dx500 < input_CMD >tmp;
%=======Read Data form NEC output file===============================
[freq,Z,gain_t,E_theta,E_phi,n_freq_meas,run_time] = nec_read('NEC.out');