exponenta event banner

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

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

f (x) =∑i=1n-1 ((x2) (xi + 12 + 1) + (xi + 12) (x2 + 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

Связанные темы