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

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

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

Проблема имеет границы: $$-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. Когда проблема, столь большая, как это имеет очевидную структуру разреженности, не устанавливая HessPattern опция требует огромной суммы ненужной памяти и расчета. Это вызвано тем, что 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

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

disp(fval2 - fval)

  -2.9005e-07

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

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

Этот код создает tbroy функция.

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