Exponenta Banner
exponenta event banner

Увеличение монохроматических поляризованных легких образцов вмешательства Используя 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)$Поэтому давайте ограничим пространство жизнеспособных решений, добавив к проблеме ограничения.

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

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

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

$$ \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 сможет быстро идентифицировать узкие пики. Сокращение 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

См. также

|

Связанные темы

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

Поддержка