Пример для Optimization Toolbox™

В этом примере показано, как использовать два нелинейных решателя оптимизации и как установить опции. Нелинейными решателями, которые мы используем в этом примере, является fminunc и fmincon.

Все принципы, обрисованные в общих чертах в этом примере, применяются к другим нелинейным решателям, таким как fgoalattain, fminimax, lsqnonlin, lsqcurvefit, и fsolve.

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

Неограниченный пример оптимизации

Считайте задачу нахождения минимумом функции:

xexp(-(x2+y2))+(x2+y2)/20.

Постройте функцию, чтобы понять то, где она минимизирована

f = @(x,y) x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;
fsurf(f,[-2,2],'ShowContours','on')

График показывает, что минимум около точки (-1/2,0).

Обычно вы задаете целевую функцию как файл MATLAB. На данный момент эта функция достаточно проста задать как анонимная функция:

fun = @(x) f(x(1),x(2));

Возьмите предположение в решении:

x0 = [-.5; 0];

Установите опции оптимизации не использовать fminuncкрупномасштабный алгоритм по умолчанию, поскольку тот алгоритм требует, чтобы градиент целевой функции был обеспечен:

options = optimoptions('fminunc','Algorithm','quasi-newton');

Просмотрите итерации, когда решатель вычисляет:

options.Display = 'iter';

Вызовите fminunc, неограниченный нелинейный минимизатор:

[x, fval, exitflag, output] = fminunc(fun,x0,options);
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           3          -0.3769                         0.339
     1           6        -0.379694              1          0.286  
     2           9        -0.405023              1         0.0284  
     3          12        -0.405233              1        0.00386  
     4          15        -0.405237              1       3.17e-05  
     5          18        -0.405237              1       3.35e-08  

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

Решатель нашел решение в:

uncx = x
uncx = 2×1

   -0.6691
    0.0000

Значение функции в решении:

uncf = fval
uncf = -0.4052

Мы будем использовать количество функциональных оценок как мера КПД в этом примере. Общее количество функциональных оценок:

output.funcCount
ans = 18

Неограниченная оптимизация дополнительными параметрами

Мы теперь передадим дополнительные параметры в качестве дополнительных аргументов к целевой функции. Мы показываем два различных способа сделать это - использование файла MATLAB или использование вложенной функции.

Рассмотрите целевую функцию от предыдущего раздела:

f(x,y)=xexp(-(x2+y2))+(x2+y2)/20.

Мы параметризуем функцию с (a, b, c) следующим образом:

f(x,y,a,b,c)=(x-a)exp(-((x-a)2+(y-b)2))+((x-a)2+(y-b)2)/c.

Эта функция является переключенной и масштабированной версией исходной целевой функции.

Метод 1: Функция файла MATLAB

Предположим, что у нас есть целевая функция файла MATLAB под названием bowlpeakfun заданный как:

type bowlpeakfun
function y = bowlpeakfun(x, a, b, c)
%BOWLPEAKFUN Objective function for parameter passing in TUTDEMO.

%   Copyright 2008 The MathWorks, Inc.

y = (x(1)-a).*exp(-((x(1)-a).^2+(x(2)-b).^2))+((x(1)-a).^2+(x(2)-b).^2)/c;

Задайте параметры:

a = 2;
b = 3;
c = 10;

Создайте указатель анонимной функции на файл MATLAB:

f = @(x)bowlpeakfun(x,a,b,c)
f = function_handle with value:
    @(x)bowlpeakfun(x,a,b,c)

Вызовите fminunc найти минимум:

x0 = [-.5; 0];
options = optimoptions('fminunc','Algorithm','quasi-newton');
[x, fval] = fminunc(f,x0,options)
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
x = 2×1

    1.3639
    3.0000

fval = -0.3840

Метод 2: вложенная функция

Рассмотрите следующую функцию, которая реализует цель как вложенную функцию

type nestedbowlpeak
function [x,fval] =  nestedbowlpeak(a,b,c,x0,options)
%NESTEDBOWLPEAK Nested function for parameter passing in TUTDEMO.

%   Copyright 2008 The MathWorks, Inc.

[x,fval] = fminunc(@nestedfun,x0,options);  
    function y = nestedfun(x)
      y = (x(1)-a).*exp(-((x(1)-a).^2+(x(2)-b).^2))+((x(1)-a).^2+(x(2)-b).^2)/c;    
    end
end

В этом методе параметры (a, b, c) отображаются к вложенной целевой функции под названием nestedfun. Внешняя функция, nestedbowlpeak, вызовы fminunc и передает целевую функцию, nestedfun.

Задайте параметры, исходное предположение и опции:

a = 2;
b = 3;
c = 10;
x0 = [-.5; 0];
options = optimoptions('fminunc','Algorithm','quasi-newton');

Запустите оптимизацию:

[x,fval] =  nestedbowlpeak(a,b,c,x0,options)
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.
x = 2×1

    1.3639
    3.0000

fval = -0.3840

Вы видите, что оба метода произвели идентичные ответы, таким образом используйте, какой бы ни один вы находите самыми удобными.

Ограниченный пример оптимизации: Неравенства

Рассмотрите вышеупомянутую проблему с ограничением:

минимизировать xexp(-(x2+y2))+(x2+y2)/20,

удовлетворяющий  xy/2+(x+2)2+(y-2)2/22.

Ограничительное множество является внутренней частью наклоненного эллипса. Посмотрите на контуры целевой функции, построенной вместе с наклоненным эллипсом

f = @(x,y) x.*exp(-x.^2-y.^2)+(x.^2+y.^2)/20;
g = @(x,y) x.*y/2+(x+2).^2+(y-2).^2/2-2;
fimplicit(g)
axis([-6 0 -1 7])
hold on
fcontour(f)
plot(-.9727,.4685,'ro');
legend('constraint','f contours','minimum');
hold off

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

x0 = [-2 1];

Установите опции оптимизации: используйте алгоритм внутренней точки и включите отображение результатов в каждой итерации:

options = optimoptions('fmincon','Algorithm','interior-point','Display','iter');

Решатели требуют, чтобы нелинейные ограничительные функции дали два выходных параметров: один для нелинейных неравенств, второго для нелинейных равенств. Таким образом, мы пишем ограничение с помощью deal функция, чтобы дать оба выходных параметров:

gfun = @(x) deal(g(x(1),x(2)),[]);

Вызовите нелинейный ограниченный решатель. Нет никаких линейных равенств или неравенств или границ, так передача [] для тех аргументов:

[x,fval,exitflag,output] = fmincon(fun,x0,[],[],[],[],[],[],gfun,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       3    2.365241e-01    0.000e+00    1.972e-01
    1       6    1.748504e-01    0.000e+00    1.734e-01    2.260e-01
    2      10   -1.570560e-01    0.000e+00    2.608e-01    9.347e-01
    3      14   -6.629160e-02    0.000e+00    1.241e-01    3.103e-01
    4      17   -1.584082e-01    0.000e+00    7.934e-02    1.826e-01
    5      20   -2.349124e-01    0.000e+00    1.912e-02    1.571e-01
    6      23   -2.255299e-01    0.000e+00    1.955e-02    1.993e-02
    7      26   -2.444225e-01    0.000e+00    4.293e-03    3.821e-02
    8      29   -2.446931e-01    0.000e+00    8.100e-04    4.035e-03
    9      32   -2.446933e-01    0.000e+00    1.999e-04    8.126e-04
   10      35   -2.448531e-01    0.000e+00    4.004e-05    3.289e-04
   11      38   -2.448927e-01    0.000e+00    4.036e-07    8.156e-05

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.

Решение этой проблемы было найдено в:

x
x = 1×2

   -0.9727    0.4686

Значение функции в решении:

fval
fval = -0.2449

Общее количество функциональных оценок было:

Fevals = output.funcCount
Fevals = 38

Ограничению неравенства удовлетворяют в решении.

[c, ceq] = gfun(x)
c = -2.4608e-06
ceq =

     []

С тех пор c (x) близко к 0, ограничение "активно", означая, что ограничение влияет на решение. Вспомните, что неограниченное решение было найдено в

uncx
uncx = 2×1

   -0.6691
    0.0000

и неограниченная целевая функция, как находили, была

uncf
uncf = -0.4052

Ограничение переместило решение и увеличило цель на

fval-uncf
ans = 0.1603

Ограниченный пример оптимизации: предоставленные пользователями градиенты

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

минимизировать xexp(-(x2+y2))+(x2+y2)/20,

удовлетворяющий  xy/2+(x+2)2+(y-2)2/22.

Обеспечить градиент f (x) к fmincon, мы пишем целевую функцию в форме файла MATLAB:

type onehump
function [f,gf] = onehump(x)
% ONEHUMP Helper function for Tutorial for the Optimization Toolbox demo

%   Copyright 2008-2009 The MathWorks, Inc.

r = x(1)^2 + x(2)^2;
s = exp(-r);
f = x(1)*s+r/20;

if nargout > 1
   gf = [(1-2*x(1)^2)*s+x(1)/10;
       -2*x(1)*x(2)*s+x(2)/10];
end

Ограничение и его градиент содержатся в файле MATLAB tiltellipse:

type tiltellipse
function [c,ceq,gc,gceq] = tiltellipse(x)
% TILTELLIPSE Helper function for Tutorial for the Optimization Toolbox demo

%   Copyright 2008-2009 The MathWorks, Inc.

c = x(1)*x(2)/2 + (x(1)+2)^2 + (x(2)-2)^2/2 - 2;
ceq = [];

if nargout > 2
   gc = [x(2)/2+2*(x(1)+2);
       x(1)/2+x(2)-2];
   gceq = [];
end

Выскажите предположение в решении:

x0 = [-2; 1];

Установите опции оптимизации: мы продолжаем использовать тот же алгоритм в целях сравнения.

options = optimoptions('fmincon','Algorithm','interior-point');

Мы также устанавливаем опции использовать информацию о градиенте в ограничительных функциях и цели. Примечание: эти опции MUST быть включенными или информация о градиенте будут проигнорированы.

options = optimoptions(options,...
    'SpecifyObjectiveGradient',true,...
    'SpecifyConstraintGradient',true);

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

options.Display = 'iter';

Вызовите решатель:

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    2.365241e-01    0.000e+00    1.972e-01
    1       2    1.748504e-01    0.000e+00    1.734e-01    2.260e-01
    2       4   -1.570560e-01    0.000e+00    2.608e-01    9.347e-01
    3       6   -6.629161e-02    0.000e+00    1.241e-01    3.103e-01
    4       7   -1.584082e-01    0.000e+00    7.934e-02    1.826e-01
    5       8   -2.349124e-01    0.000e+00    1.912e-02    1.571e-01
    6       9   -2.255299e-01    0.000e+00    1.955e-02    1.993e-02
    7      10   -2.444225e-01    0.000e+00    4.293e-03    3.821e-02
    8      11   -2.446931e-01    0.000e+00    8.100e-04    4.035e-03
    9      12   -2.446933e-01    0.000e+00    1.999e-04    8.126e-04
   10      13   -2.448531e-01    0.000e+00    4.004e-05    3.289e-04
   11      14   -2.448927e-01    0.000e+00    4.036e-07    8.156e-05

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.

fmincon предполагаемые градиенты хорошо в предыдущем примере, таким образом, итерации в текущем примере подобны.

Решение этой проблемы было найдено в:

xold = x
xold = 2×1

   -0.9727
    0.4686

Значение функции в решении:

minfval = fval
minfval = -0.2449

Общее количество функциональных оценок было:

Fgradevals = output.funcCount
Fgradevals = 14

Сравните это с количеством функциональных оценок без градиентов:

Fevals
Fevals = 38

Изменение допусков завершения по умолчанию

На этот раз мы решаем ту же ограниченную задачу

минимизировать xexp(-(x2+y2))+(x2+y2)/20,

удовлетворяющий  xy/2+(x+2)2+(y-2)2/22,

более точно путем переопределения критериев завершения по умолчанию (опции. StepTolerance и опции. OptimalityTolerance). Мы продолжаем использовать градиенты. Значения по умолчанию для fminconалгоритм внутренней точки является опциями. StepTolerance = 1e-10, опции. OptimalityTolerance = 1e-6.

Замените два критерия завершения по умолчанию: погрешности на x для завершения и fval.

options = optimoptions(options,...
    'StepTolerance',1e-15,...
    'OptimalityTolerance',1e-8);

Вызовите решатель:

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    2.365241e-01    0.000e+00    1.972e-01
    1       2    1.748504e-01    0.000e+00    1.734e-01    2.260e-01
    2       4   -1.570560e-01    0.000e+00    2.608e-01    9.347e-01
    3       6   -6.629161e-02    0.000e+00    1.241e-01    3.103e-01
    4       7   -1.584082e-01    0.000e+00    7.934e-02    1.826e-01
    5       8   -2.349124e-01    0.000e+00    1.912e-02    1.571e-01
    6       9   -2.255299e-01    0.000e+00    1.955e-02    1.993e-02
    7      10   -2.444225e-01    0.000e+00    4.293e-03    3.821e-02
    8      11   -2.446931e-01    0.000e+00    8.100e-04    4.035e-03
    9      12   -2.446933e-01    0.000e+00    1.999e-04    8.126e-04
   10      13   -2.448531e-01    0.000e+00    4.004e-05    3.289e-04
   11      14   -2.448927e-01    0.000e+00    4.036e-07    8.156e-05
   12      15   -2.448931e-01    0.000e+00    4.000e-09    8.230e-07

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.

Мы теперь принимаем решение видеть больше десятичных чисел в решении, для того, чтобы видеть более точно значение, которое имеют новые допуски.

format long

Оптимизатор нашел решение в:

x
x = 2×1

  -0.972742227363546
   0.468569289098342

Сравните это с предыдущим значением:

xold
xold = 2×1

  -0.972742694488360
   0.468569966693330

Изменение

x - xold
ans = 2×1
10-6 ×

   0.467124813385844
  -0.677594988729435

Значение функции в решении:

fval
fval = 
  -0.244893137879894

Решение, улучшенное

fval - minfval
ans = 
    -3.996450220755676e-07

(это отрицательно, поскольку новое решение меньше),

Общее количество функциональных оценок было:

output.funcCount
ans = 
    15

Сравните это с количеством функциональных оценок, когда задача будет решена с обеспеченными пользователями градиентами, но с допусками по умолчанию:

Fgradevals
Fgradevals = 
    14

Ограниченный пример оптимизации с предоставленным пользователями гессианом

Если вы даете не только градиент, но также и Гессиан, решатели еще более точны и эффективны.

fminconрешатель внутренней точки берет матрицу Гессиана в качестве отдельной функции (не часть целевой функции). Функция Гессиана H (x, lambda) должна оценить Гессиан функции Лагранжа; см. Руководство пользователя для определения этого термина.

Решатели вычисляют значения lambda.ineqnonlin и lambda.eqlin; ваша функция Гессиана говорит решатели, как использовать эти значения.

В этой проблеме мы имеем, но одно ограничение неравенства, таким образом, Гессиан:

type hessfordemo
function H = hessfordemo(x,lambda)
% HESSFORDEMO Helper function for Tutorial for the Optimization Toolbox demo

%   Copyright 2008-2009 The MathWorks, Inc.

s = exp(-(x(1)^2+x(2)^2));
H = [2*x(1)*(2*x(1)^2-3)*s+1/10, 2*x(2)*(2*x(1)^2-1)*s;
        2*x(2)*(2*x(1)^2-1)*s, 2*x(1)*(2*x(2)^2-1)*s+1/10];
hessc = [2,1/2;1/2,1];
H = H + lambda.ineqnonlin(1)*hessc;

Для того, чтобы использовать Гессиан, необходимо установить опции соответственно:

options = optimoptions('fmincon',...
    'Algorithm','interior-point',...
    'SpecifyConstraintGradient',true,...
    'SpecifyObjectiveGradient',true,...
    'HessianFcn',@hessfordemo);

Допуски были задержаны к значениям по умолчанию. На этот раз должно быть меньше функциональных количеств.

options.Display = 'iter';

Вызовите решатель:

[x,fval,exitflag,output] = fmincon(@onehump,x0,[],[],[],[],[],[], ...
                                   @tiltellipse,options);
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    2.365241e-01    0.000e+00    1.972e-01
    1       3    5.821325e-02    0.000e+00    1.443e-01    8.728e-01
    2       5   -1.218829e-01    0.000e+00    1.007e-01    4.927e-01
    3       6   -1.421167e-01    0.000e+00    8.486e-02    5.165e-02
    4       7   -2.261916e-01    0.000e+00    1.989e-02    1.667e-01
    5       8   -2.433609e-01    0.000e+00    1.537e-03    3.486e-02
    6       9   -2.446875e-01    0.000e+00    2.057e-04    2.727e-03
    7      10   -2.448911e-01    0.000e+00    2.068e-06    4.191e-04
    8      11   -2.448931e-01    0.000e+00    2.001e-08    4.218e-06

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.

Были меньше и различные итерации на этот раз.

Решение этой проблемы было найдено в:

x
x = 2×1

  -0.972742246093537
   0.468569316215571

Значение функции в решении:

fval
fval = 
  -0.244893121872758

Общее количество функциональных оценок было:

output.funcCount
ans = 
    11

Сравните это с номером с помощью только оценки градиента с теми же допусками по умолчанию:

Fgradevals
Fgradevals = 
    14

Похожие темы