Main Content

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

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

f(x)=i=1n-1((xi2)(xi+12+1)+(xi+12)(xi2+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)
     3
disp(fval)
  205.9313
disp(output.constrviolation)
   2.2427e-13

終了フラグの値が 3 であるということは、目的関数の値が許容誤差 FunctionTolerance より小さく変動したために fmincon が停止したことを示します。最終目的関数値は fval によって与えられます。output.constrviolation の数値は非常に小さく、制約が満たされていることがわかります。

自分で制約違反を計算するには、次のコードを実行します。

norm(Aeq*x-beq,Inf)
ans = 2.2427e-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

関連するトピック