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

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

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

  1. Цель оптимизации и ограничительные функции должны быть заданы в терминах вектора, сказать x. Однако символьные переменные являются скаляром или с комплексным знаком, не с векторным знаком. Это требует, чтобы вы перевели между векторами и скалярами.

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

  3. Вычисление градиентов и Гессианов символически может быть длительным. Поэтому необходимо выполнить это вычисление только однажды и сгенерировать код через matlabFunction, вызывать во время выполнения решателя.

  4. Выполнение символьных выражений с subs функция длительна. Намного более эффективно использовать matlabFunction.

  5. matlabFunction генерирует код, который зависит от ориентации входных векторов. Начиная с fmincon вызывает целевую функцию с вектор-столбцами, необходимо стараться вызвать matlabFunction с вектор-столбцами символьных переменных.

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

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

f(x1,x2)=log(1+3(x2-(x13-x1))2+(x1-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 = 

log(x1-432+3-x13+x1+x22+1)журнал ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1)

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

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

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

(-63x12-1-x13+x1+x2-2x1+83σ1-6x13+6x1+6x2σ1)where  σ1=x1-432+3-x13+x1+x22+1[-(6* (3*x1^2 - 1) * (-x1^3 + x1 + x2) - 2*x1 + sym (8/3)) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1); (-6*x1^3 + 6*x1 + 6*x2) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1)]

hessf = jacobian(gradf,x)
hessf = 

(63x12-12-36x1-x13+x1+x2+2σ2-σ32σ22σ1σ16σ2--6x13+6x1+6x22σ22)where  σ1=-6x13+6x1+6x2σ3σ22-18x12-6σ2  σ2=x1-432+3-x13+x1+x22+1  σ3=63x12-1-x13+x1+x2-2x1+83[(6* (3*x1^2 - 1) ^2 - 36*x1* (-x1^3 + x1 + x2) + 2) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1) - (6* (3*x1^2 - 1) * (-x1^3 + x1 + x2) - 2*x1 + sym (8/3)) ^2 / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1) ^2, ((-6*x1^3 + 6*x1 + 6*x2) * (6* (3*x1^2 - 1) * (-x1^3 + x1 + x2) - 2*x1 + sym (8/3))) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1) ^2 - (18*x1^2 - 6) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1); ((-6*x1^3 + 6*x1 + 6*x2) * (6* (3*x1^2 - 1) * (-x1^3 + x1 + x2) - 2*x1 + sym (8/3))) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1) ^2 - (18*x1^2 - 6) / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1), 6 / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1) - (-6*x1^3 + 6*x1 + 6*x2) ^2 / ((x1 - sym (4/3)) ^2 + 3* (-x1^3 + x1 + x2) ^2 + 1) ^2]

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 Алгоритм Внутренней точки

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

5sinh(x2/5)x14

5tanh(x1/5)x22-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
 

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте