해석적 헤세 행렬을 사용하는 fmincon Interior-Point 알고리즘

이 예제에서는 도함수 정보를 사용하여 풀이 과정을 좀 더 빠르고 견고하게 만드는 방법을 보여줍니다. fmincon interior-point 알고리즘은 헤세 행렬 함수를 입력값으로 받을 수 있습니다. 헤세 행렬을 제공하는 경우 제약 조건이 있는 최소화 문제에 대해 더 정확한 해를 더 빠르게 얻을 수 있습니다.

헬퍼 함수 bigtoleft는 x(1) 좌표가 음수가 되면 음의 방향으로 빠르게 증가하는 목적 함수입니다. 기울기는 요소를 3개 가진 벡터입니다. bigtoleft 헬퍼 함수의 코드는 이 예제의 마지막 부분에 나와 있습니다.

이 예제에 대한 제약 조건 세트는 하나는 위쪽으로 뾰족하고 하나는 아래쪽으로 뾰족한 모양의 두 원뿔 내부의 교점입니다. 제약 조건 함수는 두 개 성분을 갖는 벡터로, 각 원뿔에 대해 성분이 하나씩 있습니다. 이 예제는 3차원이므로 제약 조건의 기울기는 3×2 행렬입니다. twocone 헬퍼 함수의 코드는 이 예제의 마지막 부분에 나와 있습니다.

목적 함수를 사용하여 색을 지정하는 제약 조건을 적용한 Figure를 생성합니다.

% Create figure
figure1 = figure;

% Create axes
axes1 = axes('Parent',figure1);
view([-63.5 18]);
grid('on');
hold('on');

% 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 솔버에서 2계 도함수 정보를 사용하려면 헤세 행렬을 라그랑주의 헤세 행렬로 생성해야 합니다. 라그랑주의 헤세 행렬은 다음 방정식으로 주어집니다.

$${\nabla }_{xx}^{2}L(x,\lambda )={\nabla }^{2}f(x)+\sum {\lambda }_i{\nabla }^2{c}_i(x)+\sum {\lambda }_i{\nabla }^{2}ce{q}_i(x).$$

여기서 $$f(x)$$는 bigtoleft 함수이고, $$c_i(x)$$는 두 원뿔의 제약 조건 함수입니다. 이 예제의 마지막 부분에 있는 hessinterior 헬퍼 함수는 점 x에서 라그랑주 승수 구조체 lambda로 라그랑주의 헤세 행렬을 계산합니다. 함수는 먼저 $${\nabla }^{2}f(x)$$를 계산합니다. 그런 다음 두 제약 조건 헤세 행렬 $${\nabla }^2{c}_1(x)$$와 $${\nabla }^2{c}_2(x)$$를 계산하고 그 두 값을 해당하는 라그랑주 승수 lambda.ineqnonlin(1)과 lambda.ineqnonlin(2)에 곱한 다음 값을 더합니다. twocone 제약 조건 함수의 정의를 통해 $${\nabla }^2{c}_1(x)={\nabla }^2{c}_2(x)$$임을 알 수 있습니다. 따라서 계산이 간단해집니다.

도함수를 사용하기 위한 옵션 생성하기

fmincon이 목적 함수 기울기, 제약 조건 기울기 및 헤세 행렬을 사용하도록 설정하려면 관련 옵션을 설정해야 합니다. 라그랑주의 헤세 행렬을 사용하는 HessianFcn 옵션은 interior-point 알고리즘에만 사용할 수 있습니다.

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.

<stopping criteria details>

해와 풀이 과정 검토하기

해, 목적 함수 값, 종료 플래그, 함수 실행 횟수, 반복 횟수를 검토합니다.

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.

<stopping criteria details>

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.

<stopping criteria details>

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

참고 항목

도움말 항목

• 2차 제약 조건이 있는 선형 목적 함수 또는 2차 목적 함수
• Calculate Gradients and Hessians Using Symbolic Math Toolbox