基于问题的工作流程中的供给导数- MATLAB和Simulink万博1manbetxgydF4y2Ba

基于问题的工作流程中的供给导数gydF4y2Ba

为什么包含衍生品?gydF4y2Ba

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

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

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

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

当您不使用时，将应用自动区分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不是基于所支持的函数进行优化表达万博1manbetx式;看到gydF4y2Ba万博1manbetx支持优化变量和表达式的操作gydF4y2Ba．因此，要包括gydF4y2Bafun2gydF4y2Ba在优化问题中，必须将其转换为使用的优化表达式gydF4y2Bafcn2optimexprgydF4y2Ba．gydF4y2Ba

要在受支持的函数上使用AD，请在没万博1manbetx有不受支持的函数的情况下设置问题gydF4y2Bafun2gydF4y2Ba，并包括gydF4y2Bafun2gydF4y2Ba以后。gydF4y2Ba

Prob =优化问题;X = optimvar(gydF4y2Ba“x”gydF4y2Ba）;Y = optimvar(gydF4y2Ba“y”gydF4y2Ba）;Fun1 = 100*(y - x²)²+ (1 - x)²;概率。目标= 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 = inputVariables(xidx);y = inputVariables(yidx);gydF4y2Ba

在此代码中，控制变量的名称为gydF4y2Ba数据源gydF4y2Ba．gydF4y2Ba

或者，您可以使用gydF4y2BavarindexgydF4y2Ba．gydF4y2Ba

Idx = varindex(prob);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);Arg2 = 100;Arg3 = arg1.^2;Arg4 = (1 - x);Obj = ((arg2 .* arg3) + arg4.^2);gydF4y2BaSSQ = x²+ y²;Fun2 = besselj(1,ssq);Obj = Obj + fun2;gydF4y2Ba

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

计算目标梯度。arg5 = 1;Arg6 = 0 ([2,1]);arg6 (xidx:) = (- (arg5。* 2 * (arg4 (:)))) + ((-(( arg5。*最长(:))。* 2 * (__arg1(:))))。* 2 * (x (:)));Arg6 (yidx，:) = (arg5.*arg2(:)).*2.*(arg1(:)));Grad = 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

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

运行问题有和没有梯度gydF4y2Ba

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

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

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

X0 = [1;1];问题。X0 = X0;gydF4y2Ba

调用gydF4y2BafmincongydF4y2Ba在问题结构上。gydF4y2Ba

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

找到满足约束条件的局部最小值。优化完成是因为目标函数在可行方向上不递减，在最优性容差值范围内，约束条件满足在约束容差值范围内。 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: ' internal -point' firstorderopt: 2.2058e-08 cgiterations: 7 message: '(发现满足约束的局部最小值)。(()优化完成，因为目标函数在()可行方向上不下降，直到最优性公差的值之内，()且约束条件满足到约束公差的值之内。(() <停止条件详细信息>())(()优化完成:相对一阶最优性度量2.125635e-08，()小于选项。OptimalityTolerance = 1.000000e-06，相对最大约束(违例)为0.000000e+00，小于选项。约束容忍= 1.000000e-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 constrviolation: 0 stepsize: 1.5571e-06 algorithm: ' internal -point' firstorderopt: 6.0139e-07 cgiterations: 7 message: '(发现满足约束的局部最小值)。(()优化完成，因为目标函数在()可行方向上不下降，直到最优性公差的值之内，()且约束条件满足到约束公差的值之内。(() <停止条件详细信息>())(()优化完成:相对一阶最优性度量，5.795368e-07，()小于选项。OptimalityTolerance = 1.000000e-06，相对最大约束(违例)为0.000000e+00，小于选项。约束容忍= 1.000000e-06。(() '最佳可行:[1×1 struct]objectivederivative: "finite-differences" constraintderivative: "closed-form" solver: 'fmincon'

这两个解决万博尤文图斯方案报告相同的函数值以显示精度，并且都需要相同的迭代次数。然而，有梯度的解决方案只需要32个函数计算，而没有梯度的解决方案需要64个。gydF4y2Ba

包括黑森gydF4y2Ba

要包括黑森，你必须使用gydF4y2Baprob2structgydF4y2Ba，即使所有函数都支持优化表达式。万博1manbetx这个例子展示了如何使用黑森gydF4y2BafmincongydF4y2Ba内点gydF4y2Ba算法。的gydF4y2BafminuncgydF4y2Ba信赖域gydF4y2Ba算法和gydF4y2BafmincongydF4y2Batrust-region-reflectivegydF4y2Ba算法使用不同的语法;看到gydF4y2BaHessian for fminunc信任区域或fmincon信任区域反射算法gydF4y2Ba．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，并创建一个HessiangydF4y2BaHessianFcngydF4y2Ba选择gydF4y2BafmincongydF4y2Ba使用黑森。为了在这种情况下创建黑森，分别创建目标约束和非线性约束的二阶导数。gydF4y2Ba

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

函数gydF4y2BaHf = hessfun(in1)gydF4y2Ba% HESSFUNgydF4y2Ba% hf = hessfun (in1)gydF4y2Ba这个函数由符号数学工具箱8.6版本生成。gydF4y2Ba% 10-Aug-2020 10:50:44gydF4y2BaX = 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.*.．.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

黑头= 2*眼(2);gydF4y2Ba

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

H = @(x,lam)(hessianfun(x) + hessianc*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

fcn2optimexprgydF4y2Ba|gydF4y2Baprob2structgydF4y2Ba|gydF4y2BavarindexgydF4y2Ba