使用符号数学工具箱计算梯度和黑森斯
这个例子展示了如何获得更快和更健壮的解决非线性优化问题使用fmincon以及符号数学工具箱™函数。如果你有符号数学工具箱许可,你可以使用这些符号数学工具箱函数轻松计算目标函数和约束函数的分析梯度和Hessians:
雅可比矩阵
(符号数学工具箱)生成标量函数的梯度,并生成向量函数的偏导数矩阵。例如,你可以通过应用得到Hessian矩阵(目标函数的二阶导数)雅可比矩阵梯度。这个例子展示了如何使用雅可比矩阵生成目标函数和约束函数的符号梯度和Hessians。
matlabFunction
(符号数学工具箱)生成计算符号表达式值的匿名函数或文件。这个例子展示了如何使用matlabFunction生成计算目标函数和约束函数及其在任意点上的导数的文件。
与优化工具箱™函数相比,符号数学工具箱函数具有不同的语法和结构。特别是,符号变量是实数或复标量,而优化工具箱函数传递向量参数。因此,您需要采取几个步骤,以一种适合于的内点算法的形式符号化地生成目标函数、约束条件及其所有必要的导数fmincon.
基于问题的优化可以自动计算和使用梯度;看到优化工具箱中的自动区分.有关使用自动区分来解决此问题的基于问题的方法,请参见基于优化变量的约束静电非线性优化.
问题描述
考虑将10个电子放入导电体的静电问题。受限于处于体内,电子以一种使其总势能最小的方式排列自己。所有的电子都在物体的边界上。电子是不可区分的,所以这个问题不存在唯一的最小值(排列一个解中的电子会给出另一个有效解)。这个例子的灵感来自多兰、Moré和曼森[58].
这个例子是一个由不平等定义的指导机构
.
这个身体看起来像球体上的金字塔。
(X, Y) = meshgrid (1: .01:1);Z1 = -abs(X) -abs(Y);Z2 = -1 -√(1 - x ^2 - y ^2);Z2 =实际(Z2);W1 = Z1;W2 = Z2;W1(Z1 < Z2) = nan;只在Z1 > Z2处画点W2(Z1 < Z2) = nan;只在Z1 > Z2处画点手=图;%该图形的句柄,用于稍后的更多绘图集(gcf,“颜色”,' w ')%白色背景冲浪(X, Y, W1,“线型”,“没有”);持有在冲浪(X, Y W2,“线型”,“没有”);视图(-44年,18)
这幅画的上下表面之间有一点缝隙。这个间隙是用于创建图形的一般绘图例程的产物。这个程序擦除一个表面上与另一个表面接触的任何矩形补丁。
创建变量
生成一个符号向量x由实数符号变量组成的30 × 1向量xij,我1到10之间,和j1到3之间。这些变量代表了电子的三个坐标我:ξ1对应于
坐标,ξ2)对应于
坐标,xi3对应于
坐标。
X = cell(3,10);为i = 1:10为J = 1:3 x{J,i} = sprintf(“x % d % d ',我,j);结束结束x = x (:);现在x是一个30 × 1的向量x =符号(x,“真实”的);
显示向量x.
x
x =
包括线性约束
写出线性约束条件
每一个电子都有四个线性不等式:
Xi3 - xi1 - xi2≤0 Xi3 - xi1 + xi2≤0 Xi3 + xi1 - xi2≤0 Xi3 + xi1 + xi2≤0
这个问题总共有40个线性不等式。
用结构化的方式写出不等式。
B = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1);一个= 0(40岁,30);为i=1:10 A(4*i-3:4*i,3*i-2:3*i) = B;结束b = 0(40岁,1);
你们可以看到* x≤b代表着不平等。
disp (A * x)
建立非线性约束及其梯度和黑森
非线性约束也是结构化的。
.
生成约束及其梯度和hessian。
c =符号(0 (10));我= 1:10;c = (x(3 *我2)。^ 2 + x(3张)。^ 2 +(x(3*i)+1).^2 - 1).'; gradc = jacobian(c,x).';% .'执行转置Hessc = cell(1,10);为I = 1:10 hessc{I} =雅可比矩阵(gradc(:, I),x);结束
约束向量c是行向量,和的梯度c(我)表示为我矩阵的第Th列gradc.此表格是正确的,如非线性约束.
海赛矩阵,hessc {1}、……hessc {10},是正方形和对称的。此示例将它们存储在单元格数组中,这比将它们存储在诸如hessc1、……hessc10.
使用.”语法的转置。的”语法意味着共轭转置,它有不同的符号导数。
建立目标函数及其梯度和黑森函数
目标函数,势能,是每对电子之间距离的倒数之和。
.
距离是矢量分量差的平方和的平方根。
计算能量(目标函数)及其梯度和黑森函数。
能量=信谊(0);为i = 1:3:25为J = i+3:3:28 dist = x(i:i+2) - x(J: J +2);Energy = Energy + 1/√(dist.'*dist);结束结束gradenergy =雅可比矩阵(能源、x) ';hessenergy =雅可比矩阵(gradenergy x);
创建目标函数文件
目标函数有两个输出,能源而且gradenergy.调用时将两个函数放在一个向量中matlabFunction以减少子表达式的数量matlabFunction生成,并且仅当调用函数(fmincon在本例中)请求两个输出。您可以将生成的文件放在MATLAB路径上的任何文件夹中。在本例中,将它们放在当前文件夹中。
curdir = [pwd filesep];您可能需要使用currdir = pwd文件名= [currdir,“demoenergy.m”];matlabFunction(能源、gradenergy“文件”文件名,“var”, {x});
matlabFunction返回能源作为第一个输出和gradenergy第二。该函数也接受单个输入向量{x}而不是输入列表x11、……x103.
结果文件demoenergy.m包含以下行或类似行:
函数(能源、gradenergy) = demoenergy (in)% DEMOENERGY% [energy, gradenergy] = demoenergy (in1)...: x101 = in(28日);...能量= 1./t140.^(1./2) +...;如果nargout > 1...gradenergy = [(t174。* (t185 - 2。* x11))。/ 2 -...];结束
该函数具有具有梯度的目标函数的正确形式;看到编写标量目标函数.
创建约束函数文件
生成非线性约束函数,并将其格式化。
文件名= [currdir,“democonstr.m”];gradc matlabFunction (c, [], [],“文件”文件名,“var”{x},...“输出”,{“c”,“量表”,“gradc”,“gradceq”});
结果文件democonstr.m包含以下行或类似行:
函数[c,测查,gradc, gradceq] = democonstr (in)% DEMOCONSTR% [c, ceq, gradc, gradceq] = democonstr (in1)...: x101 = in(28日);...c = [t417。^ 2 +...];如果Nargout > 1 ceq = [];结束如果Nargout > 2 gradc = [2.*x11,...];结束如果Nargout > 3 gradceq = [];结束
该函数具有具有梯度的约束函数的正确形式;看到非线性约束.
生成黑森文件
为了生成问题的拉格朗日的黑森,首先生成能量黑森和约束黑森的文件。
目标函数的黑森函数,hessenergy,是一个包含超过150000个符号的大符号表达式,由求值可知大小(char (hessenergy)).因此,运行matlabFunction (hessenergy)这需要大量的时间。
生成文件hessenergy.m.
文件名= [currdir,“hessenergy.m”];matlabFunction (hessenergy“文件”文件名,“var”, {x});
相比之下,约束函数的hessian值小,计算速度快。
为I = 1:10 ii = num2str(I);name = [“hessc”二);文件名= [currdir, name,“m”];我matlabFunction (hessc {},“文件”文件名,“var”, {x});结束
在为目标和约束生成所有文件之后,将它们与辅助函数中的适当拉格朗日乘子放在一起hessfinal.m,其代码出现在此示例结束.
运行优化
开始优化时,电子随机分布在半径为1/2、中心为[0,0,-1]的球体上。
rng默认的%的再现性Xinitial = randn (10);列是普通的三维向量为j = 1:10 Xinitial (:, j) = Xinitial (:, j) /规范(Xinitial (:, j)) / 2;这归一化为1/2球结束Xinitial(3,:) = Xinitial(3,:) - 1;%在[0,0,-1]中心Xinitial = Xinitial (:);转换为列向量
设置选项以使用内点算法,并使用梯度和黑森。
选择= optimoptions (@fmincon,“算法”,“内点”,“SpecifyObjectiveGradient”,真的,...“SpecifyConstraintGradient”,真的,“HessianFcn”@hessfinal,“显示”,“最后一次”);
调用fmincon.
[xfinal, fval exitflag、输出]= fmincon (@demoenergy Xinitial,...A、b ,[],[],[],[],@ democonstr选项);
找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少,优化完成,在最优性公差的值内,约束满足在约束公差的值内。> <停止标准细节
查看目标函数值、退出标志、迭代次数和函数计算次数。
disp (fval)
34.1365
disp (exitflag)
1
disp (output.iterations)
19
disp (output.funcCount)
28
尽管电子的初始位置是随机的,但最终的位置几乎是对称的。
为我= 1:10 plot3 (xfinal(3 *我2),xfinal(3张),xfinal(3 *我),“r”。,“MarkerSize”25);结束
要检查整个3d几何图形,请旋转图形。
rotate3d
图(手)
与无梯度和黑森优化相比
梯度和Hessians的使用使优化运行得更快、更准确。要比较相同的优化不使用梯度或Hessian信息,设置选项为不使用梯度和Hessian。
选择= optimoptions (@fmincon,“算法”,“内点”,...“显示”,“最后一次”);[xfinal2, fval2 exitflag2 output2] = fmincon (@demoenergy Xinitial,...A、b ,[],[],[],[],@ democonstr选项);
找到目标函数值较低的可行点。找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少,优化完成,在最优性公差的值内,约束满足在约束公差的值内。> <停止标准细节
查看目标函数值、退出标志、迭代次数和函数计算次数。
disp (fval2)
34.1365
disp (exitflag2)
1
disp (output2.iterations)
77
disp (output2.funcCount)
2434
比较有和没有导数信息的迭代次数和函数求值次数。
台=表([output.iterations; output2.iterations], [output.funcCount; output2.funcCount],...“RowNames”,{“衍生品”,“没有衍生品”},“VariableNames”,{“迭代”,函数宏指令的})
台=2×2表迭代Fevals __________ ______有导数19 28无导数77 2434
清除符号变量假设
本例中的符号变量假设它们在符号引擎工作空间中是真实的。删除变量并不能从符号引擎工作区中清除这个假设。使用清除变量假设信谊.
信谊x
验证假设是否为空。
假设(x)
ans =空sym: 1 × 0
Helper函数
下面的代码创建hessfinalhelper函数。
函数H = hessfinal (X,λ)%调用函数hessenergy来启动H = hessenergy (X);加拉格朗日乘子*约束黑森。H = H + hessc1(X) * lambda.ineqnonlin(1);H = H + hessc2(X) * lambda.ineqnonlin(2);H = H + hessc3(X) * lambda.ineqnonlin(3);H = H + hessc4(X) * lambda.ineqnonlin(4);H = H + hessc5(X) * lambda.ineqnonlin(5);H = H + hessc6(X) * lambda.ineqnonlin(6);H = H + hessc7(X) * lambda.ineqnonlin(7);H = H + hessc8(X) * lambda.ineqnonlin(8);H = H + hessc9(X) * lambda.ineqnonlin(9);H = H + hessc10(X) * lambda.ineqnonlin(10);结束
这个例子展示了如何获得更快和更健壮的解决非线性优化问题使用 与优化工具箱™函数相比,符号数学工具箱函数具有不同的语法和结构。特别是,符号变量是实数或复标量,而优化工具箱函数传递向量参数。因此,您需要采取几个步骤,以一种适合于的内点算法的形式符号化地生成目标函数、约束条件及其所有必要的导数 基于问题的优化可以自动计算和使用梯度;看到 考虑将10个电子放入导电体的静电问题。受限于处于体内,电子以一种使其总势能最小的方式排列自己。所有的电子都在物体的边界上。电子是不可区分的,所以这个问题不存在唯一的最小值(排列一个解中的电子会给出另一个有效解)。这个例子的灵感来自多兰、Moré和曼森 这个例子是一个由不平等定义的指导机构
. 这个身体看起来像球体上的金字塔。 这幅画的上下表面之间有一点缝隙。这个间隙是用于创建图形的一般绘图例程的产物。这个程序擦除一个表面上与另一个表面接触的任何矩形补丁。 生成一个符号向量 显示向量
写出线性约束条件
每一个电子都有四个线性不等式: 这个问题总共有40个线性不等式。 用结构化的方式写出不等式。 你们可以看到
非线性约束也是结构化的。
. 生成约束及其梯度和hessian。 约束向量 海赛矩阵, 使用 目标函数,势能,是每对电子之间距离的倒数之和。
. 距离是矢量分量差的平方和的平方根。 计算能量(目标函数)及其梯度和黑森函数。 目标函数有两个输出, 结果文件 该函数具有具有梯度的目标函数的正确形式;看到 生成非线性约束函数,并将其格式化。 结果文件 该函数具有具有梯度的约束函数的正确形式;看到 为了生成问题的拉格朗日的黑森,首先生成能量黑森和约束黑森的文件。 目标函数的黑森函数, 生成文件 相比之下,约束函数的hessian值小,计算速度快。 在为目标和约束生成所有文件之后,将它们与辅助函数中的适当拉格朗日乘子放在一起 开始优化时,电子随机分布在半径为1/2、中心为[0,0,-1]的球体上。 设置选项以使用 调用 查看目标函数值、退出标志、迭代次数和函数计算次数。 尽管电子的初始位置是随机的,但最终的位置几乎是对称的。 要检查整个3d几何图形,请旋转图形。 梯度和Hessians的使用使优化运行得更快、更准确。要比较相同的优化不使用梯度或Hessian信息,设置选项为不使用梯度和Hessian。 查看目标函数值、退出标志、迭代次数和函数计算次数。 比较有和没有导数信息的迭代次数和函数求值次数。 本例中的符号变量假设它们在符号引擎工作空间中是真实的。删除变量并不能从符号引擎工作区中清除这个假设。使用清除变量假设 验证假设是否为空。 下面的代码创建
雅可比矩阵
(符号数学工具箱)matlabFunction
(符号数学工具箱)问题描述
(X, Y) = meshgrid (1: .01:1);Z1 = -abs(X) -abs(Y);Z2 = -1 -√(1 - x ^2 - y ^2);Z2 =实际(Z2);W1 = Z1;W2 = Z2;W1(Z1 < Z2) = nan;
创建变量
X = cell(3,10);
x
x =
包括线性约束
Xi3 - xi1 - xi2≤0 Xi3 - xi1 + xi2≤0 Xi3 + xi1 - xi2≤0 Xi3 + xi1 + xi2≤0
B = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1);一个= 0(40岁,30);
disp (A * x)
建立非线性约束及其梯度和黑森
c =符号(0 (10));我= 1:10;c = (x(3 *我2)。^ 2 + x(3张)。^ 2 +(x(3*i)+1).^2 - 1).'; gradc = jacobian(c,x).';% .'执行转置
建立目标函数及其梯度和黑森函数
能量=信谊(0);
创建目标函数文件
curdir = [pwd filesep];
matlabFunction
函数
创建约束函数文件
文件名= [currdir,
函数
生成黑森文件
文件名= [currdir,
为
运行优化
rng
选择= optimoptions (@fmincon,
[xfinal, fval exitflag、输出]= fmincon (@demoenergy Xinitial,
找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少,优化完成,在最优性公差的值内,约束满足在约束公差的值内。> <停止标准细节
disp (fval)
34.1365
disp (exitflag)
1
disp (output.iterations)
19
disp (output.funcCount)
28
为
rotate3d
图(手)
与无梯度和黑森优化相比
选择= optimoptions (@fmincon,
找到目标函数值较低的可行点。找到了满足约束条件的局部极小值。由于目标函数在可行方向上不减少,优化完成,在最优性公差的值内,约束满足在约束公差的值内。> <停止标准细节
disp (fval2)
34.1365
disp (exitflag2)
1
disp (output2.iterations)
77
disp (output2.funcCount)
2434
台=表([output.iterations; output2.iterations], [output.funcCount; output2.funcCount],
台=
清除符号变量假设
信谊
假设(x)
ans =空sym: 1 × 0
Helper函数
函数