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);

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

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

制約を含める

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

load browneq

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

prob.Constraints = Aeq*x == beq;

問題の確認と求解

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

show(prob.Objective)
  brownfgh(x)

この問題には 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: 2.3698e-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: 32
               funcCount: 33065
         constrviolation: 2.4869e-14
                stepsize: 7.0542e-06
               algorithm: 'interior-point'
           firstorderopt: 4.2570e-06
            cgiterations: 3
                 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 を使用した勾配とヘッシアンの計算を参照してください。

参考

関連するトピック