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

Введение

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

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

где фаза в нуле времени для источника, является скоростью света, является частотой света, является амплитудой источника и является расстоянием от источника до.

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

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

% 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')

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

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

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

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

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

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

подвергающийся

где и координаты и апертурный радиус точечного источника, соответственно. Каждому источнику дают апертуру с радиусом 3. Данные границы охватывают выполнимую область.

Цель () и нелинейные ограничительные функции задана в отдельных файлах 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 default value of the optimality tolerance,
and constraints are satisfied to within the default 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.624779 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.341500 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.482697 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.432552 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.622113 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 2.641997 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