Сравните суррогатную оптимизацию с другими решателями

Этот пример сравнивает surrogateopt к двум другим решателям: fmincon, рекомендуемый решатель для сглаженных проблем и patternsearch, рекомендуемый решатель для несглаженных проблем. Пример использует несглаженную функцию на двумерной области.

type nonSmoothFcn

function [f, g] = nonSmoothFcn(x)
%NONSMOOTHFCN is a non-smooth objective function

%   Copyright 2005 The MathWorks, Inc.

for i = 1:size(x,1)
    if  x(i,1) < -7
        f(i) = (x(i,1))^2 + (x(i,2))^2 ;
    elseif x(i,1) < -3
        f(i) = -2*sin(x(i,1)) - (x(i,1)*x(i,2)^2)/10 + 15 ;
    elseif x(i,1) < 0
        f(i) = 0.5*x(i,1)^2 + 20 + abs(x(i,2))+ patho(x(i,:));
    elseif x(i,1) >= 0
        f(i) = .3*sqrt(x(i,1)) + 25 +abs(x(i,2)) + patho(x(i,:));
    end
end

%Calculate gradient
g = NaN;
if x(i,1) < -7
    g = 2*[x(i,1); x(i,2)];
elseif x(i,1) < -3
    g = [-2*cos(x(i,1))-(x(i,2)^2)/10; -x(i,1)*x(i,2)/5];
elseif x(i,1) < 0
    [fp,gp] = patho(x(i,:));
    if x(i,2) > 0
        g = [x(i,1)+gp(1); 1+gp(2)];
    elseif x(i,2) < 0
        g =  [x(i,1)+gp(1); -1+gp(2)];
    end
elseif x(i,1) >0
    [fp,gp] = patho(x(i,:));
    if x(i,2) > 0
        g = [.15/sqrt(x(i,1))+gp(1); 1+ gp(2)];
    elseif x(i,2) < 0
        g = [.15/sqrt(x(i,1))+gp(1); -1+ gp(2)];
    end
end

function [f,g] = patho(x)
Max = 500;
f = zeros(size(x,1),1);
g = zeros(size(x));
for k = 1:Max  %k 
   arg = sin(pi*k^2*x)/(pi*k^2);
   f = f + sum(arg,2);
   g = g + cos(pi*k^2*x);
end

mplier = 0.1; % Scale the control variable
Objfcn = @(x)nonSmoothFcn(mplier*x); % Handle to the objective function
range = [-6 6;-6 6]/mplier; % Range used to plot the objective function
rng default % Reset the global random number generator
showNonSmoothFcn(Objfcn,range);
title('Nonsmooth Objective Function')
view(-151,44)

drawnow

Смотрите как хорошо surrogateopt делает в определении местоположения глобального минимума в количестве по умолчанию итераций.

lb = -6*ones(1,2)/mplier;
ub = -lb;
[xs,fvals,eflags,outputs] = surrogateopt(Objfcn,lb,ub);

Surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

fprintf("Lowest found value = %g.\r",fvals)

Lowest found value = 13.

figure
showNonSmoothFcn(Objfcn,range);
view(-151,44)
hold on
p1 = plot3(xs(1),xs(2),fvals,'om','MarkerSize',15,'MarkerFaceColor','m');
legend(p1,{'Solution'})
hold off

Сравните с `patternsearch`

Установите patternsearch опции, чтобы использовать то же количество функциональных оценок, начинающих со случайной точки в границах.

rng default
x0 = lb + rand(size(lb)).*(ub - lb);
optsps = optimoptions('patternsearch','MaxFunctionEvaluations',200,'PlotFcn','psplotbestf');
[xps,fvalps,eflagps,outputps] = patternsearch(Objfcn,x0,[],[],[],[],lb,ub,[],optsps);

Optimization terminated: mesh size less than options.MeshTolerance.

figure
showNonSmoothFcn(Objfcn,range);
view(-151,44)
hold on
p1 = plot3(x0(1),x0(2),Objfcn(x0),'ob','MarkerSize',12,'MarkerFaceColor','b');
p2 = plot3(xps(1),xps(2),fvalps,'om','MarkerSize',15,'MarkerFaceColor','m');
legend([p1,p2],{'Start Point','Solution'})
hold off

patternsearch найденный тем же решением как surrogateopt.

Ограничьте количество функциональных оценок и попробуйте еще раз.

optsurr = optimoptions('surrogateopt','MaxFunctionEvaluations',40);
[xs,fvals,eflags,outputs] = surrogateopt(Objfcn,lb,ub,optsurr);

Surrogateopt stopped because it exceeded the function evaluation limit set by 
'options.MaxFunctionEvaluations'.

optsps.MaxFunctionEvaluations = 40;
[xps,fvalps,eflagps,outputps] = patternsearch(Objfcn,x0,[],[],[],[],lb,ub,[],optsps);

Maximum number of function evaluations exceeded: increase options.MaxFunctionEvaluations.

Снова, оба решателя нашли глобальное решение быстро.

Сравните с `fmincon`

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

Установите fmincon опции, чтобы использовать функцию построения графика, то же количество функциональных оценок как предыдущие решатели и та же стартовая точка как patternsearch.

opts = optimoptions('fmincon','PlotFcn','optimplotfval','MaxFunctionEvaluations',200);
[fmsol,fmfval,eflag,fmoutput] = fmincon(Objfcn,x0,[],[],[],[],lb,ub,[],opts);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

figure
showNonSmoothFcn(Objfcn,range);
view(-151,44)
hold on
p1 = plot3(x0(1),x0(2),Objfcn(x0),'ob','MarkerSize',12,'MarkerFaceColor','b');
p2 = plot3(fmsol(1),fmsol(2),fmfval,'om','MarkerSize',15,'MarkerFaceColor','m');
legend([p1,p2],{'Start Point','Solution'})
hold off

fmincon застревает в локальном минимуме около стартовой точки.

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