Main Content

Algoritmo interior-point de fmincon con matriz hessiana analítica

Este ejemplo muestra cómo utilizar la información de la derivada para hacer que el proceso de resolución sea más rápido y robusto. El algoritmo interior-point de fmincon puede aceptar una función de matriz hessiana como una entrada. Cuando proporciona una matriz hessiana, puede obtener una solución más rápida y precisa a un problema de minimización restringido.

La función auxiliar bigtoleft es una función objetivo que se hace negativa rápidamente a medida que la coordenada x(1) se convierte en negativa. Su gradiente es un vector de tres elementos. El código para la función auxiliar bigtoleft aparece al final de este ejemplo.

La restricción establecida para este ejemplo es la intersección de los interiores de dos conos, uno que señala hacia arriba y el otro, hacia abajo. La función de restricción es un vector de dos componentes que contiene un componente para cada cono. Puesto que este ejemplo es tridimensional, el gradiente de la restricción es una matriz de 3 por 2. El código para la función auxiliar twocone aparece al final de este ejemplo.

Cree una figura de las restricciones, coloreada utilizando la función objetivo.

% 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.

Crear una función de matriz hessiana

Para utilizar información de la derivada de segundo orden en el solver fmincon, debe crear una matriz hessiana que sea la matriz hessiana del lagrangiano. La ecuación proporciona la matriz hessiana del lagrangiano

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

En este caso, f(x) es la función bigtoleft y ci(x) son las funciones de restricción de dos conos. La función auxiliar hessinterior al final de este ejemplo calcula la matriz hessiana del lagrangiano en un punto x con la estructura de multiplicadores de Lagrange lambda. La primera función calcula 2f(x). Después, calcula las dos matrices hessianas de restricción 2c1(x) y 2c2(x), las multiplica por sus correspondientes multiplicadores de Lagrange lambda.ineqnonlin(1) y lambda.ineqnonlin(2) y las suma. En la definición de la función de restricción twocone puede ver que 2c1(x)=2c2(x), lo que simplifica el cálculo.

Crear opciones para utilizar derivadas

Para permitir que fmincon utilice el gradiente objetivo, gradientes de restricción y la matriz hessiana, debe establecer las opciones adecuadas. La opción HessianFcn que utiliza la matriz hessiana del lagrangiano solo está disponible para el algoritmo interior-point.

options = optimoptions('fmincon','Algorithm','interior-point',...
    "SpecifyConstraintGradient",true,"SpecifyObjectiveGradient",true,...
    'HessianFcn',@hessinterior);

Minimizar utilizando toda la información de las derivadas

Establezca el punto inicial x0 = [-1,-1,-1].

x0 = [-1,-1,-1];

El problema no tiene restricciones lineales ni límites. Establezca esos argumentos en [].

A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];

Resuelva el problema.

[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.

Examinar la solución y el proceso de resolución

Examine la solución, el valor de la función objetivo, el indicador de salida y el número de evaluaciones de función y de iteraciones.

disp(x)
   -6.5000   -0.0000   -3.5000
disp(fval)
  -2.8941e+03
disp(eflag)
     1
disp([output.funcCount,output.iterations])
     7     6

Si no utiliza una función de matriz hessiana, fmincon necesita más iteraciones para convergir y requiere más evaluaciones de función.

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

Si tampoco incluye la información de gradiente, fmincon necesita el mismo número de iteraciones, pero requiere muchas más evaluaciones de función.

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

Funciones auxiliares

Este código crea la función auxiliar 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

Este código crea la función auxiliar 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

Este código crea la función auxiliar 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

Temas relacionados