Main Content

使用解析黑塞函数的 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 求解器中使用二阶导数信息,您必须创建一个黑塞矩阵,它是拉格朗日乘数的黑塞矩阵。拉格朗日函数的黑塞矩阵由以下方程给出

xx2L(x,λ)=2f(x)+λi2ci(x)+λi2ceqi(x).

此处,f(x)bigtoleft 函数,ci(x) 是两个锥约束函数。位于此示例末尾hessinterior 辅助函数使用拉格朗日乘数结构体 lambda 计算在点 x 处的拉格朗日函数的黑塞矩阵。该函数首先计算 2f(x)。然后,它计算两个约束:黑塞 2c1(x)2c2(x),将它们乘以对应的拉格朗日乘数 lambda.ineqnonlin(1)lambda.ineqnonlin(2),并将它们相加。从 twocone 约束函数的定义可以看出 2c1(x)=2c2(x),这可以简化计算。

设置选项以使用导数

要使 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

相关主题