Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

問題ベースでの線形制約付き非線形最小化

この例では、問題ベースのアプローチを使用して線形等式制約に従って非線形関数を最小化する方法を説明します。ここでは、最適化変数で制約を定式化します。この例では、fcn2optimexpr を使用して目的関数ファイルを最適化式に変換する方法も説明します。

線形等式制約を使用した最小化 (信頼領域 Reflective 法アルゴリズム)の例では、勾配とヘッシアンを使用するソルバーベースのアプローチを使用します。問題ベースのアプローチを使用すると同じ問題を簡単に解くことができますが、問題ベースのアプローチでは現在、サポートされない演算子が問題に含まれる場合は勾配情報が使用されないため、求解にかかる時間が長くなります。また、問題ベースのアプローチではヘッセ情報をサポートしていません。

問題および目的の作成

問題は、次を最小化することです。

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

一連の線形等式制約 Aeq*x = beq に従います。まず、最適化問題と変数を作成します。

prob = optimproblem;
N = 1000;
x = optimvar('x',N);

目的関数は、この例を実行する際に使用できる brownfgh.m ファイルにあります。最適化変数は指数内に出現しないため、fcn2optimexpr を使用して関数を最適化式に変換しなければなりません。詳細については、最適化変数および式でサポートされる演算非線形関数から最適化式への変換を参照してください。

prob.Objective = fcn2optimexpr(@brownfgh,x,'OutputSize',[1,1]);

制約を含める

ワークスペースに Aeq および beq 行列を取得するには、次のコマンドを実行します。

load browneq

線形制約を問題に含めます。

prob.Constraints = Aeq*x == beq;

問題の確認と求解

問題の目的関数を確認します。

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

%   Copyright 1990-2008 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

この問題には 100 個の線形等式制約があるため、結果として得られる制約式は長すぎて例に含めることができません。制約を表示するには、コメントを解除して次の行を実行します。

% show(prob.Constraints)

初期点をフィールド x を持つ構造体として設定します。

x0.x = -ones(N,1);
x0.x(2:2:N) = 1;

solve を呼び出して問題を解きます。

[sol,fval,exitflag,output] = solve(prob,x0)
Solving problem using fmincon.

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.
sol = struct with fields:
    x: [1000x1 double]

fval = 207.5463
exitflag = 
    SolverLimitExceeded

output = struct with fields:
              iterations: 2
               funcCount: 3007
         constrviolation: 1.6875e-13
                stepsize: 1.9303
               algorithm: 'interior-point'
           firstorderopt: 2.6432
            cgiterations: 0
                 message: 'Solver stopped prematurely....'
            bestfeasible: [1x1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

関数の評価制限を超えるため、ソルバーが途中で停止します。最適化を続行するには、最終点から最適化を再開し、より多くの関数評価を行えるようにします。

options = optimoptions(prob,'MaxFunctionEvaluations',1e5);
[sol,fval,exitflag,output] = solve(prob,sol,'Options',options)
Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
sol = struct with fields:
    x: [1000x1 double]

fval = 205.9313
exitflag = 
    OptimalSolution

output = struct with fields:
              iterations: 31
               funcCount: 32056
         constrviolation: 2.4869e-14
                stepsize: 1.6847e-05
               algorithm: 'interior-point'
           firstorderopt: 6.0482e-06
            cgiterations: 0
                 message: 'Local minimum found that satisfies the constraints....'
            bestfeasible: [1x1 struct]
     objectivederivative: "finite-differences"
    constraintderivative: "closed-form"
                  solver: 'fmincon'

ソルバーベースの解との比較

線形等式制約を使用した最小化 (信頼領域 Reflective 法アルゴリズム)で説明するようにソルバーベースのアプローチを使用して問題を解くには、初期点をベクトルに変換します。その後、brownfgh に指定された勾配とヘッセ情報を使用するようにオプションを設定します。

xstart = x0.x;
fun = @brownfgh;
opts = optimoptions('fmincon','SpecifyObjectiveGradient',true,'HessianFcn','objective',...
    'Algorithm','trust-region-reflective');
[x,fval,exitflag,output] = ...
   fmincon(fun,xstart,[],[],Aeq,beq,[],[],[],opts);
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.
fprintf("Fval = %g\nNumber of iterations = %g\nNumber of function evals = %g.\n",...
    fval,output.iterations,output.funcCount)
Fval = 205.931
Number of iterations = 22
Number of function evals = 23.

線形等式制約を使用した最小化 (信頼領域 Reflective 法アルゴリズム)に示すソルバーベースの解では、目的関数に指定された勾配とヘッシアンが使用されます。ソルバー fmincon は、導関数情報を使用することにより、関数評価を 23 回のみ使用して 22 回の反復で解に収束します。ソルバーベースの解の最終目的関数値は、この問題ベースの解と同じになります。

ただし、シンボリック数式を使用せずに勾配関数とヘッセ関数を作成するのは困難で、エラーが発生しやすくなります。シンボリック数式を使用して導関数を計算する方法を示す例については、Symbolic Math Toolbox を使用した勾配とヘッシアンの計算を参照してください。

参考

関連するトピック