Целевая функция с градиентом

Для данного $$n$$ что является положительным, кратным 4, целевая функция является

где $$p=7/3$$, $$x_0=0$$, и $$x_{n+1}=0$$. The tbroyfg Функция helper в конце этого примера реализует целевую функцию, включая ее градиент.

Задача имеет границы $$-10\le x_i\le 10$$ для всех $$i$$. Использовать $$n$$ = 800.

n = 800;
lb = -10*ones(n,1);
ub = -lb;

Гессианский шаблон

Шаблон разреженности матрицы Гессия предопределен и сохранен в файле tbroyhstr.mat. Структура разреженности для Гессиана этой задачи перемычка, как видно на следующем шпионском графике.

load tbroyhstr
spy(Hstr)

На этом графике центральная полоса сама является пятиполосной матрицей. Следующий график показывает матрицу более четко.

spy(Hstr(1:20,1:20))

Опции задачи

Установите опции, чтобы использовать trust-region-reflective алгоритм. Этот алгоритм требует, чтобы вы установили SpecifyObjectiveGradient опция для true.

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

options = optimoptions('fmincon','SpecifyObjectiveGradient',true,'HessPattern',Hstr,...
    'Algorithm','trust-region-reflective');

Решите задачу

Установите начальную точку -1 для нечетных индексов и + 1 для четных индексов.

x0 = -ones(n,1);
x0(2:2:n) = 1;

Задача не имеет линейных или нелинейных ограничений, поэтому установите эти параметры равными [].

A = [];
b = [];
Aeq = [];
beq = [];
nonlcon = [];

Функции fmincon чтобы решить проблему.

[x,fval,exitflag,output] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

Local minimum possible.

fmincon stopped because the final change in function value relative to 
its initial value is less than the value of the function tolerance.

Исследуйте процесс решения и решения

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

disp(exitflag);

     3

disp(fval)

  270.4790

disp(output.firstorderopt)

    0.0163

disp(output.iterations)

     7

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

Улучшите решение

Попробуйте использовать пятиполосный предкондиционер вместо диагонального предкондиционера по умолчанию. Использование optimoptions, установите PrecondBandWidth опция 2 и снова решить проблему. (Шумовая полоса - это количество верхних или нижних диагоналей, не считая основной диагонали.)

options.PrecondBandWidth = 2;
[x2,fval2,exitflag2,output2] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

Local minimum possible.

fmincon stopped because the final change in function value relative to 
its initial value is less than the value of the function tolerance.

disp(exitflag2);

     3

disp(fval2)

  270.4790

disp(output2.firstorderopt)

   7.5340e-05

disp(output2.iterations)

     9

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

disp(fval2 - fval)

  -2.9005e-07

Значение целевой функции уменьшается на крошечную величину. Решение в основном улучшает меру оптимальности первого порядка, а не целевую функцию.

Функция помощника

Этот код создает tbroyfg вспомогательная функция.

function [f,grad] = tbroyfg(x,dummy)
%TBROYFG Objective function and its gradients for nonlinear minimization
% with bound constraints and banded preconditioner.
% Documentation example.

%   Copyright 1990-2008 The MathWorks, Inc.

n = length(x);  % n should be a multiple of 4
p = 7/3;
y=zeros(n,1);
i = 2:(n-1);
y(i) = abs((3-2*x(i)).*x(i) - x(i-1) - x(i+1)+1).^p;
y(n) = abs((3-2*x(n)).*x(n) - x(n-1)+1).^p;
y(1) = abs((3-2*x(1)).*x(1) - x(2)+1).^p;
j = 1:(n/2);
z = zeros(length(j),1);
z(j) = abs(x(j) + x(j+n/2)).^p;
f = 1 + sum(y) + sum(z);
%
% Evaluate the gradient.
if nargout > 1
   g = zeros(n,1);
   t = zeros(n,1);
   i = 2:(n-1);
   t(i) = (3-2*x(i)).*x(i) - x(i-1) - x(i+1) + 1;
   g(i) = p*abs(t(i)).^(p-1).*sign(t(i)).*(3-4*x(i));
   g(i-1) = g(i-1) - p*abs(t(i)).^(p-1).*sign(t(i));
   g(i+1) = g(i+1) - p*abs(t(i)).^(p-1).*sign(t(i));
   tt = (3-2*x(n)).*x(n) - x(n-1) + 1;
   g(n) = g(n) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(n));
   g(n-1) = g(n-1) - p*abs(tt).^(p-1).*sign(tt);
   tt = (3-2*x(1)).*x(1)-x(2)+1;
   g(1) = g(1) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(1));
   g(2) = g(2) - p*abs(tt).^(p-1).*sign(tt);
   j = 1:(n/2);
   t(j) = x(j) + x(j+n/2);
   g(j) = g(j) + p*abs(t(j)).^(p-1).*sign(t(j));
   jj = j + (n/2);
   g(jj) = g(jj) + p*abs(t(j)).^(p-1).*sign(t(j));
   grad = g;
end
end