主要内容

使用符号数学与优化工具箱™求解器

此示例显示了如何使用符号数学工具箱™功能雅可比矩阵matlabFunction为优化求解人员提供分析导数。当你提供目标函数和约束函数的梯度和Hessians时,最优化工具箱™求解器通常更准确和高效。

基于问题的优化可以自动计算和使用梯度;看到优化工具箱中的自动微分.有关使用自动区分的基于问题的示例,请参见受限静电非线性优化,基于问题

在使用优化函数的符号计算时,有几个注意事项:

  1. 优化目标和约束函数应该用向量来定义,例如x.然而,符号变量是标量或复数值,而不是向量值。这需要你在向量和标量之间进行转换。

  2. 优化梯度,有时是Hessians,应该在目标或约束函数的主体内计算。这意味着必须将符号梯度或黑森州人放置在目标或约束函数文件或功能手柄中的适当位置。

  3. 象征性地计算梯度和Hessians是很耗时的。因此,您应该只执行一次这个计算,并通过matlabFunction,在求解器执行期间调用。

  4. 的符号表达式求值subs函数是耗时的。它使用起来更有效率matlabFunction

  5. matlabFunction生成依赖于输入向量方向的代码。自粉刺用列向量调用目标函数时,必须小心调用matlabFunction用符号变量的列向量。

第一个例子:与黑森州的无约束最小化

最小化的目标函数为:

f x 1 x 2 日志 1 + 3. x 2 - x 1 3. - x 1 2 + x 1 - 4 / 3. 2

这个函数是正的,在点处有一个唯一的最小值为零x1= 4/3,x2=(4/3)^ 3 - 4/3 = 1.0370 ......

我们把自变量写成x1x2因为在这种形式中,它们可以用作符号变量。作为向量的组件x他们会被写下来x (1)x (2).函数有一个曲折的山谷,如下图所示。

信谊x1x2真实的x = [x1; x2];%符号变量列向量日志(f = 1 + 3 * (x2 - (x1 ^ 3 - x1)) ^ 2 + (x1 - 4/3) ^ 2)
f =

日志 x 1 - 4 3. 2 + 3. - x 1 3. + x 1 + x 2 2 + 1 日志((x1 -信谊(4/3))^ 2 + 3 * (- x1 ^ 3 + x1 + x2) ^ 2 + 1)

fsurf (f (2 - 2),“ShowContours”“上”)查看(127,38)

图中包含一个坐标轴。坐标轴包含一个函数曲面类型的对象。

计算f的梯度和Hessian:

gradf =雅可比矩阵(f (x)。%列gradf
gradf =

- 6 3. x 1 2 - 1 - x 1 3. + x 1 + x 2 - 2 x 1 + 8 3. σ 1 - 6 x 1 3. + 6 x 1 + 6 x 2 σ 1 在哪里 σ 1 x 1 - 4 3. 2 + 3. - x 1 3. + x 1 + x 2 2 + 1 (- (6 * (3 * x1 ^ 2 - 1) * (- x1 ^ 3 + x1 + x2) - 2 * x1 +符号(8/3))/ ((x1 -信谊(4/3))^ 2 + 3 * (- x1 ^ 3 + x1 + x2) ^ 2 + 1);(- 6 * x1 ^ 3 + 6 * x1 + 6 * x2) / ((x1 -信谊(4/3))^ 2 + 3 * (- x1 ^ 3 + x1 + x2) ^ 2 + 1))

hessf =雅可比矩阵(gradf, x)
Hessf =

6 3. x 1 2 - 1 2 - 36. x 1 - x 1 3. + x 1 + x 2 + 2 σ 2 - σ 3. 2 σ 2 2 σ 1 σ 1 6 σ 2 - - 6 x 1 3. + 6 x 1 + 6 x 2 2 σ 2 2 在哪里 σ 1 - 6 x 1 3. + 6 x 1 + 6 x 2 σ 3. σ 2 2 - 18. x 1 2 - 6 σ 2 σ 2 x 1 - 4 3. 2 + 3. - x 1 3. + x 1 + x 2 2 + 1 σ 3. 6 3. x 1 2 - 1 - x 1 3. + x 1 + x 2 - 2 x 1 + 8 3. [(6 *(3 * x1 ^ 2 - 1)^ 2 - 36 * x1 *( - x1 ^ 3 + x1 + x2)+ 2)/((x1 - sym(4/3))^ 2 + 3 *(-x1^3 + x1 + x2)^2 + 1) - (6*(3*x1^2 - 1)*(- x1^3 + x1 + x2) - 2*x1 + sym(8/3))^2/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1)^2, ((- 6*x1^3 + 6*x1 + 6*x2)*(6*(3*x1^2 - 1)*(- x1^3 + x1 + x2) - 2*x1 + sym(8/3)))/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1)^2 - (18*x1^2 - 6)/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1); ((- 6*x1^3 + 6*x1 + 6*x2)*(6*(3*x1^2 - 1)*(- x1^3 + x1 + x2) - 2*x1 + sym(8/3)))/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1)^2 - (18*x1^2 - 6)/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1), 6/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1) - (- 6*x1^3 + 6*x1 + 6*x2)^2/((x1 - sym(4/3))^2 + 3*(- x1^3 + x1 + x2)^2 + 1)^2]

Fminunc.解算器希望传入一个向量x,然后SpecifyObjectiveGradient选项设置为真的赫索斯·福克选项设置为“目标”,期望一个包含三个输出的列表:(f (x), gradf (x) hessf (x))

matlabFunction从三个输入列表中确切地生成三个输出列表。此外,使用vars.选项,matlabFunction接受向量输入。

跳频= matlabFunction (f, gradf hessf,'vars', {x});

现在求解从点[-1,2]开始的最小化问题:

选择= optimoptions ('fminunc'......“SpecifyObjectiveGradient”,真的,......“HessianFcn”“目标”......“算法”“信赖域”......“显示”“最后一次”);[xfinal, fval exitflag、输出]= fminunc (fh、[1,2]选项)
地方最低可能。Fminunc停止,因为功能值相对于其初始值的最终变化小于功能公差的值。
Xfinal =2×11.3333 - 1.0370
fval = 7.6623 e-12
exitflag = 3
输出=结构体字段:迭代:14 Funccount:15步骤:0.0027 Cgiterations:11:3.4391E-05算法:'信任区域'消息:'...'constrviolation:[]

将其与不使用梯度或Hessian信息的迭代次数进行比较。这需要“拟牛顿”算法。

选择= optimoptions ('fminunc'“显示”“最后一次”“算法”“拟牛顿”);fh2中= matlabFunction (f,'vars', {x});% fh2 =客观,无梯度或黑森[Xfinal,Fval,ExitFlag,Output2] = Fminunc(FH2,[ -  1; 2],选项)
发现本地最低限度。优化完成,因为梯度的大小小于最优耐受性的值。
Xfinal =2×11.3333 - 1.0371
fval = 2.1985e-11
exitflag = 1
Output2 =结构体字段:firstderopt: 2.4587e-06 algorithm: '拟牛顿' message: '…'

使用渐变和Hessians时,迭代次数较低,并且函数评估较少:

Sprintf([“使用渐变有%d次迭代”......“和黑森,但没有他们。”],......output.Iltations,Output2.Irtations)
ans = '使用梯度和Hessian有14次迭代,但没有它们有18次。'
Sprintf([“使用梯度进行了%d函数评估”......“和黑森,但没有他们。”],......output.funccount,output2.funccount)
ans = '使用梯度和Hessian进行了15次函数评估,但没有使用梯度和Hessian的有81次。'

第二个例子:fmincon内点算法的约束最小化

我们考虑相同的目标函数和起点,但现在有两个非线性约束:

5 sinh x 2 / 5 x 1 4

5 双曲正切 x 1 / 5 x 2 2 - 1

约束使优化远离全局最小值点[1.333,1.037]。想象这两个约束条件:

(X, Y) = meshgrid (2: .01:3);Z = (5*sinh(y /5) >= X.^4);%z = 1满足第一个约束,z = 0否则Z = Z+ 2*(5*tanh(X./5)) >= y ^2 - 1;% Z=2,其中第二个满足,Z=3,其中两个都满足冲浪(X, Y, Z,“线型”“没有”);图= GCF;图.Color =.' w '%白色背景视图(0,90)plot3(。43.96, .0373, 4,“o”“MarkerEdgeColor”“r”'Markersize'8);%最佳点包含(“x”)ylabel(“y”) 抓住

图中包含一个坐标轴。轴包含面、线两个对象。

我们在最优点周围画了一个红色的小圆。

这是目标函数在可行区域上的图,这个区域满足两个约束条件,上图用暗红色表示,还有一个围绕最优点的红色小圆:

W = log(1 + 3) *(Y - (X ^3 - X))^2 + (x - 4/3) ^2;%w =目标函数w(z <3)= nan;%只在满足约束的地方绘图冲浪(X, Y, W,“线型”“没有”);视图(68年,20)plot3(。43.96, .0373, .8152,“o”“MarkerEdgeColor”“r”......'Markersize'8);%最佳点包含(“x”)ylabel(“y”) zlabel (“z”) 抓住

图中包含一个坐标轴。轴包含面、线两个对象。

非线性约束必须以表单写入c (x) < = 0.我们计算所有符号约束及其衍生产品,并将它们放在函数句柄中使用matlabFunction

约束的梯度应该是列向量;它们必须以矩阵的形式置于目标函数中,矩阵的每一列表示一个约束函数的梯度。这是这个形式的转置雅可比矩阵我们取下面的转置。

我们将非线性约束放入函数句柄中。粉刺期望非线性约束和梯度按顺序输出[c ceq gradc gradceq].由于没有非线性等式约束,我们输出[]量表信gradceq

C1 = x1^4 - 5*sinh(x2/5);C2 = x2^2 - 5*tanh(x1/5) - 1;C = [c1 c2];gradc =雅可比矩阵(c、x) ';把%转置成正确的形式约束= matlabFunction (gradc, c, [] [],'vars', {x});

内部点算法要求其Hessian函数被写入单独的函数,而不是作为目标函数的一部分。这是因为非线性约束的函数需要在其黑森州中包含这些限制。它的黑森州是拉格朗日的幽灵;有关更多信息,请参阅用户指南。

Hessian函数需要两个输入参数:位置矢量x和拉格朗日乘法器结构lambda。用于非线性约束的Lambda结构的部分是lambda.ineqnonlin.lambda.eqnonlin.对于电流约束,没有线性等式,所以我们使用两个乘数lambda.ineqnonlin (1)lambda.ineqnonlin (2)

在第一个例子中,我们计算了目标函数的Hessian。现在我们计算这两个约束函数的Hessians,并使函数句柄版本为matlabFunction

hessc1 =雅可比矩阵(gradc (: 1), x);%约束=第一个c列hessc2 =雅可比矩阵(gradc (:, 2), x);hessfh = matlabFunction (hessf,'vars', {x});hessc1h = matlabFunction (hessc1,'vars', {x});hessc2h = matlabFunction (hessc2,'vars', {x});

为了得到最后的Hessian,我们把三个Hessian放在一起,在约束函数中加入适当的拉格朗日乘数。

myhess = @(x,lambda)(hessfh(x)+......lambda.ineqnonlin (1) * hessc1h (x) +......lambda.ineqnonlin (2) * hessc2h (x));

设置使用内点算法、梯度和Hessian的选项,使目标函数同时返回目标和梯度,并运行求解器:

选择= optimoptions (“fmincon”......“算法”“内点”......“SpecifyObjectiveGradient”,真的,......'specifyconstraintgradient',真的,......“HessianFcn”myhess,......“显示”“最后一次”);% fh2 =不含黑森纤维的客观物fh2中= matlabFunction (f gradf'vars', {x});[xfinal, fval exitflag、输出]= fmincon (fh2中,[1,2],......[],[],[],[],[],[], 约束,选项)
找到满足约束条件的局部最小值。优化完成,因为目标函数在可行的方向上不降低,到在最优性公差的值内,并且对约束公差的值满足约束。
Xfinal =2×10.4396 - 0.0373
fval = 0.8152
exitflag = 1
输出=结构体字段:迭代:10 Funccount:13 CorroudViration:0步骤:1.9160E-06算法:'内部点'Firstordopt:1.9217E-08 Cgiterations:0消息:'...'最佳:[1x1 struct]

再次,求解器使许多更少的迭代和函数计算梯度和Hessian提供比当他们不是:

选择= optimoptions (“fmincon”“算法”“内点”......“显示”“最后一次”);% fh3 =客观无梯度或黑森fh3 = matlabfunction(f,'vars', {x});%无梯度约束:约束= matlabFunction (c []'vars', {x});[xfinal, fval exitflag output2] = fmincon (fh3 [1, 2],......[],[],[],[],[],[], 约束,选项)
找到了较低目标函数值的可行点。找到满足约束条件的局部最小值。优化完成,因为目标函数在可行的方向上不降低,到在最优性公差的值内,并且对约束公差的值满足约束。
Xfinal =2×10.4396 - 0.0373
fval = 0.8152
exitflag = 1
Output2 =结构体字段:stepsize: 9.6632e-07 algorithm: ' internal -point' firstderopt: 3.8435e-07 cgiterations: 0 message: '…'最佳可行:[1x1 struct]
Sprintf([“使用渐变有%d次迭代”......“和黑森,但没有他们。”],......output.Iltations,Output2.Irtations)
ans ='使用渐变和黑森州有10个迭代,但是没有它们。“
Sprintf([“使用梯度进行了%d函数评估”......“和黑森,但没有他们。”],......output.funccount,output2.funccount)
ans ='使用渐变和黑森州有13个功能评估,但没有它们。“

清理符号变量

本例中使用的符号变量被假定为实数。要从符号引擎工作区中清除这个假设,仅删除变量是不够的。您必须使用语法清除变量的假设

假设((x1, x2),“清楚”

当以下命令的输出为空时,所有假设都被清除

假设([x1,x2])
ans =空符号:1 × 0

相关话题