Максимизация монохроматических поляризованных легких интерференционных шаблонов Используя GlobalSearch и MultiStart

Открыть скрипт

В этом примере показано, как использовать функции GlobalSearch и MultiStart.

Введение

В этом примере показано, как Global Optimization Toolbox функционирует, особенно GlobalSearch и MultiStart, может помочь определить местоположение максимума электромагнитного интерференционного шаблона. Для простоты моделирования шаблон является результатом монохроматического поляризованного легкого распространения из точечных источников.

Электрическое поле, должное получать, i измерилась в направлении поляризации в точке x, и время t

$$E_i = \frac{A_i}{d_i(x)} \sin(\phi_i + \omega (t - d_i(x)/c) ),$$

где$\phi_i$ фаза в начальный момент времени для источника$i$,$c$ является скоростью света,$\omega$ является частотой света,$A_i$ является амплитудой источника $i$и$d_i(x)$ является расстоянием от источника$i$ до$x$.

Для фиксированной точки$x$ интенсивность света является средним во времени квадрата сетевого электрического поля. Сетевое электрическое поле является суммой электрических полей из-за всех источников. Среднее во времени зависит только от размеров и относительных фаз электрических полей в$x$. Чтобы вычислить сетевое электрическое поле, сложите отдельные вклады с помощью метода фазовращателя. Для фазовращателей каждый источник вносит вектор. Длина вектора является амплитудой, разделенной на расстояние от источника и угол вектора,$\phi_i - \omega d_i(x)/c$ фаза в точке.

В данном примере мы задаем три точечных источника с той же частотой ($\omega$) и амплитуда ($A$), но варьировались начальная фаза ($\phi_i$). Мы располагаем эти источники на фиксированной плоскости.

% Frequency is proportional to the number of peaks
relFreqConst = 2*pi*2.5;
amp = 2.2;
phase = -[0; 0.54; 2.07];

numSources = 3;
height = 3;

% All point sources are aligned at [x_i,y_i,z]
xcoords = [2.4112
           0.2064
           1.6787];
ycoords = [0.3957
           0.3927
           0.9877];
zcoords = height*ones(numSources,1);

origins = [xcoords ycoords zcoords];

Визуализируйте интерференционный шаблон

Теперь давайте визуализировать срез интерференционного шаблона на плоскости z = 0.

Как вы видите из графика ниже, существует много peaks и оврагов, указывающих на конструктивную и разрушительную интерференцию.

% Pass additional parameters via an anonymous function:
waveIntensity_x = @(x) waveIntensity(x,amp,phase, ...
    relFreqConst,numSources,origins);
% Generate the grid
[X,Y] = meshgrid(-4:0.035:4,-4:0.035:4);
% Compute the intensity over the grid
Z = arrayfun(@(x,y) waveIntensity_x([x y]),X,Y);
% Plot the surface and the contours
figure
surf(X,Y,Z,'EdgeColor','none')
xlabel('x')
ylabel('y')
zlabel('intensity')

Изложение задачи оптимизации

Мы интересуемся местоположением, где эта интенсивность волны достигает своего самого высокого пика.

Интенсивность волны ($I$) уменьшается, когда мы переезжаем от источника, пропорционального$1/d_i(x)$. Поэтому давайте ограничим пробел эффективных решений путем добавления ограничений в проблему.

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

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

Теперь наша проблема стала:

$$ \max_{x,y} I(x,y) $$

при ограничениях

$$ (x - x_{c1})^2 + (y - y_{c1})^2 \le r_1^2 $$

$$ (x - x_{c2})^2 + (y - y_{c2})^2 \le r_2^2 $$

$$ (x - x_{c3})^2 + (y - y_{c3})^2 \le r_3^2 $$

$$-0.5 \leq x \leq 3.5$$

$$-2 \leq y \leq 3$$

где$(x_{cn},y_{cn})$ и$r_n$ координаты и апертурный радиус$n^{th}$ точечного источника, соответственно. Каждому источнику дают апертуру с радиусом 3. Данные границы охватывают выполнимую область.

Цель ($I(x,y)$) и нелинейные ограничительные функции задана в отдельных файлах MATLAB®, waveIntensity.m и apertureConstraint.m, соответственно, которые перечислены в конце этого примера.

Визуализация с ограничениями

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

% Visualize the contours of our interference surface
domain = [-3 5.5 -4 5];
figure;
ezcontour(@(X,Y) arrayfun(@(x,y) waveIntensity_x([x y]),X,Y),domain,150);
hold on

% Plot constraints
g1 = @(x,y)  (x-xcoords(1)).^2 + (y-ycoords(1)).^2 - 9;
g2 = @(x,y)  (x-xcoords(2)).^2 + (y-ycoords(2)).^2 - 9;
g3 = @(x,y)  (x-xcoords(3)).^2 + (y-ycoords(3)).^2 - 9;
h1 = ezplot(g1,domain);
h1.Color = [0.8 0.7 0.1]; % yellow
h1.LineWidth = 1.5;
h2 = ezplot(g2,domain);
h2.Color = [0.3 0.7 0.5]; % green
h2.LineWidth = 1.5;
h3 = ezplot(g3,domain);
h3.Color = [0.4 0.4 0.6]; % blue
h3.LineWidth = 1.5;

% Plot bounds
lb = [-0.5 -2];
ub = [3.5 3];
line([lb(1) lb(1)],[lb(2) ub(2)],'LineStyle','--')
line([ub(1) ub(1)],[lb(2) ub(2)],'LineStyle','--')
line([lb(1) ub(1)],[lb(2) lb(2)],'LineStyle','--')
line([lb(1) ub(1)],[ub(2) ub(2)],'LineStyle','--')
title('Pattern Contours with Constraint Boundaries')

Подготовка и решение задачи с локальным решателем

Учитывая нелинейные ограничения, нам нужны ограниченный нелинейный решатель, а именно, fmincon.

Давайте настроим структуру задачи, описывающую нашу задачу оптимизации. Мы хотим максимизировать функцию интенсивности, таким образом, мы инвертируем возвращенные значения, формируют waveIntensity. Давайте выберем произвольную стартовую точку, которая, оказывается, около выполнимой области.

Для этой небольшой проблемы мы будем использовать fminconалгоритм SQP.

% Pass additional parameters via an anonymous function:
apertureConstraint_x = @(x) apertureConstraint(x,xcoords,ycoords);

% Set up fmincon's options
x0 = [3 -1];
opts = optimoptions('fmincon','Algorithm','sqp');
problem = createOptimProblem('fmincon','objective', ...
    @(x) -waveIntensity_x(x),'x0',x0,'lb',lb,'ub',ub, ...
    'nonlcon',apertureConstraint_x,'options',opts);

% Call fmincon
[xlocal,fvallocal] = fmincon(problem)

Local 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.


xlocal =

   -0.5000    0.4945


fvallocal =

   -1.4438

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

[~,maxIdx] = max(Z(:));
xmax = [X(maxIdx),Y(maxIdx)]
figure
contour(X,Y,Z)
hold on

% Show bounds
line([lb(1) lb(1)],[lb(2) ub(2)],'LineStyle','--')

% Create textarrow showing the location of xlocal
annotation('textarrow',[0.25 0.21],[0.86 0.60],'TextEdgeColor',[0 0 0],...
    'TextBackgroundColor',[1 1 1],'FontSize',11,'String',{'Single Run Result'});

% Create textarrow showing the location of xglobal
annotation('textarrow',[0.44 0.50],[0.63 0.58],'TextEdgeColor',[0 0 0],...
    'TextBackgroundColor',[1 1 1],'FontSize',12,'String',{'Global Max'});

axis([-1 3.75 -3 3])

xmax =

    1.2500    0.4450

Используя GlobalSearch и MultiStart

Высказанный произвольное исходное предположение, fmincon застревает в соседнем локальном максимуме. Решатели Global Optimization Toolbox, особенно GlobalSearch и MultiStart, дайте нам лучший шанс при нахождении глобального максимума, поскольку они попробуют fmincon от нескольких сгенерированных начальных точек (или наших собственных точек, если мы выбираем).

Наша проблема была уже настроена в problem структура, поэтому теперь мы создаем наши объекты решателя и запускаем их. Первый выход от run местоположение лучшего найденного результата.

% Construct a GlobalSearch object
gs = GlobalSearch;
% Construct a MultiStart object based on our GlobalSearch attributes
ms = MultiStart;

rng(4,'twister') % for reproducibility

% Run GlobalSearch
tic;
[xgs,~,~,~,solsgs] = run(gs,problem);
toc
xgs

% Run MultiStart with 15 randomly generated points
tic;
[xms,~,~,~,solsms] = run(ms,problem,15);
toc
xms

GlobalSearch stopped because it analyzed all the trial points.

All 14 local solver runs converged with a positive local solver exit flag.
Elapsed time is 0.229525 seconds.

xgs =

    1.2592    0.4284


MultiStart completed the runs from all start points.

All 15 local solver runs converged with a positive local solver exit flag.
Elapsed time is 0.109984 seconds.

xms =

    1.2592    0.4284

Исследование результатов

Давайте исследуем результаты, которые возвратили оба решателя. Важная вещь отметить состоит в том, что результаты будут варьироваться на основе случайных стартовых точек, созданных для каждого решателя. Другой пробегает этот пример, может дать различные результаты. Координаты лучших результатов xgs и xms распечатанный к командной строке. Мы покажем уникальные результаты, возвращенные GlobalSearch и MultiStart и подсветите лучшие результаты каждого решателя, в терминах близости к глобальному решению.

Пятый выход каждого решателя является вектором, содержащим отличные минимумы (или максимумы, в этом случае) найденный. Мы построим (x, y) пары результатов, solsgs и solsms, против нашего контурного графика мы использовали прежде.

% Plot GlobalSearch results using the '*' marker
xGS = cell2mat({solsgs(:).X}');
scatter(xGS(:,1),xGS(:,2),'*','MarkerEdgeColor',[0 0 1],'LineWidth',1.25)

% Plot MultiStart results using a circle marker
xMS = cell2mat({solsms(:).X}');
scatter(xMS(:,1),xMS(:,2),'o','MarkerEdgeColor',[0 0 0],'LineWidth',1.25)
legend('Intensity','Bound','GlobalSearch','MultiStart','Location','best')

title('GlobalSearch and MultiStart Results')

Ослабление границ

С трудными границами на проблеме, обоих GlobalSearch и MultiStart смогли определить местоположение глобального максимума в этом запуске.

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

% Relax the bounds to spread out the start points
problem.lb = -5*ones(2,1);
problem.ub = 5*ones(2,1);

% Run GlobalSearch
tic;
[xgs,~,~,~,solsgs] = run(gs,problem);
toc
xgs

% Run MultiStart with 15 randomly generated points
tic;
[xms,~,~,~,solsms] = run(ms,problem,15);
toc
xms

GlobalSearch stopped because it analyzed all the trial points.

All 4 local solver runs converged with a positive local solver exit flag.
Elapsed time is 0.173760 seconds.

xgs =

    0.6571   -0.2096


MultiStart completed the runs from all start points.

All 15 local solver runs converged with a positive local solver exit flag.
Elapsed time is 0.134150 seconds.

xms =

    2.4947   -0.1439


% Show the contours
figure
contour(X,Y,Z)
hold on

% Create textarrow showing the location of xglobal
annotation('textarrow',[0.44 0.50],[0.63 0.58],'TextEdgeColor',[0 0 0],...
    'TextBackgroundColor',[1 1 1],'FontSize',12,'String',{'Global Max'});
axis([-1 3.75 -3 3])

% Plot GlobalSearch results using the '*' marker
xGS = cell2mat({solsgs(:).X}');
scatter(xGS(:,1),xGS(:,2),'*','MarkerEdgeColor',[0 0 1],'LineWidth',1.25)

% Plot MultiStart results using a circle marker
xMS = cell2mat({solsms(:).X}');
scatter(xMS(:,1),xMS(:,2),'o','MarkerEdgeColor',[0 0 0],'LineWidth',1.25)

% Highlight the best results from each:
% GlobalSearch result in red, MultiStart result in blue
plot(xgs(1),xgs(2),'sb','MarkerSize',12,'MarkerFaceColor',[1 0 0])
plot(xms(1),xms(2),'sb','MarkerSize',12,'MarkerFaceColor',[0 0 1])
legend('Intensity','GlobalSearch','MultiStart','Best GS','Best MS','Location','best')

title('GlobalSearch and MultiStart Results with Relaxed Bounds')

Лучший результат GlobalSearch показан красным квадратом и лучшим результатом MultiStart показан синим квадратом.

Настройка GlobalSearch Параметры

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

Заметьте, что GlobalSearch только запустил fmincon несколько раз. Чтобы увеличить шанс нахождения глобального максимума, мы хотели бы запустить больше точек. Чтобы ограничить набор стартовой точки кандидатами, скорее всего, чтобы найти глобальный максимум, мы дадим каждому решателю команду игнорировать стартовые точки, которые не удовлетворяют ограничениям путем установки StartPointsToRun свойство к bounds-ineqs. Кроме того, мы установим MaxWaitCycle и BasinRadiusFactor свойства так, чтобы GlobalSearch сможет идентифицировать узкий peaks быстро. Сокращение MaxWaitCycle причины GlobalSearch уменьшить радиус области притяжения BasinRadiusFactor чаще, чем с настройкой по умолчанию.

% Increase the total candidate points, but filter out the infeasible ones
gs = GlobalSearch(gs,'StartPointsToRun','bounds-ineqs', ...
    'MaxWaitCycle',3,'BasinRadiusFactor',0.3);
% Run GlobalSearch
tic;
xgs = run(gs,problem);
toc
xgs

GlobalSearch stopped because it analyzed all the trial points.

All 10 local solver runs converged with a positive local solver exit flag.
Elapsed time is 0.242955 seconds.

xgs =

    1.2592    0.4284

Использование MultiStartпараллельные возможности

Грубая сила способ улучшить наши возможности нахождения глобального максимума должна просто попробовать больше стартовых точек. Снова, это не может быть практично во всех ситуациях. В нашем случае мы только попробовали маленький набор до сих пор, и время выполнения не было ужасно долго. Так, разумно попробовать больше стартовых точек. Чтобы ускорить расчет, мы запустим MultiStart параллельно, если Parallel Computing Toolbox™ доступен.

% Set the UseParallel property of MultiStart
ms = MultiStart(ms,'UseParallel',true);

try
    demoOpenedPool = false;
    % Create a parallel pool if one does not already exist
    % (requires Parallel Computing Toolbox)
    if max(size(gcp)) == 0 % if no pool
        parpool
        demoOpenedPool = true;
    end
catch ME
    warning(message('globaloptim:globaloptimdemos:opticalInterferenceDemo:noPCT'));
end

% Run the solver
tic;
xms = run(ms,problem,100);
toc
xms

if demoOpenedPool
    % Make sure to delete the pool if one was created in this example
    delete(gcp) % delete the pool
end

MultiStart completed the runs from all start points.

All 100 local solver runs converged with a positive local solver exit flag.
Elapsed time is 0.956671 seconds.

xms =

    1.2592    0.4284

Объективные и нелинейные ограничения

Здесь мы перечисляем функции, которые задают задачу оптимизации:

function p = waveIntensity(x,amp,phase,relFreqConst,numSources,origins)
% WaveIntensity Intensity function for opticalInterferenceDemo.

%   Copyright 2009 The MathWorks, Inc.  

d = distanceFromSource(x,numSources,origins);
ampVec = [sum(amp./d .* cos(phase - d*relFreqConst));
    sum(amp./d .* sin(phase - d*relFreqConst))];

% Intensity is ||AmpVec||^2
p = ampVec'*ampVec;


function [c,ceq] = apertureConstraint(x,xcoords,ycoords)
% apertureConstraint Aperture constraint function for opticalInterferenceDemo.

%   Copyright 2009 The MathWorks, Inc.  

ceq = []; 
c = (x(1) - xcoords).^2 + (x(2) - ycoords).^2 - 9;


function d = distanceFromSource(v,numSources,origins)
% distanceFromSource Distance function for opticalInterferenceDemo.

%   Copyright 2009 The MathWorks, Inc.  

d = zeros(numSources,1);
for k = 1:numSources
    d(k) = norm(origins(k,:) - [v 0]);
end

Смотрите также

|

Похожие темы

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

Поддержка

Для просмотра документации необходимо авторизоваться на сайте