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

Этот пример использует:

Открыть Live Script

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

Основанная на проблеме оптимизация может вычислять и использовать градиенты автоматически; см. Раздел «Автоматическая дифференциация» в Optimization Toolbox (Optimization Toolbox). Для основанного на проблеме примера, использующего автоматическую дифференциацию, см. «Ограниченная электростатическая нелинейная оптимизация», «Основанная на проблеме» (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)

Figure contains an axes. The axes contains an object of type functionsurface.

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

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

(-6 3 x12-1 -x13+x1+x2-2 x1+83σ1-6 x13+6 x1+6 x2σ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 = 

(6 3 x12-12-36 x1 -x13+x1+x2+2σ2-σ32σ22σ1σ16σ2--6 x13+6 x1+6 x22σ22)where  σ1=-6 x13+6 x1+6 x2 σ3σ22-18 x12-6σ2  σ2=x1-432+3 -x13+x1+x22+1  σ3=6 3 x12-1 -x13+x1+x2-2 x1+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]

The 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

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

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

Figure contains an axes. The axes contains 2 objects of type surface, line.

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

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

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

Figure contains an axes. The axes contains 2 objects of type surface, line.

Нелинейные ограничения должны быть записаны в форме 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

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

Поддержка

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