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

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

где $$p\mathrm{=}\mathrm{}$$7/3$${,}_{\mathrm{}}\mathrm{x0}$$= 0 $${}_{и\mathrm{}}\mathrm{xn}$$+ 1 = 0. tbroyfg вспомогательная функция в конце этого примера реализует целевую функцию, включая ее градиент.

Задача имеет границы $$-\mathrm{}\mathrm{}{}_{}\mathrm{}\mathrm{}$$10≤xi≤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