fmincon Алгоритм внутренней точки с аналитическим гессианом

Скрипт Open Live Script

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

Функция помощника bigtoleft целевая функция, которая становится быстро отрицательной как x(1) координата становится отрицательной. Его градиент является трехэлементным вектором. Код для bigtoleft функция помощника появляется в конце этого примера.

Ограничительное множество для этого примера является пересечением внутренних частей двух конусов — один, обращенный вверх и один, обращенный вниз. Ограничительная функция является двухкомпонентным вектором, содержащим один компонент для каждого конуса. Поскольку этот пример 3D, градиент ограничения является 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)+∑λi∇2ci(x)+∑λi∇2ceqi(x).

Здесь, f(x) bigtoleft функция, и ci(x) две конических ограничительных функции. hessinterior функция помощника в конце этого примера вычисляет Гессиан функции Лагранжа в точке x со структурой множителя Лагранжа lambda. Функция сначала вычисляет ∇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

Похожие темы

Документация Optimization Toolbox
Поддержка
Для просмотра документации необходимо авторизоваться на сайте