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

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

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

Проблема имеет границы $$-10\le x_i\le 10$$ \forall $$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