Этот пример показывает, как использовать функции GlobalSearch и MultiStart.
Этот пример показывает, как 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