Создайте функцию гессиана

Использовать производную информацию второго порядка в fmincon решатель, необходимо создать Гессиан, который является Гессианом функции Лагранжа. Гессиан функции Лагранжа дан уравнением

Здесь, $$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 опция с помощью Гессиана функции Лагранжа доступна только для алгоритма внутренней точки.

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