基于问题的工作流中的供应衍生品

为何要包括衍生品？

此示例说明了当自动导数不适用或当您想要包含 Hessian（使用自动微分无法提供）时，如何在基于非线性问题的优化中包含导数信息。在优化中加入梯度或 Hessian 可以带来以下好处：

更为稳健的结果。有限差分步骤有时会达到目标或非线性约束函数未定义、非有限或复杂的点。
提高了准确性。解析梯度比有限差分估计更准确。
收敛速度更快。包含 Hessian 可以加快收敛速度，这意味着更少的迭代次数。
提高了性能。与有限差分估计相比，解析梯度的计算速度更快，特别是对于稀疏结构的问题。然而，对于复杂的表达式，解析梯度的计算速度可能会更慢。

自动微分应用于优化

从 R2020b 开始，solve 函数可以对支持的函数使用自动微分，以便为求解器提供梯度。这些自动导数仅适用于梯度（一阶导数），而不适用于 Hessians（二阶导数）。

从 R2022b 开始，无论您是否使用 fcn2optimexpr 创建目标或约束函数，自动微分均适用。以前，使用 fcn2optimexpr 会导致自动微分不适用。如果您需要使用 fcn2optimexpr，此示例显示如何包含衍生信息。

在没有自动微分的基于问题的优化中使用导数的方法是使用 prob2struct 转换问题，然后编辑生成的目标函数和约束函数。此例展示的是一种混合方法，其中自动微分为部分目标函数提供导数。

创建优化问题

使用控制变量 x 和 y，使用目标函数

$\begin{array}{l} fun1 = 100 {(y - x^{2})}^{2} + {(1 - x)}^{2} \\ fun2 = besselj (1, x^{2} + y^{2}) \\ objective = 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 计算梯度和 Hessian。要通过与有限差分近似进行比较来验证手动生成的梯度函数，请使用 checkGradients。

基于求解器的方法有一个控制变量。每个优化变量（本例中为 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 次。

包括 Hessian

要包含 Hessian，您必须使用 prob2struct，即使您的所有函数都支持优化表达式。此示例展示了如何将 Hessian 用于 fmincon interior-point 算法。fminunc trust-region 算法和 fmincon trust-region-reflective 算法使用不同的语法；请参阅适用于 fminunc 信赖域或 fmincon 信赖域反射算法的黑塞函数。

正如适用于 fmincon 内点算法的黑塞函数中所描述的，Hessian 是拉格朗日的 Hessian。

\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)

通过创建一个函数文件来计算 Hessian，并创建一个 HessianFcn 选项让 fmincon 使用 Hessian，从而包含这个 Hessian。在这种情况下，要创建 Hessian，请分别创建目标和非线性约束的二阶导数。

本例的目标函数的二阶导数矩阵有些复杂。其目标函数列表 hessianfun(x) 由 Symbolic Math Toolbox 使用与使用 Symbolic Math Toolbox 计算梯度和 Hessian 中描述的相同方法创建。

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]);

相比之下，非线性不等式约束的 Hessian 很简单；它是单位矩阵的两倍。

hessianc = 2*eye(2);

创建拉格朗日函数的 Hessian 作为函数句柄。

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

创建选项以使用此 Hessian。

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 次）。总而言之，提供解析的 Hessian 计算可以提高解过程的效率，但是开发一个计算 Hessian 的函数可能会很困难。