Первый пример: безусловная минимизация с гессианом

Целевая функция, чтобы минимизировать:

Эта функция положительна с уникальным минимальным значением нуля, достигнутого в x1 = 4/3, x2 =(4/3)^3 - 4/3 = 1.0370...

Мы пишем независимые переменные как x1 и x2 потому что в этой форме они могут использоваться в качестве символьных переменных. Как компоненты векторного x они были бы записанным x(1) и x(2). Функция имеет извилистый овраг, как изображено в графике ниже.

syms x1 x2 real
x = [x1;x2]; % column vector of symbolic variables
f = log(1 + 3*(x2 - (x1^3 - x1))^2 + (x1 - 4/3)^2)

f = 
$$\log \left({\left(x_1-\frac{4}{3}\right)}^2+3\hspace{0.17em}{\left(-{x_1}^3+x_1+x_2\right)}^2+1\right)$$

fsurf(f,[-2 2],'ShowContours','on')
view(127,38)

Вычислите градиент и Гессиан f:

gradf = jacobian(f,x).' % column gradf

gradf = 
$$\begin{array}{l}\left(\begin{array}{c}-\frac{6\hspace{0.17em}\left(3\hspace{0.17em}{x_1}^2-1\right)\hspace{0.17em}\left(-{x_1}^3+x_1+x_2\right)-2\hspace{0.17em}{x}_1+\frac{8}{3}}{{\sigma }_1}\\ \frac{-6\hspace{0.17em}{x_1}^3+6\hspace{0.17em}{x}_1+6\hspace{0.17em}{x}_2}{{\sigma }_1}\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1={\left(x_1-\frac{4}{3}\right)}^2+3\hspace{0.17em}{\left(-{x_1}^3+x_1+x_2\right)}^2+1\end{array}$$

hessf = jacobian(gradf,x)

hessf = 
$$\begin{array}{l}\left(\begin{array}{cc}\frac{6\hspace{0.17em}{\left(3\hspace{0.17em}{x_1}^2-1\right)}^2-36\hspace{0.17em}{x}_1\hspace{0.17em}\left(-{x_1}^3+x_1+x_2\right)+2}{{\sigma }_2}-\frac{{{\sigma }_3}^2}{{{\sigma }_2}^2}& {\sigma }_1\\ {\sigma }_1& \frac{6}{{\sigma }_2}-\frac{{\left(-6\hspace{0.17em}{x_1}^3+6\hspace{0.17em}{x}_1+6\hspace{0.17em}{x}_2\right)}^2}{{{\sigma }_2}^2}\end{array}\right)\\ \\ \mathrm{where}\\ \\ \mathrm{ }{\sigma }_1=\frac{\left(-6\hspace{0.17em}{x_1}^3+6\hspace{0.17em}{x}_1+6\hspace{0.17em}{x}_2\right)\hspace{0.17em}{\sigma }_3}{{{\sigma }_2}^2}-\frac{18\hspace{0.17em}{x_1}^2-6}{{\sigma }_2}\\ \\ \mathrm{ }{\sigma }_2={\left(x_1-\frac{4}{3}\right)}^2+3\hspace{0.17em}{\left(-{x_1}^3+x_1+x_2\right)}^2+1\\ \\ \mathrm{ }{\sigma }_3=6\hspace{0.17em}\left(3\hspace{0.17em}{x_1}^2-1\right)\hspace{0.17em}\left(-{x_1}^3+x_1+x_2\right)-2\hspace{0.17em}{x}_1+\frac{8}{3}\end{array}$$

fminunc решатель ожидает передавать в векторе x, и, с SpecifyObjectiveGradient набор опции к true и HessianFcn набор опции к 'objective', ожидает список трех выходных параметров: [f(x),gradf(x),hessf(x)].

matlabFunction генерирует точно этот список трех выходных параметров из списка трех входных параметров. Кроме того, использование vars опция, matlabFunction принимает векторные входные параметры.

fh = matlabFunction(f,gradf,hessf,'vars',{x});

Теперь решите задачу минимизации, запускающуюся в точке [-1,2]:

options = optimoptions('fminunc', ...
    'SpecifyObjectiveGradient', true, ...
    'HessianFcn', 'objective', ...
    'Algorithm','trust-region', ...
    'Display','final');
[xfinal,fval,exitflag,output] = fminunc(fh,[-1;2],options)

Local minimum possible.

fminunc stopped because the final change in function value relative to 
its initial value is less than the value of the function tolerance.

xfinal = 2×1

    1.3333
    1.0370

fval = 7.6623e-12

exitflag = 3

output = struct with fields:
         iterations: 14
          funcCount: 15
           stepsize: 0.0027
       cgiterations: 11
      firstorderopt: 3.4391e-05
          algorithm: 'trust-region'
            message: '...'
    constrviolation: []

Сравните это с количеством итераций, не используя градиента или информации о Гессиане. Это требует 'quasi-newton' алгоритм.

options = optimoptions('fminunc','Display','final','Algorithm','quasi-newton');
fh2 = matlabFunction(f,'vars',{x}); 
% fh2 = objective with no gradient or Hessian
[xfinal,fval,exitflag,output2] = fminunc(fh2,[-1;2],options)

Local minimum found.

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

xfinal = 2×1

    1.3333
    1.0371

fval = 2.1985e-11

exitflag = 1

output2 = struct with fields:
       iterations: 18
        funcCount: 81
         stepsize: 2.1164e-04
     lssteplength: 1
    firstorderopt: 2.4587e-06
        algorithm: 'quasi-newton'
          message: '...'

Количество итераций ниже при использовании градиентов и Гессианов, и существует существенно меньше вычислений функции:

sprintf(['There were %d iterations using gradient' ...
    ' and Hessian, but %d without them.'], ...
    output.iterations,output2.iterations)

ans = 
'There were 14 iterations using gradient and Hessian, but 18 without them.'

sprintf(['There were %d function evaluations using gradient' ...
    ' and Hessian, but %d without them.'], ...
    output.funcCount,output2.funcCount)

ans = 
'There were 15 function evaluations using gradient and Hessian, but 81 without them.'

Второй Пример: Ограниченная Минимизация Используя fmincon Алгоритм Внутренней точки

Мы рассматриваем ту же целевую функцию и начальную точку, но теперь имеем два нелинейных ограничения:

Ограничения держат оптимизацию отдельно от глобальной минимальной точки [1.333,1.037]. Визуализируйте эти два ограничения:

[X,Y] = meshgrid(-2:.01:3);
Z = (5*sinh(Y./5) >= X.^4); 
% Z=1 where the first constraint is satisfied, Z=0 otherwise
Z = Z+ 2*(5*tanh(X./5) >= Y.^2 - 1); 
% Z=2 where the second is satisfied, Z=3 where both are
surf(X,Y,Z,'LineStyle','none');
fig = gcf;
fig.Color = 'w'; % white background
view(0,90)
hold on
plot3(.4396, .0373, 4,'o','MarkerEdgeColor','r','MarkerSize',8); 
% best point
xlabel('x')
ylabel('y')
hold off

Мы построили маленький красный круг вокруг оптимальной точки.

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

W = log(1 + 3*(Y - (X.^3 - X)).^2 + (X - 4/3).^2); 
% W = the objective function
W(Z < 3) = nan; % plot only where the constraints are satisfied
surf(X,Y,W,'LineStyle','none');
view(68,20)
hold on
plot3(.4396, .0373, .8152,'o','MarkerEdgeColor','r', ...
    'MarkerSize',8); % best point
xlabel('x')
ylabel('y')
zlabel('z')
hold off

Нелинейные ограничения должны быть написаны в форме c(x) <= 0. Мы вычисляем все символьные ограничения и их производные, и размещаем их в указатель на функцию с помощью matlabFunction.

Градиенты ограничений должны быть вектор-столбцами; они должны быть размещены в целевую функцию как матрица с каждым столбцом матрицы, представляющей градиент одной ограничительной функции. Это - транспонирование формы, сгенерированной jacobian, таким образом, мы берем транспонирование ниже.

Мы помещаем нелинейные ограничения в указатель на функцию. fmincon ожидает, что нелинейные ограничения и градиенты будут выведены в порядке [c ceq gradc gradceq]. С тех пор нет никаких нелинейных ограничений равенства, мы выход [] для ceq и gradceq.

c1 = x1^4 - 5*sinh(x2/5);
c2 = x2^2 - 5*tanh(x1/5) - 1;
c = [c1 c2];
gradc = jacobian(c,x).'; % transpose to put in correct form
constraint = matlabFunction(c,[],gradc,[],'vars',{x});

Алгоритм внутренней точки требует, чтобы его функция Гессиана была написана как отдельная функция, вместо того, чтобы быть частью целевой функции. Это вызвано тем, что нелинейно ограниченная функция должна включать те ограничения в свой Гессиан. Его Гессианом является Гессиан функции Лагранжа; см. Руководство пользователя для получения дополнительной информации.

Функция Гессиана берет два входных параметра: радиус-вектор x, и lambda структуры множителя Лагранжа. Частями структуры lambda, которую вы используете в нелинейных ограничениях, является lambda.ineqnonlin и lambda.eqnonlin. Для текущего ограничения нет никаких линейных равенств, таким образом, мы используем эти два множителя lambda.ineqnonlin(1) и lambda.ineqnonlin(2).

Мы вычислили Гессиан целевой функции в первом примере. Теперь мы вычисляем Гессианы двух ограничительных функций и делаем версии указателя на функцию с matlabFunction.

hessc1 = jacobian(gradc(:,1),x); % constraint = first c column
hessc2 = jacobian(gradc(:,2),x);

hessfh = matlabFunction(hessf,'vars',{x});
hessc1h = matlabFunction(hessc1,'vars',{x});
hessc2h = matlabFunction(hessc2,'vars',{x});

Чтобы сделать итоговый Гессиан, мы соединяем эти три Гессиана, добавляя соответствующие множители Лагранжа в ограничительные функции.

myhess = @(x,lambda)(hessfh(x) + ...
    lambda.ineqnonlin(1)*hessc1h(x) + ...
    lambda.ineqnonlin(2)*hessc2h(x));

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

options = optimoptions('fmincon', ...
    'Algorithm','interior-point', ...
    'SpecifyObjectiveGradient',true, ...
    'SpecifyConstraintGradient',true, ...
    'HessianFcn',myhess, ...
    'Display','final');
% fh2 = objective without Hessian
fh2 = matlabFunction(f,gradf,'vars',{x});
[xfinal,fval,exitflag,output] = fmincon(fh2,[-1;2],...
    [],[],[],[],[],[],constraint,options)

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.

xfinal = 2×1

    0.4396
    0.0373

fval = 0.8152

exitflag = 1

output = struct with fields:
         iterations: 10
          funcCount: 13
    constrviolation: 0
           stepsize: 1.9160e-06
          algorithm: 'interior-point'
      firstorderopt: 1.9217e-08
       cgiterations: 0
            message: '...'

Снова, решатель делает много меньшим количеством итераций и вычислений функции с градиентом и Гессианом предоставленный чем тогда, когда они не:

options = optimoptions('fmincon','Algorithm','interior-point',...
    'Display','final');
% fh3 = objective without gradient or Hessian
fh3 = matlabFunction(f,'vars',{x});
% constraint without gradient:
constraint = matlabFunction(c,[],'vars',{x});
[xfinal,fval,exitflag,output2] = fmincon(fh3,[-1;2],...
    [],[],[],[],[],[],constraint,options)

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.

xfinal = 2×1

    0.4396
    0.0373

fval = 0.8152

exitflag = 1

output2 = struct with fields:
         iterations: 18
          funcCount: 57
    constrviolation: 0
           stepsize: 9.6632e-07
          algorithm: 'interior-point'
      firstorderopt: 3.8435e-07
       cgiterations: 0
            message: '...'

sprintf(['There were %d iterations using gradient' ...
    ' and Hessian, but %d without them.'],...
    output.iterations,output2.iterations)

ans = 
'There were 10 iterations using gradient and Hessian, but 18 without them.'

sprintf(['There were %d function evaluations using gradient' ...
    ' and Hessian, but %d without them.'], ...
    output.funcCount,output2.funcCount)

ans = 
'There were 13 function evaluations using gradient and Hessian, but 57 without them.'

Чистка символьных переменных

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

assume([x1,x2],'clear')

Все предположения очищены, когда выход следующей команды пуст

assumptions([x1,x2])

 
ans =
 
Empty sym: 1-by-0