線形等式制約を使用した最小化 (信頼領域 Reflective 法アルゴリズム)

fmincon trust-region-reflective アルゴリズムでは、線形等式制約のみに従って (範囲または他の制約なしに) 非線形目的関数を最小化できます。たとえば、次を最小化します。

$f (x) = \sum_{i = 1}^{n - 1} ({(x_{i}^{2})}^{(x_{i + 1}^{2} + 1)} + {(x_{i + 1}^{2})}^{(x_{i}^{2} + 1)}),$

いくつかの線形等式制約に従います。この例では $n = 1000$ を取ります。

問題の作成

browneq.mat ファイルには、線形制約 Aeq*x = beq を表す行列 Aeq および beq が含まれます。Aeq 行列には、100 個の線形制約を表す 100 行が含まれています (つまり Aeq は 100 行 1,000 列の行列)。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)

disp(fval)

  205.9313

disp(output.constrviolation)

   2.2604e-13

終了フラグの値が 3 であるということは、目的関数の値が許容誤差 FunctionTolerance より小さく変動したために fmincon が停止したことを示します。最終目的関数値は 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