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

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

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

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

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

Начиная с 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}) \\ цель = фун1 + 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) функция, которая помогает вычислить производные. См. «Вычисление градиентов и гессианов с использованием символьных математических 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 trust-region или fmincon trust- области reflective алгоритмов.

Как описано в Гессиане для алгоритма внутренней точки 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 с использованием того же подхода, как описано в Calculate Gradients and Hessians Using 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. Сводные данные, предоставление аналитического вычисления Гессиана может улучшить эффективность процесса решения, но разработка функции для вычисления Гессиана может оказаться трудной.

См. также

fcn2optimexpr | prob2struct | varindex

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