Минимизация с линейными ограничениями равенства, Доверительная область отражающий алгоритм

fmincon trust-region-reflective алгоритм может минимизировать нелинейную целевую функцию, удовлетворяющую линейным ограничениям равенства только (никакие границы или любые другие ограничения). Например, минимизировать

f(x)=i=1n-1((xi2)(xi+12+1)+(xi+12)(xi2+1)),

подвергните некоторым линейным ограничениям равенства. Этот пример берет n=1000.

Создайте проблему

browneq.mat файл содержит матрицы Aeq и beq, которые представляют линейные ограничения   Aeq*x = beq. Aeq матрица имеет 100 строк, представляющих 100 линейных ограничений (так Aeq 100 1000 матрица). Загрузите browneq.mat данные.

load browneq.mat

brownfgh функция помощника в конце этого примера реализует целевую функцию, включая ее градиент и Гессиан.

Установка опций

trust-region-reflective алгоритм требует, чтобы целевая функция включала градиент. Алгоритм принимает Гессиан в целевой функции. Установите опции включать всю производную информацию.

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

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

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

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

Проблема не имеет никаких границ, линейных ограничений неравенства или нелинейных ограничений, таким образом, устанавливает те параметры на [].

A = [];
b = [];
lb = [];
ub = [];
nonlcon = [];

Вызовите fmincon решать задачу.

[x,fval,exitflag,output] = ...
   fmincon(@brownfgh,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)
  205.9313
disp(output.constrviolation)
   2.2604e-13

exitflag значение 3 указывает на тот fmincon остановки, потому что изменение в значении целевой функции меньше допуска FunctionTolerance. Итоговое значение целевой функции дано fval. Ограничениям удовлетворяют, как показано в output.constrviolation, который отображает очень небольшое число.

Чтобы вычислить нарушение ограничений самостоятельно, выполните следующий код.

norm(Aeq*x-beq,Inf)
ans = 2.2604e-13

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

Следующий код создает brownfgh функция помощника.

function [f,g,H] = brownfgh(x)
%BROWNFGH  Nonlinear minimization problem (function, its gradients
% and Hessian).
% Documentation example.         

%   Copyright 1990-2019 The MathWorks, Inc.

% Evaluate the function.
  n = length(x);
  y = zeros(n,1);
  i = 1:(n-1);
  y(i) = (x(i).^2).^(x(i+1).^2+1)+(x(i+1).^2).^(x(i).^2+1);
  f = sum(y);

% Evaluate the gradient.
  if nargout > 1
     i=1:(n-1);
     g = zeros(n,1);
     g(i) = 2*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
             2*x(i).*((x(i+1).^2).^(x(i).^2+1)).*log(x(i+1).^2);
     g(i+1) = g(i+1)+...
              2*x(i+1).*((x(i).^2).^(x(i+1).^2+1)).*log(x(i).^2)+...
              2*(x(i).^2+1).*x(i+1).*((x(i+1).^2).^(x(i).^2));
  end

% Evaluate the (sparse, symmetric) Hessian matrix
  if nargout > 2
     v = zeros(n,1);
     i = 1:(n-1);
     v(i) = 2*(x(i+1).^2+1).*((x(i).^2).^(x(i+1).^2))+...
            4*(x(i+1).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i).^2).^((x(i+1).^2)-1))+...
            2*((x(i+1).^2).^(x(i).^2+1)).*(log(x(i+1).^2));
     v(i) = v(i)+4*(x(i).^2).*((x(i+1).^2).^(x(i).^2+1)).*((log(x(i+1).^2)).^2);
     v(i+1) = v(i+1)+...
              2*(x(i).^2).^(x(i+1).^2+1).*(log(x(i).^2))+...
              4*(x(i+1).^2).*((x(i).^2).^(x(i+1).^2+1)).*((log(x(i).^2)).^2)+...
              2*(x(i).^2+1).*((x(i+1).^2).^(x(i).^2));
     v(i+1) = v(i+1)+4*(x(i).^2+1).*(x(i+1).^2).*(x(i).^2).*((x(i+1).^2).^(x(i).^2-1));
     v0 = v;
     v = zeros(n-1,1);
     v(i) = 4*x(i+1).*x(i).*((x(i).^2).^(x(i+1).^2))+...
            4*x(i+1).*(x(i+1).^2+1).*x(i).*((x(i).^2).^(x(i+1).^2)).*log(x(i).^2);
     v(i) = v(i)+ 4*x(i+1).*x(i).*((x(i+1).^2).^(x(i).^2)).*log(x(i+1).^2);
     v(i) = v(i)+4*x(i).*((x(i+1).^2).^(x(i).^2)).*x(i+1);
     v1 = v;
     i = [(1:n)';(1:(n-1))'];
     j = [(1:n)';(2:n)'];
     s = [v0;2*v1];
     H = sparse(i,j,s,n,n);
     H = (H+H')/2;
  end
end

Похожие темы