基于问题的工作流程中的供应导数- MATLAB & SimulinkgydF4y2Ba

基于问题的工作流中的供应派生gydF4y2Ba

为什么包括衍生品吗?gydF4y2Ba

这个例子展示了当自动导数不适用时，或者当您想包含Hessian时，如何在基于非线性问题的优化中包含导数信息，而Hessian在使用自动微分时是不提供的。在优化中包含梯度或Hessian可以带来以下好处:gydF4y2Ba

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

自动微分在优化中的应用gydF4y2Ba

从R2020b开始gydF4y2Ba解决gydF4y2Ba函数可以使用支持函数的自动微分，以便为求解器提供梯度。这些自动导数只适用于梯度(一阶导数)，而不适用于海森(二阶导数)。gydF4y2Ba

从R2022b开始，无论您是否使用，自动区分都适用gydF4y2Bafcn2optimexprgydF4y2Ba创建一个目标或约束函数。之前,使用gydF4y2Bafcn2optimexprgydF4y2Ba导致自动区分不适用。如果你需要使用gydF4y2Bafcn2optimexprgydF4y2Ba，这个例子展示了如何包含派生信息。gydF4y2Ba

在基于问题的优化中使用导数而不需要自动微分的方法是用gydF4y2Baprob2structgydF4y2Ba，然后编辑生成的目标函数和约束函数。这个例子展示了一种混合方法，其中自动微分为部分目标函数提供导数。gydF4y2Ba

创建优化问题gydF4y2Ba

与控制变量gydF4y2BaxgydF4y2Ba而且gydF4y2BaygydF4y2Ba，使用目标函数gydF4y2Ba

$\begin{array}{l} fun1gydF4y2Ba ＝gydF4y2Ba One hundred.gydF4y2Ba {（gydF4y2Ba ygydF4y2Ba -gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {（gydF4y2Ba 1gydF4y2Ba -gydF4y2Ba xgydF4y2Ba ）gydF4y2Ba}^{2gydF4y2Ba} \\ fun2gydF4y2Ba ＝gydF4y2Ba besseljgydF4y2Ba （gydF4y2Ba 1gydF4y2Ba ，gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba \\ 目标= fun1gydF4y2Ba +gydF4y2Ba fun2gydF4y2Ba ．gydF4y2Ba \end{array}$

包括平方和的约束条件gydF4y2BaxgydF4y2Ba而且gydF4y2BaygydF4y2Ba不超过4。gydF4y2Ba

fun2gydF4y2Ba不基于支持函数的优化表达式;看到gydF4y2Ba优化变量和表达式的支持操作gydF4y2Ba．因此,包括gydF4y2Bafun2gydF4y2Ba在优化问题中，必须使用gydF4y2Bafcn2optimexprgydF4y2Ba．gydF4y2Ba

若要在受支持的函数上使用AD，请在不受支持的函数的情况下设置问题gydF4y2Bafun2gydF4y2Ba,包括gydF4y2Bafun2gydF4y2Ba以后。gydF4y2Ba

概率= optimproblem;x = optimvar (gydF4y2Ba“x”gydF4y2Ba）;y = optimvar (gydF4y2Ba“y”gydF4y2Ba）;Fun1 = 100*(y - x^2)^2 + (1 - x)^2;概率。目标= fun1;probe . constraints .cons = x^2 + y^2 <= 4;gydF4y2Ba

将问题转化为基于求解器的形式gydF4y2Ba

包括的导数gydF4y2Bafun2gydF4y2Ba，首先转换问题没有gydF4y2Bafun2gydF4y2Ba结构中使用gydF4y2Baprob2structgydF4y2Ba．gydF4y2Ba

问题= prob2struct(概率,gydF4y2Ba.．.gydF4y2Ba“ObjectiveFunctionName”gydF4y2Ba，gydF4y2Ba“generatedObjectiveBefore”gydF4y2Ba）;gydF4y2Ba

在转换期间,gydF4y2Baprob2structgydF4y2Ba创建表示目标和非线性约束函数的函数文件。默认情况下，这些函数都有名称gydF4y2Ba“generatedObjective.m”gydF4y2Ba而且gydF4y2Ba“generatedConstraints.m”gydF4y2Ba．目标函数文件没有gydF4y2Bafun2gydF4y2Ba是gydF4y2Ba“generatedObjectiveBefore.m”gydF4y2Ba．gydF4y2Ba

生成的目标函数和约束函数包括梯度。gydF4y2Ba

计算导数和跟踪变量gydF4y2Ba

计算相关的导数gydF4y2Bafun2gydF4y2Ba．如果您有符号数学工具箱™许可证，则可以使用gydF4y2Ba梯度gydF4y2Ba(符号数学工具箱)gydF4y2Ba或gydF4y2Ba雅可比矩阵gydF4y2Ba(符号数学工具箱)gydF4y2Ba函数来帮助计算导数。看到gydF4y2Ba使用符号数学工具箱计算梯度和黑森斯gydF4y2Ba．gydF4y2Ba

基于求解器的方法有一个控制变量。每个优化变量(gydF4y2BaxgydF4y2Ba而且gydF4y2BaygydF4y2Ba，在本例中)是控制变量的一部分。您可以在生成的目标函数文件中找到优化变量和单个控制变量之间的映射gydF4y2Ba“generatedObjectiveBefore.m”gydF4y2Ba．文件的开头包含这些代码行或类似的代码行。gydF4y2Ba

% %变量指标。gydF4y2Baxidx = 1;yidx = 2;gydF4y2Ba将基于求解器的变量映射到基于问题的变量。gydF4y2Bax =数据源(xidx);y =数据源(yidx);gydF4y2Ba

在这段代码中，控制变量具有名称gydF4y2Ba数据源gydF4y2Ba．gydF4y2Ba

或者，您可以使用gydF4y2BavarindexgydF4y2Ba．gydF4y2Ba

idx = varindex(概率);disp (idx.x)gydF4y2Ba

1gydF4y2Ba

disp (idx.y)gydF4y2Ba

2gydF4y2Ba

完整的目标函数包括gydF4y2Bafun2gydF4y2Ba：gydF4y2Ba

Fun2 = besselj(1,x²+ y²);gydF4y2Ba

用标准微积分，计算gydF4y2Bagradfun2gydF4y2Ba的梯度gydF4y2Bafun2gydF4y2Ba．gydF4y2Ba

$gradfun2gydF4y2Ba ＝gydF4y2Ba ［gydF4y2Ba \begin{matrix} 2gydF4y2Ba xgydF4y2Ba （gydF4y2Ba besseljgydF4y2Ba （gydF4y2Ba 0gydF4y2Ba ，gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba -gydF4y2Ba besseljgydF4y2Ba （gydF4y2Ba 1gydF4y2Ba ，gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba /gydF4y2Ba （gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba ）gydF4y2Ba \\ 2gydF4y2Ba ygydF4y2Ba （gydF4y2Ba besseljgydF4y2Ba （gydF4y2Ba 0gydF4y2Ba ，gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba -gydF4y2Ba besseljgydF4y2Ba （gydF4y2Ba 1gydF4y2Ba ，gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba /gydF4y2Ba （gydF4y2Ba {xgydF4y2Ba}^{2gydF4y2Ba} +gydF4y2Ba {ygydF4y2Ba}^{2gydF4y2Ba} ）gydF4y2Ba ）gydF4y2Ba \end{matrix} ］gydF4y2Ba ．gydF4y2Ba$

编辑目标和约束文件gydF4y2Ba

编辑gydF4y2Ba“generatedObjectiveBefore.m”gydF4y2Ba包括gydF4y2Bafun2gydF4y2Ba．gydF4y2Ba

计算目标函数。Arg1 = (y - x.^2);最长= 100;长度= __arg1。^ 2;Arg4 = (1 - x);Obj = ((arg2 .* arg3) + arg4.^2);gydF4y2BaSSQ = x²+ y²;fun2 = besselj (ssq);Obj = Obj + fun2;gydF4y2Ba

通过编辑将计算出的梯度包含到目标函数文件中gydF4y2Ba“generatedObjectiveBefore.m”gydF4y2Ba．如果你有一个软件版本不执行梯度计算，包括所有这些行。如果您的软件执行梯度计算，则只包括粗体行，用于计算的梯度gydF4y2Bafun2gydF4y2Ba．gydF4y2Ba

计算目标梯度。如果nargout > 1 arg5 = 1;Arg6 = 0 ([2,1]);arg6 (xidx:) = (- (arg5。* 2 * (arg4 (:)))) + ((-(( arg5。*最长(:))。* 2 * (__arg1(:))))。* 2 * (x (:)));arg6 (yidx:) = ((arg5。*最长(:))。* 2。* (__arg1 (:)));研究生= arg6 (:);gydF4y2Baarg7 = besselj(0、ssq);Arg8 = arg7 - fun2/ssq;gfun = (2 * x * arg8;…2 * y * arg8];Grad = Grad + gfun;gydF4y2Ba结束gydF4y2Ba

回想一下，非线性约束是gydF4y2BaX²+ y²<= 4gydF4y2Ba．这个约束函数的梯度为gydF4y2Ba2 * (x, y)gydF4y2Ba．如果您的软件计算了约束梯度并将其包含在生成的约束文件中，那么您就不需要再做任何事情了。如果你的软件没有计算约束梯度，那么通过编辑包含非线性约束的梯度gydF4y2Ba“generatedConstraints.m”gydF4y2Ba包括这些行。gydF4y2Ba

在这里插入梯度计算。如果你包含一个渐变，通过设置% specificconstraintgradient选项为true来通知求解器。如果nargout > 2gydF4y2BacineqGrad = 2 * (x, y);gydF4y2BaceqGrad = [];结束gydF4y2Ba

运行有或没有梯度的问题gydF4y2Ba

使用基于求解器和基于问题(无梯度)的方法来运行问题，以查看差异。要使用派生信息运行基于求解器的问题，请在问题结构中创建适当的选项。gydF4y2Ba

选择= optimoptions (gydF4y2Ba“fmincon”gydF4y2Ba，gydF4y2Ba“SpecifyObjectiveGradient”gydF4y2Ba,真的,gydF4y2Ba.．.gydF4y2Ba“SpecifyConstraintGradient”gydF4y2Ba,真正的);问题。选择=选项;gydF4y2Ba

非线性问题需要一个非空的gydF4y2Bax0gydF4y2Ba问题结构中的字段。gydF4y2Ba

x0 = (1, 1);问题。x0 = x0;gydF4y2Ba

调用gydF4y2BafmincongydF4y2Ba在问题结构上。gydF4y2Ba

[xsolver, fvalsolver exitflagsolver outputsolver] = fmincon(问题)gydF4y2Ba

找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少，优化完成，在最优性公差的值内，约束满足在约束公差的值内。<停止条件详细信息> xsolver = 1.2494 1.5617 fvalsolver = -0.0038 exitflagsolver = 1 outputsolver = struct with fields: iterations: 15 funcCount: 32 constr违规:0 stepsize: 1.5569e-06 algorithm: ' internal -point' firstorderopt: 2.2058e-08 cgiterations: 7 message: ' ' '输入本地最小值发现满足约束。“enter”优化完成，因为目标函数在“enter”可行方向上不减少，在最优性容限范围内，“enter”和约束条件满足在约束容限范围内。“enter”“enter”优化完成:相对一阶优化度量2.125635e-08，“enter”小于“options”。OptimalityTolerance = 1.000000e-06，相对最大约束“违反”，0.000000e+00小于options。ConstraintTolerance = 1.000000 e-06。[1×1 struct]gydF4y2Ba

把这个解与得到的解比较一下gydF4y2Ba解决gydF4y2Ba没有导数信息。gydF4y2Ba

init。x＝x0（1）;init。y＝x0（2）;F2 = @(x,y)besselj(1,x²+ y²);fun2 = fcn2optimexpr (f2, x, y); prob.Objective = prob.Objective + fun2; [xproblem,fvalproblem,exitflagproblem,outputproblem] = solve(prob,init,.．.gydF4y2Ba“ConstraintDerivative”gydF4y2Ba，gydF4y2Ba“有限差分”gydF4y2Ba，gydF4y2Ba.．.gydF4y2Ba“ObjectiveDerivative”gydF4y2Ba，gydF4y2Ba“有限差分”gydF4y2Ba）gydF4y2Ba

使用fmincon解决问题。找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少，优化完成，在最优性公差的值内，约束满足在约束公差的值内。<停止条件详细信息> xproblem = struct with fields: x: 1.2494 y: 1.5617 fvalproblem = -0.0038 exitflagproblem = OptimalSolution outputproblem = struct with fields: iterations: 15 funcCount: 64 constr违规:0 stepsize: 1.5571e-06 algorithm: 'interior-point' firstorderopt: 6.0139e-07 cgiterations: 7 message: ' ' ' ' '本地最小值发现，满足约束。“enter”优化完成，因为目标函数在“enter”可行方向上不减少，在最优性容限范围内，“enter”和约束条件满足在约束容限范围内。“enter”“enter”优化完成:相对一阶优化度量5.795368e-07，“enter”小于“options”。OptimalityTolerance = 1.000000e-06，相对最大约束“违反”，0.000000e+00小于options。ConstraintTolerance = 1.000000 e-06。[1×1 struct]objectivederivative: "finite-differences" constraintderivative: "closed-form" solver: 'fmincon'

两个解决方案报告相同的函数值以显示精度，并且都需要相同的迭代次数。然而，具有梯度的解只需要32个函数求值，而没有梯度的解则需要64个函数求值。gydF4y2Ba

包括黑森gydF4y2Ba

要包含黑森，你必须使用gydF4y2Baprob2structgydF4y2Ba，即使优化表达式支持所有函数。这个例子展示了如何使用HessiangydF4y2BafmincongydF4y2Ba内点gydF4y2Ba算法。的gydF4y2BafminuncgydF4y2Ba信赖域gydF4y2Ba算法和gydF4y2BafmincongydF4y2Batrust-region-reflectivegydF4y2Ba算法使用不同的语法;看到gydF4y2Bafminunc信任区域或fmincon信任区域反射算法的HessiangydF4y2Ba．gydF4y2Ba

中描述的gydF4y2BaHessian for fmincon内点算法gydF4y2Ba黑森量就是拉格朗日量中的黑森量。gydF4y2Ba

{\nablagydF4y2Ba}_{xgydF4y2Ba xgydF4y2Ba}^{2gydF4y2Ba} lgydF4y2Ba （gydF4y2Ba xgydF4y2Ba ，gydF4y2Ba λgydF4y2Ba ）gydF4y2Ba ＝gydF4y2Ba {\nablagydF4y2Ba}^{2gydF4y2Ba} fgydF4y2Ba （gydF4y2Ba xgydF4y2Ba ）gydF4y2Ba +gydF4y2Ba \sumgydF4y2Ba {λgydF4y2Ba}_{ggydF4y2Ba ，gydF4y2Ba 我gydF4y2Ba} {\nablagydF4y2Ba}^{2gydF4y2Ba} {ggydF4y2Ba}_{我gydF4y2Ba} （gydF4y2Ba xgydF4y2Ba ）gydF4y2Ba +gydF4y2Ba \sumgydF4y2Ba {λgydF4y2Ba}_{hgydF4y2Ba ，gydF4y2Ba 我gydF4y2Ba} {\nablagydF4y2Ba}^{2gydF4y2Ba} {hgydF4y2Ba}_{我gydF4y2Ba} （gydF4y2Ba xgydF4y2Ba ）gydF4y2Ba ．gydF4y2Ba

(1）gydF4y2Ba

通过创建一个计算Hessian的函数文件来包含这个Hessian，并创建一个gydF4y2BaHessianFcngydF4y2Ba选择gydF4y2BafmincongydF4y2Ba使用黑森。要在本例中创建Hessian，需要分别创建客观约束和非线性约束的二阶导数。gydF4y2Ba

这个例子中目标函数的二阶导数矩阵有点复杂。其目标函数清单gydF4y2Bahessianfun (x)gydF4y2Ba是由符号数学工具箱使用相同的方法创建的gydF4y2Ba使用符号数学工具箱计算梯度和黑森斯gydF4y2Ba．gydF4y2Ba

函数gydF4y2Ba高频= hessfun (in)gydF4y2Ba% HESSFUNgydF4y2Ba% hf = hessfun (in1)gydF4y2Ba此函数是由符号数学工具箱8.6版本生成的。gydF4y2Ba% 10 - 8月- 2020 10:50:44gydF4y2Bax = in (1:);y = in (2:);t2 = x ^ 2;t4 = y ^ 2;t6 = x。* 4.0 e + 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.*.．.gydF4y2Ba(t2 * -3.0 + t4 + t2。* t5。* 2.0 + t3, t4 *。* 4.0 + t2。^ 3。* 2.0)* 2.0 + 2.0,gydF4y2Ba.．.gydF4y2Bat22 t22,gydF4y2Ba.．.gydF4y2Bat19-t12。* * (t2-t4 t15。* 3.0 + t2 * t5。* 4.0 +gydF4y2Ba.．.gydF4y2Bat3, t4 *。* 2.0 + t4。^ 3 * 2.0)* 2.0 + 2.0 e + 2], (2, 2));gydF4y2Ba

相比之下，非线性不等式约束的Hessian方法比较简单;它是单位矩阵的两倍。gydF4y2Ba

hessianc = 2 *眼(2);gydF4y2Ba

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

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

创建选项来使用这个黑森。gydF4y2Ba

problem.options.HessianFcn = H;gydF4y2Ba

解决问题并显示迭代次数和函数求值次数。这个解和以前差不多。gydF4y2Ba

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

找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少，优化完成，在最优性公差的值内，约束满足在约束公差的值内。8 10gydF4y2Ba

这一次,gydF4y2BafmincongydF4y2Ba只需要8次迭代，而不是15次;只需要10次函数计算，而不是32次。总之，提供一个解析的Hessian计算可以提高求解过程的效率，但是开发一个计算Hessian的函数可能是困难的。gydF4y2Ba

另请参阅gydF4y2Ba

prob2structgydF4y2Ba|gydF4y2BavarindexgydF4y2Ba|gydF4y2Bafcn2optimexprgydF4y2Ba