Первый пример: неограниченная минимизация с помощью гессена

Задача минимизации заключается в следующем:

Эта функция является положительной, с уникальным минимальным значением, равным нулю, полученным при 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 Interior-Point

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

Ограничения удерживают оптимизацию от глобальной минимальной точки [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.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: '...'
       bestfeasible: [1x1 struct]

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

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)

Feasible point with lower objective function value found.


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: '...'
       bestfeasible: [1x1 struct]

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