Суррогатный проект оптимизации шестиэлементной антенны Яги-Уда

Этот пример использует:

Этот пример оптимизирует 6-элементную антенну Yagi-Uda как для направленности, так и для 50 $Ω$ входное совпадение с использованием метода глобальной оптимизации, называемого суррогатной оптимизацией. Диаграммы направленности излучения и входное сопротивление антенн чувствительны к параметрам, которые определяют их формы. Многомерная поверхность, над которой должны быть выполнены такие оптимизации, имеет несколько локальных оптимумов. Это делает задачу нахождения правильного набора параметров, удовлетворяющих целям оптимизации, особенно сложной и требует использования методов глобальной оптимизации. Эта антенна Яги-Уда предназначена для операции в составе ретрансляционной станции в радиостанции HAM настройки.

Расчётные параметры

Выберите исходные расчётные параметры в центре полосы УКВ.

fc = 144.5e6;
wirediameter = 12e-3;
c = physconst('lightspeed');
lambda = c/fc;
Z0 = 50;
BW = 0.015*fc;
fmin = fc - 2*(BW);
fmax = fc + 2*(BW);
Nf = 101;
freq = linspace(fmin,fmax,Nf);

Создайте антенну Яги-Уда

Ведомым элементом для этой антенны Яги-Уды является диполь. Это стандартный возбудитель для такой антенны. Настройте параметры длины и ширины диполя. Поскольку мы моделируем цилиндрические структуры как эквивалентные металлические полосы, ширина вычисляется с помощью служебной функции, доступной в Antenna Toolbox™. Длина выбрана приблизительно такой, чтобы $λ / 2$ на проект частоте.

d = dipole;
d.Length = 0.982.*(lambda/2);
d.Width = cylinder2strip(wirediameter/2);
d.Tilt = 90;
d.TiltAxis = 'Y';

Создайте антенну Яги-Уда с возбудителем в качестве диполя. Установите количество директоров равным четырем. Выбор длин элементов и интервалов между элементами является начальным предположением и будет служить начальной точкой для процедуры оптимизации. Показать начальный проект.

Numdirs = 4;
refLength = 0.25; 
dirLength = [0.940 0.910 0.850 0.830]; 
refSpacing = 0.35;
dirSpacing = [0.15 0.2 0.3 0.3 ];
initialdesign = [refLength refSpacing dirSpacing].*lambda;
yagidesign = yagiUda;
yagidesign.Exciter = d;
yagidesign.NumDirectors = Numdirs;
yagidesign.ReflectorLength = refLength; 
yagidesign.DirectorLength = dirLength; 
yagidesign.ReflectorSpacing = refSpacing*lambda;
yagidesign.DirectorSpacing = dirSpacing*lambda;
show(yagidesign)

График диаграммы направленности излучения на расчетной частоте

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

pattern(yagidesign,fc);

Эта исходная антенна Яги-Уда не имеет более высокой направленности в предпочтительном направлении, что означает в зените (повышение = 90 o), и поэтому является плохо спроектированным излучателем.

Настройка оптимизации

Используйте следующие переменные в качестве управляющих переменных для оптимизации:

Длина отражателя (1 переменная)
Интервал между рефлекторами (1 переменная)
Интервалы директора (4 переменные)

В терминах одного векторного параметра controlVals, set

Длина отражателя = controlVals(1)
Интервал между рефлекторами = controlVals(2)
Интервалы между директорами = controlVals(3:6)

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

type yagi_objective_function_surrogate.m

function objectivevalue = yagi_objective_function_surrogate(y,controlVals,fc,BW,ang,Z0,constraints)
% YAGI_OBJECTIVE_FUNCTION returns the objective for a 6 element Yagi
% OBJECTIVE_VALUE =
% YAGI_OBJECTIVE_FUNCTION(Y,CONTROLVALS,FREQ,ANG,Z0,constraints), assigns
% the appropriate parasitic dimensions, CONTROLVALS to the Yagi antenna Y,
% and uses the frequency FREQ, angle pair,ANG, reference impedance Z0 and
% the constraints to calculate the objective function value.

% The YAGI_OBJECTIVE_FUNCTION function is used for an internal example.
% Its behavior may change in subsequent releases, so it should not be
% relied upon for programming purposes.

% Copyright 2018 The MathWorks, Inc.

y.ReflectorLength = controlVals(1);
y.ReflectorSpacing = controlVals(2);
y.DirectorSpacing = controlVals(3:end);


% Unpack constraints
Gmin = constraints.Gmin;
Gdev = constraints.Gdeviation;
FBmin = constraints.FBmin;
S11min = constraints.S11min;
K = constraints.Penalty;

% Calculate antenna port and field parameters
output = analyzeAntenna(y,fc,BW,ang,Z0);

% Form objective function
output1 = output.MaxDirectivity+output.MismatchLoss;    % Directivity/Gain at zenith
Gain1 = output.MaxDirectivity1+output.MismatchLoss1;    % Directivity/Gain at zenith
Gain2 = output.MaxDirectivity2+output.MismatchLoss2;    % Directivity/Gain at zenith
% Gain constraint, e.g. G > 10
c1 = 0;
if output1<Gmin
    c1 = abs(Gmin-output1)+abs(Gain1-Gmin)+abs(Gain2-Gmin);
end

% Gain deviation constraint, abs(G-Gmin)<0.1;
c1_dev = 0;
c1_dev_temp = 0;
c1_dev_temp1 = 0;
c1_dev_temp2 = 0;
if abs(output1-Gmin)>Gdev
    c1_dev_temp = -Gdev + abs(output1-Gmin);
end
if abs(Gain1-Gmin)>Gdev
    c1_dev_temp1 = -Gdev + abs(Gain1-Gmin);
end
if abs(Gain2-Gmin)>Gdev
    c1_dev_temp2 = -Gdev + abs(Gain2-Gmin);
end
c1_dev = (c1_dev_temp+c1_dev_temp1+c1_dev_temp2)/3;

% Front to Back Ratio constraint, e.g. F/B > 15
c2 = 0;
% if output.FB < FBmin
%     c2 = FBmin-output.FB;
% end
c2 = (abs(FBmin-output.FB)+abs(FBmin-output.FB1)+abs(FBmin-output.FB2))/3;

% Reflection Coefficient, S11 < -10
c3 = 0;
if output.S11 > S11min
   c3 = -S11min + output.S11; 
end

% Form the objective + constraints
objectivevalue = -output1 + max(0,(c1+c1_dev+c2+c3))*K;
end

function output = analyzeAntenna(ant,fc,BW,ang,Z0)
%ANALYZEANTENNA calculate the objective function
% OUTPUT = ANALYZEANTENNA(Y,FREQ,BW,ANG,Z0) performs analysis on the
% antenna ANT at the frequency, FC, and calculates the directivity at the
% angles specified by ANG and the front-to-back ratio. The reflection
% coefficient relative to reference impedance Z0, and impedance are
% computed over the bandwidth BW around FC.

fmin = fc - (BW/2);
fmax = fc + (BW/2);
Nf = 5;
freq = unique([fc,linspace(fmin,fmax,Nf)]);
fcIdx = freq==fc;
fcIdx1 = freq==fmin;
fcIdx2 = freq==fmax;
s = sparameters(ant,freq,Z0);
Z = impedance(ant,fc);
Z1 = impedance(ant,fmin);
Z2 = impedance(ant,fmax);
az = ang(1,:);
el = ang(2,:);
Dmax = pattern(ant,fc,az(1),el(1));
Dmax1 = pattern(ant,fmin,az(1),el(1));
Dmax2 = pattern(ant,fmax,az(1),el(1));
Dback = pattern(ant,fc,az(2),el(2));
Dback1 = pattern(ant,fmin,az(2),el(2));
Dback2 = pattern(ant,fmax,az(2),el(2));
% Calculate F/B
F_by_B = Dmax-Dback;
F_by_B1 = (Dmax1-Dback1);
F_by_B2 = (Dmax2-Dback2);
% Compute S11 and mismatch loss
s11 = rfparam(s,1,1);
S11 = max(20*log10(abs(s11)));
T = mean(10*log10(1 - (abs(s11)).^2));
T1 = max(10*log10(1 - (abs(s11(fcIdx1))).^2));
T2 = max(10*log10(1 - (abs(s11(fcIdx2))).^2));
% Form the output structure
output.MaxDirectivity= Dmax;
output.BackLobeLevel = Dback;
output.MaxDirectivity1= Dmax1;
output.BackLobeLevel1 = Dback1;
output.MaxDirectivity2= Dmax2;
output.BackLobeLevel2 = Dback2;
output.FB = F_by_B;
output.FB1 = F_by_B1;
output.FB2 = F_by_B2;
output.S11 = S11;
output.MismatchLoss = T;
output.MismatchLoss1 = T1;
output.MismatchLoss2 = T2;
output.Z = Z;
output.Z1 = Z1;
output.Z2 = Z2;
end

Установите ограничения на переменные управления.

refLengthBounds = [0.1;                         % lower bound on reflector length
                   0.6];                        % upper bound on reflector spacing
dirLengthBounds = [0.3 0.3 0.3 0.3;             % lower bound on director length
                   0.7 0.7 0.7 0.7];            % upper bound on director length
refSpacingBounds = [0.25;                       % lower bound on reflector spacing
                    0.65];                      % upper bound on reflector spacing
dirSpacingBounds = [0.01 0.1 0.1 0.1;           % lower bound on director spacing 0.2 0.25 0.3 0.3
                    0.2 0.25 0.35 0.35];        % upper bound on director spacing
exciterLengthBounds = [0.45;                    % lower bound on exciter length
                       0.6];                    % upper bound on exciter length
exciterSpacingBounds = [.004;
                        .008];
LB = [refLengthBounds(1) refSpacingBounds(1) dirSpacingBounds(1,:) ].*lambda;
UB = [refLengthBounds(2) refSpacingBounds(2) dirSpacingBounds(2,:) ].*lambda;
parameterBounds.LB = LB;
parameterBounds.UB = UB;
ang = [0 0;90 -90];                             % azimuth,elevation angles for main lobe and back lobe [az;el]

Оптимизация на основе суррогатов

Global Optimization Toolbox™ обеспечивает основанную на суррогате оптимизационную функцию, называемую surrogate. Мы используем эту функцию с опциями, заданными в optimoptions функция. На каждой итерации постройте график наилучшего значения целевой функции и ограничьте общее количество итераций 300. Передайте целевую функцию в суррогатную функцию при помощи анонимной функции вместе с границами и структурой опций. Целевая функция, используемая в процессе оптимизации surrogate доступно в файле yagi_objective_function.

% Optimizer options
optimizer = 'Surrogate';

if strcmpi(optimizer,'PatternSearch')
    optimizerparams = optimoptions(@patternsearch);
    optimizerparams.UseCompletePoll = true;
    optimizerparams.PlotFcns = @psplotbestf;
    optimizerparams.UseParallel = true;
    optimizerparams.Cache = 'on';     
    optimizerparams.MaxFunctionEvaluations = 1200; 
    optimizerparams.FunctionTolerance = 1e-2;
elseif strcmpi(optimizer,'Surrogate')
    optimizerparams = optimoptions(@surrogateopt);
    optimizerparams.UseParallel = true;
    optimizerparams.MaxFunctionEvaluations = 600;   
    optimizerparams.MinSurrogatePoints = 12;        
    optimizerparams.InitialPoints = initialdesign;
else
    error('Optimizer not supported');
end

% Antenna design parameters
designparams.Antenna = yagidesign;
designparams.Bounds = parameterBounds;

% Analysis parameters
analysisparams.CenterFrequency = fc;
analysisparams.Bandwidth = BW;
analysisparams.ReferenceImpedance = Z0;
analysisparams.MainLobeDirection = ang(:,1);
analysisparams.BackLobeDirection = ang(:,2);

% Set constraints
constraints.S11min = -15;
constraints.Gmin = 10;
constraints.Gdeviation = 0.1;
constraints.FBmin = 20;
constraints.Penalty = 75;
poolobj = gcp;
optimdesign = optimizeAntennaSurrogate(designparams,analysisparams,constraints,optimizerparams);

Surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

Оптимизированный шаблон графика

Постройте график оптимизированного шаблона антенны на проект частоте.

yagidesign.ReflectorLength = optimdesign(1);
yagidesign.ReflectorSpacing = optimdesign(2);
yagidesign.DirectorSpacing = optimdesign(3:6);
pattern(yagidesign,fc)

Вырезы шаблона E и H-плоскости

Чтобы получить лучшее представление о поведении в двух ортогональных плоскостях, постройте график нормированной величины электрического поля в E и H-плоскостях, то есть азимут = 0 и 90 o соответственно. Включите метрики антенны на графиках полярного шаблона, чтобы установить направленность в зените, соотношение Фронт-Назад и ширину луча в E и H-плоскостях.

fU = fc+ BW/2;
fL = fc-BW/2;
figure;
patternElevation(yagidesign,fc,0,'Elevation',0:1:359);
pE = polarpattern('gco');
DE_fL = patternElevation(yagidesign,fL,0,'Elevation',0:1:359);
DE_fU = patternElevation(yagidesign,fU,0,'Elevation',0:1:359);
add(pE,[DE_fL, DE_fU])
pE.MagnitudeLim = [-20 15];
pE.TitleTop = 'E-plane Directivity (dBi)';
pE.LegendLabels = {[num2str(fc./1e6), ' MHz'],[num2str(fL./1e6), ' MHz'],[num2str(fU./1e6), ' MHz']};

figure;
patternElevation(yagidesign,fc,90,'Elevation',0:1:359);
pH = polarpattern('gco');
DH_fL = patternElevation(yagidesign,fL,90,'Elevation',0:1:359);
DH_fU = pattern(yagidesign,fU,90,'Elevation',0:1:359);
add(pH,[DH_fL, DH_fU])
pH.MagnitudeLim = [-20 15];
pH.TitleTop = 'H-plane Directivity (dBi)';
pH.LegendLabels = {[num2str(fc./1e6), ' MHz'],[num2str(fL./1e6), ' MHz'],[num2str(fU./1e6), ' MHz']};

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

Входной коэффициент отражения оптимизированной антенны

Входной коэффициент отражения для оптимизированной антенны Яги-Уды вычисляется и строится относительно опорного импеданса $50 Ω$ . Значение -10 дБ или ниже рассматривается как хорошее соответствие импедансу.

s = sparameters(yagidesign,freq,Z0);

figure;
rfplot(s);

Табуляция начального и оптимизированного проекта

Сведите в таблицу начальные догадки проекта и окончательные оптимизированные проектные значения.

yagiparam=  {'Reflector Length';'Reflector Spacing';'Director Spacing - 1';
             'Director Spacing - 2';'Director Spacing - 3';
             'Director Spacing - 4'};         
initialdesign = initialdesign';
optimdesign = optimdesign';
Tgeometry = table(initialdesign,optimdesign,'RowNames',yagiparam)

Tgeometry=6×2 table
                            initialdesign    optimdesign
                            _____________    ___________

    Reflector Length           0.51867          1.1523  
    Reflector Spacing          0.72614         0.52421  
    Director Spacing - 1        0.3112         0.17164  
    Director Spacing - 2       0.41494         0.46889  
    Director Spacing - 3       0.62241         0.72614  
    Director Spacing - 4       0.62241         0.69233

Изготовленные антенны

Была изготавливаем оптимизированный проект яги. Настоящий яги нуждается в опорном элементе вдоль продольной оси, чтобы придать механическую жесткость. Этот опорный элемент называется стрелой и обычно изготавливается из неметаллического материала. В этом случае для изготовления стрелы использовалась трубка pVC. Обратите внимание, что эффект этой стрелы не был смоделирован в элементе Antenna Toolbox yagiUda. Другой метод анализа входного соответствия состоит в том, чтобы вычислить и построить график VSWR (коэффициент стоячей волны напряжения). Вычислите и постройте график предсказанного VSWR из оптимизированного проекта. Изготовленные антенны VSWR также измеряли с использованием SWR-счетчика. Эти данные были сохранены в файле CSV. Наложите измеренные результаты на анализ.

vswr_measured = csvread('SWR_Values_Sep_15.csv',1,0);
figure
vswr(yagidesign,freq,Z0)
hold on
plot(vswr_measured(:,1),vswr_measured(:,2),'k-.')
legend('Analysis','Measurement')
title('VSWR comparison - No coaxial cable in analysis')
ylabel('Magnitude')

Моделируйте эффект коаксиального кабеля

Коаксиальный кабель, соединенный с изготовленной антенной яги, является RG-58/U с характерным импедансом 50 $Ω$ . Создайте модель этого коаксиального кабеля с помощью RF Toolbox.

out_radius = 3.51e-3;
in_radius = 0.91e-3;
eps_r = 2.95;
line_length = 5.05*lambda;
coax_cable = rfckt.coaxial;
coax_cable.OuterRadius = out_radius;
coax_cable.InnerRadius = in_radius;
coax_cable.EpsilonR = eps_r;
coax_cable.LossTangent = 2e-4;
coax_cable.LineLength = line_length;

Проанализируйте коаксиальный кабель в области значений частот, предназначенных для операции, и используйте импеданс яги в качестве нагрузки. Вычислите входной VSWR для коаксиального кабеля и антенны яги.

Zyagi =impedance(yagidesign,freq);
analyze(coax_cable,freq,Zyagi);
figure
hline = plot(coax_cable,'VSWRin','None');
hline.LineWidth = 2;
hold on
plot(vswr_measured(:,1),vswr_measured(:,2),'k-.')
legend('Analysis','Measurement')
title('VSWR comparison with coaxial cable model')

Кривая VSWR анализа для комбинации коаксиальной антенны и антенны яги выгодно отличается от измеренных данных.

Результаты оптимизированного проекта выгодно отличаются от изготовленной антенны. Эта антенна будет использоваться в составе ретрансляционной станции, работающей на частоте 145 МГц.

figure
yagidesign.Tilt = 90;
yagidesign.TiltAxis = [0 1 0];
show(yagidesign)

См. также

Оптимизация на основе прямого поиска шестиэлементной антенны Яги-Уды

Документация