基于问题的工作流中的供应衍生品
为何要包括衍生品?
此示例说明了当自动导数不适用或当您想要包含 Hessian(使用自动微分无法提供)时,如何在基于非线性问题的优化中包含导数信息。在优化中加入梯度或 Hessian 可以带来以下好处:
更为稳健的结果。有限差分步骤有时会达到目标或非线性约束函数未定义、非有限或复杂的点。
提高了准确性。解析梯度比有限差分估计更准确。
收敛速度更快。包含 Hessian 可以加快收敛速度,这意味着更少的迭代次数。
提高了性能。与有限差分估计相比,解析梯度的计算速度更快,特别是对于稀疏结构的问题。然而,对于复杂的表达式,解析梯度的计算速度可能会更慢。
自动微分应用于优化
从 R2020b 开始,solve
函数可以对支持的函数使用自动微分,以便为求解器提供梯度。这些自动导数仅适用于梯度(一阶导数),而不适用于 Hessians(二阶导数)。
从 R2022b 开始,无论您是否使用 fcn2optimexpr
创建目标或约束函数,自动微分均适用。以前,使用 fcn2optimexpr
会导致自动微分不适用。如果您需要使用 fcn2optimexpr
,此示例显示如何包含衍生信息。
在没有自动微分的基于问题的优化中使用导数的方法是使用 prob2struct
转换问题,然后编辑生成的目标函数和约束函数。此例展示的是一种混合方法,其中自动微分为部分目标函数提供导数。
创建优化问题
使用控制变量 x
和 y
,使用目标函数
包括 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)
1
disp(idx.y)
2
完整的目标函数包括 fun2
:
fun2 = besselj(1,x^2 + y^2);
使用标准微积分,计算 gradfun2
,即 fun2
的梯度。
编辑目标和约束文件
编辑 '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。
(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 的函数可能会很困难。
另请参阅
prob2struct
| varindex
| fcn2optimexpr
| checkGradients