Main Content

解析的ヘッシアンを使用した fmincon の内点法アルゴリズム

この例では、導関数情報を使用して、解法プロセスをより迅速でロバストにする方法を説明します。fmincon の内点法アルゴリズムは入力にヘッセ関数を受け入れることができます。ヘッシアンを与えると、制約付き最小化問題は効率的に解くことができ、より正確な解を得ることができます。

補助関数 bigtoleft は、x(1) 座標が負の場合に急速に負になる目的関数です。その勾配は 3 要素のベクトルです。補助関数 bigtoleft のコードは、この例の終わりに掲載しています。

この例の制約集合は 2 つの円錐の内部交差部分です。1 つの円錐は上向きで、もう 1 つの円錐は下向きです。制約関数は各円錐に対して 1 つの要素を含む 2 要素のベクトルです。この例は 3 次元であるため、制約の勾配は 3 行 2 列の行列です。補助関数 twocone のコードは、この例の終わりに掲載しています。

この制約の Figure を作成します。色付けには目的関数を使用します。

% 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

Figure contains an axes object. The axes object contains 2 objects of type surface.

ヘッセ関数の作成

fmincon ソルバーで 2 次導関数情報を使用するには、ラグランジュ関数のヘッシアンを作成しなければなりません。ラグランジュ関数のヘッシアンは次の式で与えられます。

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

ここで、f(x) は関数 bigtoleft であり、ci(x) は 2 つの円錐制約関数です。補助関数 hessinterior (この例の終わりに掲載) が、点 x におけるラグランジュ関数のヘッシアンをラグランジュ乗数構造体 lambda を使用して計算します。この関数は、まず 2f(x) を計算します。次いで、2 つの制約ヘッシアン 2c1(x) および 2c2(x) を計算してから、これらを対応するラグランジュ乗数 lambda.ineqnonlin(1) および lambda.ineqnonlin(2) で乗算し、その積を加算します。2c1(x)=2c2(x) となることを制約関数 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

関連するトピック