Предоставьте производные в основанном на проблеме рабочем процессе

Почему включают производные?

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

Больше устойчивых результатов. Конечные шаги дифференцирования иногда достигают точек, где цель или нелинейная ограничительная функция являются неопределенными, не конечными, или комплексными.
Увеличенная точность. Аналитические градиенты могут быть более точными, чем оценки конечной разности.
Более быстрая сходимость. Включая Гессиан может привести к более быстрой сходимости, означая меньше итераций.
Улучшенная производительность. Аналитические градиенты могут быть быстрее, чтобы вычислить, чем оценки конечной разности, особенно для проблем с разреженной структурой. Для сложных выражений, однако, аналитические градиенты могут быть медленнее, чтобы вычислить.

Автоматическое дифференцирование, примененное оптимизация

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

Автоматическое дифференцирование применяется, когда вы не используете fcn2optimexpr создать ограничительная функция или цель. Если необходимо использовать fcn2optimexpr, этот пример показывает, как включать производную информацию.

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

Создайте задачу оптимизации

С контрольными переменными x и y, используйте целевую функцию

$\begin{array}{l} fun1 = 100 {(y - x^{2})}^{2} + {(1 - x)}^{2} \\ fun2 = besselj (1, x^{2} + y^{2}) \\ цель = fun1 + fun2 . \end{array}$

Включайте ограничение что сумма квадратов x и y не больше, чем 4.

fun2 не основан на поддерживаемых функциях для выражений оптимизации; смотрите Поддерживаемые Операции для Переменных и выражений Оптимизации. Поэтому включать fun2 в задаче оптимизации необходимо преобразовать его в использование выражения оптимизации fcn2optimexpr.

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

prob = optimproblem;
x = optimvar('x');
y = optimvar('y');
fun1 = 100*(y - x^2)^2 + (1 - x)^2;
prob.Objective = fun1;
prob.Constraints.cons = x^2 + y^2 <= 4;

Преобразуйте проблему в основанную на решателе форму

Включать производные fun2, сначала преобразуйте проблему без fun2 к использованию структуры prob2struct.

problem = prob2struct(prob,...
    'ObjectiveFunctionName','generatedObjectiveBefore');

Во время преобразования, prob2struct создает файлы функции, которые представляют объективные и нелинейные ограничительные функции. По умолчанию эти функции имеют имена 'generatedObjective.m' и 'generatedConstraints.m'. Файл целевой функции без fun2 'generatedObjectiveBefore.m'.

Сгенерированная цель и ограничительные функции включают градиенты.

Вычислите производные и отслеживайте переменные

Вычислите производные, сопоставленные с fun2. Если у вас есть лицензия Symbolic Math Toolbox™, можно использовать gradient (Symbolic Math Toolbox) или jacobian (Symbolic Math Toolbox) функция, чтобы помочь вычислить производные. Смотрите Вычисляют Градиенты и Гессианы Используя Symbolic Math Toolbox™.

Основанный на решателе подход имеет одну контрольную переменную. Каждая переменная оптимизации (x и y, в этом примере), фрагмент контрольной переменной. Можно найти отображение между переменными оптимизации и одной контрольной переменной в сгенерированном файле целевой функции 'generatedObjectiveBefore.m'. Начало файла содержит эти строки кода или подобные линии.

%% Variable indices.
xidx = 1;
yidx = 2;

%% Map solver-based variables to problem-based.
x = inputVariables(xidx);
y = inputVariables(yidx);

В этом коде контрольная переменная имеет имя inputVariables.

В качестве альтернативы можно найти отображение при помощи varindex.

idx = varindex(prob);
disp(idx.x)

disp(idx.y)

Полная целевая функция включает fun2:

fun2 = besselj(1,x^2 + y^2);

Используя стандартное исчисление, вычислите gradfun2, градиент fun2.

$gradfun2 = [\begin{matrix} 2 x (besselj (0, x^{2} + y^{2}) - besselj (1, x^{2} + y^{2}) / (x^{2} + y^{2})) \\ 2 y (besselj (0, x^{2} + y^{2}) - besselj (1, x^{2} + y^{2}) / (x^{2} + y^{2})) \end{matrix}] .$

Отредактируйте файлы цели и ограничения

Отредактируйте 'generatedObjectiveBefore.m' включать fun2.

%% Compute objective function.
arg1 = (y - x.^2);
arg2 = 100;
arg3 = arg1.^2;
arg4 = (1 - x);
obj = ((arg2 .* arg3) + arg4.^2);

ssq = x^2 + y^2;
fun2 = besselj(1,ssq);
obj = obj + fun2;

Включайте расчетные градиенты в файл целевой функции путем редактирования 'generatedObjectiveBefore.m'. Если у вас есть версия программного обеспечения, которая не выполняет вычисление градиента, включает все эти линии. Если ваше программное обеспечение выполняет вычисление градиента, включайте только полужирные линии, которые вычисляют градиент fun2.

%% Compute objective gradient.
if nargout > 1
    arg5 = 1;
    arg6 = zeros([2, 1]);
    arg6(xidx,:) = (-(arg5.*2.*(arg4(:)))) + ((-((arg5.*arg2(:)).*2.*(arg1(:)))).*2.*(x(:)));
    arg6(yidx,:) = ((arg5.*arg2(:)).*2.*(arg1(:)));
    grad = arg6(:);
    
    arg7 = besselj(0,ssq);
    arg8 = arg7 - fun2/ssq;
    gfun = [2*x*arg8;...
        2*y*arg8];
    
    grad = grad + gfun;
end

Вспомните, что нелинейным ограничением является x^2 + y^2 <= 4. Градиентом этой ограничительной функции является 2*[x;y]. Если ваше программное обеспечение вычисляет ограничительный градиент и включает его в сгенерированный ограничительный файл, то вы не должны делать ничего больше. Если ваше программное обеспечение не вычисляет ограничительный градиент, то включайте градиент нелинейного ограничения путем редактирования 'generatedConstraints.m' включать эти линии.

%% Insert gradient calculation here.
% If you include a gradient, notify the solver by setting the
% SpecifyConstraintGradient option to true.
if nargout > 2
    cineqGrad = 2*[x;y];
    ceqGrad = [];
end

Запустите проблему с и без градиентов

Запустите проблему с помощью и основанного на решателе и основанного на проблеме (никакой градиент) методы, чтобы видеть различия. Чтобы запустить основанную на решателе проблему с помощью производной информации, создайте подходящие варианты в структуре задачи.

options = optimoptions('fmincon','SpecifyObjectiveGradient',true,...
    'SpecifyConstraintGradient',true);
problem.options = options;

Нелинейные проблемы требуют непустого x0 поле в структуре задачи.

x0 = [1;1];
problem.x0 = x0;

Вызвать fmincon на структуре задачи.

[xsolver,fvalsolver,exitflagsolver,outputsolver] = fmincon(problem)

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.

<stopping criteria details>

xsolver =

    1.2494
    1.5617


fvalsolver =

   -0.0038


exitflagsolver =

     1


outputsolver = 

  struct with fields:

         iterations: 15
          funcCount: 32
    constrviolation: 0
           stepsize: 1.5569e-06
          algorithm: 'interior-point'
      firstorderopt: 2.2058e-08
       cgiterations: 7
            message: '↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 2.125635e-08,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
       bestfeasible: [1×1 struct]

Сравните это решение с тем, полученным из solve без производной информации.

init.x = x0(1);
init.y = x0(2);
f2 = @(x,y)besselj(1,x^2 + y^2);
fun2 = fcn2optimexpr(f2,x,y);
prob.Objective = prob.Objective + fun2;
[xproblem,fvalproblem,exitflagproblem,outputproblem] = solve(prob,init,...
    "ConstraintDerivative","finite-differences",...
    "ObjectiveDerivative","finite-differences")

Solving problem using fmincon.

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.

<stopping criteria details>

xproblem = 

  struct with fields:

    x: 1.2494
    y: 1.5617


fvalproblem =

   -0.0038


exitflagproblem = 

    OptimalSolution


outputproblem = 

  struct with fields:

              iterations: 15
               funcCount: 64
         constrviolation: 0
                stepsize: 1.5571e-06
               algorithm: 'interior-point'
           firstorderopt: 6.0139e-07
            cgiterations: 7
                 message: '↵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.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 5.795368e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
            bestfeasible: [1×1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

Оба решения сообщают о том же значении функции точности отображения, и оба требуют того же количества итераций. Однако решение с градиентами требует только 32 вычислений функции, по сравнению с 64 для решения без градиентов.

Включайте гессиан

Чтобы включать Гессиан, необходимо использовать prob2struct, даже если все ваши функции поддерживаются для выражений оптимизации. В этом примере показано, как использовать Гессиан для fmincon interior-point алгоритм. fminunc trust-region алгоритм и fmincon trust-region-reflective использование алгоритма различный синтаксис; смотрите Гессиан для fminunc доверительной области или области доверия fmincon отражающие алгоритмы.

Как описано в Гессиане для fmincon алгоритма внутренней точки Гессианом является Гессиан функции Лагранжа.

\nabla_{x x}^{2} L (x, λ) = \nabla^{2} f (x) + \sum λ_{g, i} \nabla^{2} g_{i} (x) + \sum λ_{h, i} \nabla^{2} h_{i} (x) .

(1)

Включайте этот Гессиан путем создания файла функции, чтобы вычислить Гессиан и создания HessianFcn опция для fmincon использовать Гессиан. Чтобы создать Гессиан в этом случае, создайте вторые производные объективных и нелинейных ограничений отдельно.

Вторая производная матрица целевой функции для этого примера является несколько сложной. Его целевая функция, перечисляющая hessianfun(x) был создан Symbolic Math Toolbox с помощью того же подхода как описано в, Вычисляют Градиенты и Гессианы Используя Symbolic Math Toolbox™.

function hf = hessfun(in1)
%HESSFUN
%    HF = HESSFUN(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.6.
%    10-Aug-2020 10:50:44

x = in1(1,:);
y = in1(2,:);
t2 = x.^2;
t4 = y.^2;
t6 = x.*4.0e+2;
t3 = t2.^2;
t5 = t4.^2;
t7 = -t4;
t8 = -t6;
t9 = t2+t4;
t10 = t2.*t4.*2.0;
t11 = besselj(0,t9);
t12 = besselj(1,t9);
t13 = t2+t7;
t14 = 1.0./t9;
t16 = t3+t5+t10-2.0;
t15 = t14.^2;
t17 = t11.*t14.*x.*y.*4.0;
t19 = t11.*t13.*t14.*2.0;
t18 = -t17;
t20 = t12.*t15.*t16.*x.*y.*4.0;
t21 = -t20;
t22 = t8+t18+t21;
hf = reshape([t2.*1.2e+3-t19-y.*4.0e+2-t12.*t15.*...
    (t2.*-3.0+t4+t2.*t5.*2.0+t3.*t4.*4.0+t2.^3.*2.0).*2.0+2.0,...
    t22,t22,...
    t19-t12.*t15.*(t2-t4.*3.0+t2.*t5.*4.0+...
    t3.*t4.*2.0+t4.^3.*2.0).*2.0+2.0e+2],[2,2]);

В отличие от этого Гессиан нелинейного ограничения неравенства прост; это - дважды единичная матрица.

hessianc = 2*eye(2);

Создайте Гессиан функции Лагранжа как указатель на функцию.

H = @(x,lam)(hessianfun(x) + hessianc*lam.ineqnonlin);

Создайте опции, чтобы использовать этот Гессиан.

problem.options.HessianFcn = H;

Решите задачу и отобразите количество итераций и количество вычислений функции. Решение является приблизительно тем же самым как прежде.

[xhess,fvalhess,exitflaghess,outputhess] = fmincon(problem);
disp(outputhess.iterations)
disp(outputhess.funcCount)

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.

    8

    10

На этот раз, fmincon берет только 8 итераций по сравнению с 15 и только 10 вычислений функции по сравнению с 32. Таким образом, обеспечение аналитического вычисления Гессиана может повысить эффективность процесса решения, но разработка функции, чтобы вычислить Гессиан может затруднить.

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