计算梯度和黑森使用符号数学工具箱
这个例子展示了如何获得更快和更健壮的非线性优化问题的解决方案使用金宝搏官方网站fmincon以及符号数学工具箱™功能。如果你有Symbolic Math Toolbox许可证,你可以使用这些Symbolic Math Toolbox函数轻松地计算目标函数和约束函数的分析梯度和Hessians:
雅可比矩阵
(符号数学工具箱)生成一个标量函数的梯度,并生成一个向量函数偏导数的矩阵。例如,你可以通过应用得到黑森矩阵(目标函数的二阶导数)雅可比矩阵对梯度。这个例子展示了如何使用雅可比矩阵生成目标函数和约束函数的符号梯度和黑森函数。
matlabFunction
(符号数学工具箱)生成匿名函数或计算符号表达式值的文件。这个例子展示了如何使用matlabFunction生成评估目标函数和约束函数及其任意点导数的文件。
与优化工具箱™函数相比,符号数学工具箱函数具有不同的语法和结构。特别是,符号变量是实数或复数标量,而“优化工具箱”函数传递矢量参数。因此,您需要采取几个步骤,以一种适合的内点算法的形式,象征性地生成目标函数、约束及其所有必要的导数fmincon。
基于问题的优化可以自动计算和使用梯度;看到优化工具箱中的自动区分。有关使用自动微分来解决此问题的基于问题的方法,请参见基于优化变量的约束静电非线性优化。
问题描述
考虑在导电体中放置10个电子的静电问题。电子以一种使它们的总势能最小化的方式排列自己,受限于躺在体内。所有的电子都在物体的边界上。电子是不可区分的,所以这个问题不存在唯一的最小值(在一个解中排列电子会给出另一个有效的解)。这个例子的灵感来自Dolan, Moré和Munson[58]。
这个例子是由不等式定义的传导体
。
这个身体看起来像球体上的金字塔。
[X,Y] = meshgrid(-1:.01:1);Z1 = -abs(X) -abs(Y);Z2 = -1 -根号(1 - x ^2 - y ^2);Z2 = real(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之间,和j在1到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 = sym(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];A = 0 (40,30);为i=1:10 A(4*i-3:4*i,3*i-2:3*i) = B;结束B = 0 (40,1);
你们可以看到A*x≤b表示不等式。
disp (A * x)
创建非线性约束及其梯度和黑森
非线性约束也是结构化的。
。
生成约束及其梯度和黑森。
C = sym(0 (1,10));I = 1:10;C = (x(3*i-2)。^2 + x(3*i-1)^ 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。
使用。'转置的语法。的'句法是指共轭转置,它有不同的符号导数。
创建目标函数及其梯度和黑森
目标函数,势能,是每个电子对之间距离的倒数的和。
。
距离是矢量分量之差平方和的平方根。
计算能量(目标函数)及其梯度和黑森。
能量= sym(0);为I = 1:3:25为J = i+3:3:28 dist = x(i:i+2) - x(J: J +2);能量=能量+ 1/平方根(dist.'*dist);结束结束Gradenergy =雅可比矩阵(能量,x).';hesessenergy =雅可比矩阵(gradenergy,x);
创建目标函数文件
目标函数有两个输出,能源和gradenergy。调用时将两个函数放在一个向量中matlabFunction减少子表达式的数量matlabFunction生成,并且仅当调用函数(fmincon在本例中)请求两个输出。您可以将生成的文件放在MATLAB路径上的任何文件夹中。在这种情况下,将它们放在当前文件夹中。
Currdir = [pwd filesep];您可能需要使用currdir = pwd文件名= [currdir,“demoenergy.m”];matlabFunction(能源、gradenergy“文件”文件名,“var”, {x});
matlabFunction返回能源作为第一个输出和gradenergy作为第二个。该函数也接受一个输入向量{x}而不是输入列表x11、……x103。
生成的文件demoenergy.m包含以下行或类似行:
函数[energy,gradenergy] = demoenergy(in1)% DEMOENERGY% [energy, gradenergy] = demoenergy (in1)…X101 = in1(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,ceq,gradc,gradceq] = democonstr(in1)% DEMOCONSTR% [c, ceq, gradc, gradceq] = democonstr (in1)…X101 = in1(28,:);…C = [t417]。^ 2 +…];如果Nargout > 1 ceq = [];结束如果Nargout > 2 gradc = [2.*x11,…];结束如果Nargout > 3 gradceq = [];结束
该函数具有具有梯度的约束函数的正确形式;看到非线性约束。
生成黑森文件
为了生成问题的拉格朗日的黑森,首先生成能量黑森和约束黑森的文件。
目标函数的黑森函数,hessenergy,是一个包含超过150,000个符号的大型符号表达式,由求值可知大小(char (hessenergy))。因此,运行matlabFunction (hessenergy)这需要大量的时间。
生成文件hessenergy.m。
文件名= [currdir,“hessenergy.m”];matlabFunction (hessenergy“文件”文件名,“var”, {x});
相比之下,约束函数的黑森值较小,计算速度快。
为I = 1:10 ii = num2str(I);姓名= [“hessc”二);文件名= [currdir, name,“m”];我matlabFunction (hessc {},“文件”文件名,“var”, {x});结束
在为目标和约束生成所有文件之后,将它们与helper函数中适当的拉格朗日乘法器放在一起hessfinal.m,其代码出现在本例结束。
运行优化
从电子随机分布在半径为1/2,中心为[0,0,-1]的球体上开始优化。
rng默认的%用于再现性Xinitial = randn(3,10);列是普通的3-D向量为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,output] = 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);结束
若要检查整个3-D几何图形,请旋转图形。
rotate3d
图(手)
比较优化没有梯度和黑森
梯度和黑森的使用使优化运行得更快更准确。要比较使用无梯度或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
比较有和没有导数信息的迭代次数和函数计算次数。
tbl = table([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,lambda)%调用函数hessenergy启动H = hesessenergy (X);添加拉格朗日乘法器*约束黑森H = H + hessc1(X) * lambda.ineqnonlin(1);H = H + hessc2(X) * lambda.ineqnonlin(2);H = H + hesc3 (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个电子的静电问题。电子以一种使它们的总势能最小化的方式排列自己,受限于躺在体内。所有的电子都在物体的边界上。电子是不可区分的,所以这个问题不存在唯一的最小值(在一个解中排列电子会给出另一个有效的解)。这个例子的灵感来自Dolan, Moré和Munson 这个例子是由不等式定义的传导体
。 这个身体看起来像球体上的金字塔。 在人物的上下表面之间有轻微的间隙。此间隙是用于创建图形的一般绘图例程的工件。该例程擦除一个表面上与另一个表面接触的任何矩形补丁。 生成一个符号向量 显示矢量
写出线性约束
每一个电子都有四个线性不等式: 这个问题总共有40个线性不等式。 用结构化的方式写出不等式。 你们可以看到
非线性约束也是结构化的。
。 生成约束及其梯度和黑森。 约束向量 黑森矩阵, 使用 目标函数,势能,是每个电子对之间距离的倒数的和。
。 距离是矢量分量之差平方和的平方根。 计算能量(目标函数)及其梯度和黑森。 目标函数有两个输出, 生成的文件 该函数具有具有梯度的目标函数的正确形式;看到 生成非线性约束函数,并将其转换成正确的格式。 生成的文件 该函数具有具有梯度的约束函数的正确形式;看到 为了生成问题的拉格朗日的黑森,首先生成能量黑森和约束黑森的文件。 目标函数的黑森函数, 生成文件 相比之下,约束函数的黑森值较小,计算速度快。 在为目标和约束生成所有文件之后,将它们与helper函数中适当的拉格朗日乘法器放在一起 从电子随机分布在半径为1/2,中心为[0,0,-1]的球体上开始优化。 设置选项以使用 调用 查看目标函数值、退出标志、迭代次数和函数计算次数。 尽管电子的初始位置是随机的,但最终位置几乎是对称的。 若要检查整个3-D几何图形,请旋转图形。 梯度和黑森的使用使优化运行得更快更准确。要比较使用无梯度或Hessian信息的相同优化,请将选项设置为不使用梯度和Hessian。 查看目标函数值、退出标志、迭代次数和函数计算次数。 比较有和没有导数信息的迭代次数和函数计算次数。 本例中的符号变量假定它们在符号引擎工作区中是实数。删除变量并不能从符号引擎工作区中清除这个假设。通过使用清除变量假设 验证假设是否为空。 下面的代码创建
雅可比矩阵
(符号数学工具箱)matlabFunction
(符号数学工具箱)问题描述
[X,Y] = meshgrid(-1:.01:1);Z1 = -abs(X) -abs(Y);Z2 = -1 -根号(1 - x ^2 - y ^2);Z2 = real(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];A = 0 (40,30);
disp (A * x)
创建非线性约束及其梯度和黑森
C = sym(0 (1,10));I = 1:10;C = (x(3*i-2)。^2 + x(3*i-1)^ 2 +(x(3*i)+1).^2 - 1).'; gradc = jacobian(c,x).';% .'执行转置
创建目标函数及其梯度和黑森
能量= sym(0);
创建目标函数文件
Currdir = [pwd filesep];
matlabFunction
函数
创建约束函数文件
文件名= [currdir,
函数
生成黑森文件
文件名= [currdir,
为
运行优化
rng
选项= optimoptions(@fmincon,
[xfinal,fval,exitflag,output] = 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
tbl = table([output.iterations;output2.iterations],[output.funcCount;output2.funcCount],
台=
清除符号变量假设
信谊
假设(x)
ans =空sym: 1乘0
Helper函数
函数