Используя символьную математику с решателями Optimization Toolbox™

Этот пример использование:

Этот пример показывает, как использовать функции Symbolic Math Toolbox™ jacobian и matlabFunction, чтобы предоставить аналитические производные решателям оптимизации. Решатели Optimization Toolbox™ обычно более точны и эффективны, когда вы предоставляете градиенты и Гессианы ограничительных функций и цели.

Существует несколько факторов в использовании символьных вычислений с функциями оптимизации:

Цель оптимизации и ограничительные функции должны быть заданы с точки зрения вектора, сказать x. Однако символьные переменные являются скаляром или с комплексным знаком, не с векторным знаком. Это требует, чтобы вы перевели между векторами и скалярами.
Градиенты оптимизации, и иногда Гессианы, как предполагается, вычисляются в теле ограничительных функций или цели. Это означает, что символьный градиент или Гессиан должны быть помещены в соответствующее место в цели или ограничительном файле функции или указателе на функцию.
Вычисление градиентов и Гессианов символически может быть длительным. Поэтому необходимо выполнить это вычисление только однажды и сгенерировать код, через matlabFunction, чтобы вызвать во время выполнения решателя.
Выполнение символьных выражений с функцией subs длительно. Намного более эффективно использовать matlabFunction.
matlabFunction генерирует код, который зависит от ориентации входных векторов. Поскольку fmincon вызывает целевую функцию с вектор-столбцами, необходимо стараться вызвать matlabFunction с вектор-столбцами символьных переменных.

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

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

$f (x_{1}, x_{2}) = журнал (1 + 3 {(x_{2} - (x_{1}^{3} - x_{1}))}^{2} + (x_{1} - 4 / 3)^{2}) .$

Эта функция положительна, с уникальным минимальным значением нуля, достигнутого в 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 = 
 $журнал ({(x_{1} - \frac{4}{3})}^{2} + 3 {(- {x_{1}}^{3} + x_{1} + x_{2})}^{2} + 1)$

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

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

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

gradf = 
 $\begin{array}{l} (\begin{matrix} - \frac{6 (3 {x_{1}}^{2} - 1) (- {x_{1}}^{3} + x_{1} + x_{2}) - 2 x_{1} + \frac{8}{3}}{σ_{1}} \\ \frac{- 6 {x_{1}}^{3} + 6 x_{1} + 6 x_{2}}{σ_{1}} \end{matrix}) \\ где \\ σ_{1} = {(x_{1} - \frac{4}{3})}^{2} + 3 {(- {x_{1}}^{3} + x_{1} + x_{2})}^{2} + 1 \end{array}$

hessf = jacobian(gradf,x)

hessf = 
 $\begin{array}{l} (\begin{matrix} \frac{6 {(3 {x_{1}}^{2} - 1)}^{2} - 36 x_{1} (- {x_{1}}^{3} + x_{1} + x_{2}) + 2}{σ_{2}} - \frac{{σ_{3}}^{2}}{{σ_{2}}^{2}} & σ_{1} \\ σ_{1} & \frac{6}{σ_{2}} - \frac{{(- 6 {x_{1}}^{3} + 6 x_{1} + 6 x_{2})}^{2}}{{σ_{2}}^{2}} \end{matrix}) \\ где \\ σ_{1} = \frac{(- 6 {x_{1}}^{3} + 6 x_{1} + 6 x_{2}) σ_{3}}{{σ_{2}}^{2}} - \frac{18 {x_{1}}^{2} - 6}{σ_{2}} \\ σ_{2} = {(x_{1} - \frac{4}{3})}^{2} + 3 {(- {x_{1}}^{3} + x_{1} + x_{2})}^{2} + 1 \\ σ_{3} = 6 (3 {x_{1}}^{2} - 1) (- {x_{1}}^{3} + x_{1} + x_{2}) - 2 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 Алгоритм Внутренней точки

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

$5 \sinh (x_{2} / 5) \geq x_{1}^{4}$

$5 \tanh (x_{1} / 5) \geq x_{2}^{2} - 1 .$

Ограничения держат оптимизацию отдельно от глобальной минимальной точки [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

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