使用解析黑塞函数的 fmincon
内点算法
此示例说明如何使用导数信息来使求解过程更快、更稳健。fmincon
内点算法可以接受黑塞函数作为输入。当您提供黑塞函数时,您可以更快得到约束最小化问题更为准确的解。
辅助函数 bigtoleft
是一个目标函数,随着 x(1)
坐标变为负值,该函数迅速变为负值。它的梯度是三元素向量。bigtoleft
辅助函数的代码出现在此示例末尾。
此示例的约束集是两个圆锥体(一个尖端向上,一个尖端向下)内部的交集。约束函数是一个双分量向量,每个圆锥体包含一个分量。一个由于此示例是三维的,约束的梯度是一个 3×2 矩阵。twocone
辅助函数的代码出现在此示例末尾。
创建一个约束图窗,使用目标函数进行着色。
% Create figure figure1 = figure; % Create axes axes1 = axes('Parent',figure1); view([-63.5 18]); grid('on'); hold('all'); % Set up polar coordinates and two cones r=0:.1:6.5; th=2*pi*(0:.01:1); x=r'*cos(th); y=r'*sin(th); z=-10+sqrt(x.^2+y.^2); zz=3-sqrt(x.^2+y.^2); % Evaluate objective function on cone surfaces newxf=reshape(bigtoleft([x(:),y(:),z(:)]),66,101)/3000; newxg=reshape(bigtoleft([x(:),y(:),z(:)]),66,101)/3000; % Create lower surf with color set by objective surf(x,y,z,newxf,'Parent',axes1,'EdgeAlpha',0.25); % Create upper surf with color set by objective surf(x,y,zz,newxg,'Parent',axes1,'EdgeAlpha',0.25); axis equal
创建黑塞函数
要在 fmincon
求解器中使用二阶导数信息,您必须创建一个黑塞矩阵,它是拉格朗日乘数的黑塞矩阵。拉格朗日函数的黑塞矩阵由以下方程给出
此处, 是 bigtoleft
函数, 是两个锥约束函数。位于此示例末尾的 hessinterior
辅助函数使用拉格朗日乘数结构体 lambda
计算在点 x
处的拉格朗日函数的黑塞矩阵。该函数首先计算 。然后,它计算两个约束:黑塞 和 ,将它们乘以对应的拉格朗日乘数 lambda.ineqnonlin(1)
和 lambda.ineqnonlin(2)
,并将它们相加。从 twocone
约束函数的定义可以看出 ,这可以简化计算。
设置选项以使用导数
要使 fmincon
能够使用目标梯度、约束梯度和黑塞矩阵,您必须设置适当的选项。使用拉格朗日乘数的黑塞矩阵的 HessianFcn
选项仅适用于内点算法。
options = optimoptions('fmincon','Algorithm','interior-point',... "SpecifyConstraintGradient",true,"SpecifyObjectiveGradient",true,... 'HessianFcn',@hessinterior);
使用所有导数信息进行最小化
设置初始点 x0 = [-1,-1,-1]
。
x0 = [-1,-1,-1];
该问题没有线性约束或边界。将这些参量设置为 []
。
A = []; b = []; Aeq = []; beq = []; lb = []; ub = [];
求解。
[x,fval,eflag,output] = fmincon(@bigtoleft,x0,...
A,b,Aeq,beq,lb,ub,@twocone,options);
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
检查解和求解过程
检查解、目标函数值、退出标志以及函数计算和迭代的次数。
disp(x)
-6.5000 -0.0000 -3.5000
disp(fval)
-2.8941e+03
disp(eflag)
1
disp([output.funcCount,output.iterations])
7 6
如果您不使用黑塞函数,fmincon
需要更多次迭代才能收敛,并且需要更多函数计算。
options.HessianFcn = [];
[x2,fval2,eflag2,output2] = fmincon(@bigtoleft,x0,...
A,b,Aeq,beq,lb,ub,@twocone,options);
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
disp([output2.funcCount,output2.iterations])
13 9
如果您也未提供梯度信息,fmincon
会执行相同次数的迭代,但要求进行更多函数计算。
options.SpecifyConstraintGradient = false;
options.SpecifyObjectiveGradient = false;
[x3,fval3,eflag3,output3] = fmincon(@bigtoleft,x0,...
A,b,Aeq,beq,lb,ub,@twocone,options);
Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
disp([output3.funcCount,output3.iterations])
43 9
辅助函数
以下代码会创建 bigtoleft
辅助函数。
function [f gradf] = bigtoleft(x) % This is a simple function that grows rapidly negative % as x(1) becomes negative % f = 10*x(:,1).^3+x(:,1).*x(:,2).^2+x(:,3).*(x(:,1).^2+x(:,2).^2); if nargout > 1 gradf=[30*x(1)^2+x(2)^2+2*x(3)*x(1); 2*x(1)*x(2)+2*x(3)*x(2); (x(1)^2+x(2)^2)]; end end
以下代码会创建 twocone
辅助函数。
function [c ceq gradc gradceq] = twocone(x) % This constraint is two cones, z > -10 + r % and z < 3 - r ceq = []; r = sqrt(x(1)^2 + x(2)^2); c = [-10+r-x(3); x(3)-3+r]; if nargout > 2 gradceq = []; gradc = [x(1)/r,x(1)/r; x(2)/r,x(2)/r; -1,1]; end end
以下代码会创建 hessinterior
辅助函数。
function h = hessinterior(x,lambda) h = [60*x(1)+2*x(3),2*x(2),2*x(1); 2*x(2),2*(x(1)+x(3)),2*x(2); 2*x(1),2*x(2),0];% Hessian of f r = sqrt(x(1)^2+x(2)^2);% radius rinv3 = 1/r^3; hessc = [(x(2))^2*rinv3,-x(1)*x(2)*rinv3,0; -x(1)*x(2)*rinv3,x(1)^2*rinv3,0; 0,0,0];% Hessian of both c(1) and c(2) h = h + lambda.ineqnonlin(1)*hessc + lambda.ineqnonlin(2)*hessc; end