Максимизация монохроматических поляризованных шаблонов интерференции света с помощью 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.

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

% 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. Данные границы охватывают допустимую область.

Функции object $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

Поддержка

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