ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

Optimization Toolbox™ ソルバーによる記号数学の使用

この例では、Symbolic Math Toolbox™ の関数 jacobianmatlabFunction を使用して、導関数を最適化ソルバーに与える方法を説明します。一般に、Optimization Toolbox™ ソルバーの精度と効率は、目的関数と制約関数の勾配とヘッセ行列が与えられた場合に向上します。

追加要件:

  • Symbolic Math Toolbox

最適化関数でシンボリックな計算を使用する際は、以下のことを考慮しなければなりません。

1. 最適化の目的関数と制約関数は、x のように、ベクトルで定義する必要があります。ただし、シンボリック変数はベクトル値ではなく、スカラーまたは複素数値です。このため、ベクトルとスカラー間の変換が必要です。

2. 最適化の勾配と、場合によりますがヘッセ行列は、目的関数または制約関数の本体内で計算されることになっています。これは、シンボリックな勾配またはヘッセ行列を、目的関数または制約関数ファイル、あるいは関数ハンドル内の適切な位置に配置する必要があることを意味します。

3. 勾配とヘッセ行列のシンボリックな計算には、時間がかかることがあります。このため、ソルバーの実行中に呼び出すために、この計算を一度だけ実行し、matlabFunction を使用してコードを生成する必要があります。

4. 関数 subs を使用したシンボリック式の評価には、時間がかかります。matlabFunction を使用した方が、はるかに効率的です。

5. matlabFunction は、入力ベクトルの方向に依存するコードを生成します。fmincon は列ベクトルを含む目的関数を呼び出すため、シンボリック変数の列ベクトルを含む matlabFunction を呼び出す際は注意しなければなりません。

第 1 の例: ヘッセ行列を使用した制約なし最小化

最小化する目的関数は以下のとおりです。

$$f(x_1, x_2) = /log/left (1 + 3 /left ( x_2 - (x_1^3 - x_1)
/right )^2 + (x_1 - 4/3)^2 /right ).$$

この関数は正であり、x1 = 4/3, x2 =(4/3)^3 - 4/3 = 1.0370... において固有の最小値ゼロに到達します。

独立変数を x1 および x2 として記述します。その理由は、この形式で、独立変数をシンボリック変数として使用できるからです。ベクトル x の要素としては、x(1)x(2) と記述されます。この関数は、次のプロットで示したように、曲がりくねった谷をもっています。

syms x1 x2 real
x = [x1;x2]; % column vector of symbolic variables
f = log(1 + 3*(x2 - (x1^3 - x1))^2 + (x1 - 4/3)^2);

ezsurfc(f,[-2 2])
view(127,38)

f の勾配とヘッセ行列を計算します。

gradf = jacobian(f,x).'; % column gradf
hessf = jacobian(gradf,x);

fminunc ソルバーには、ベクトル x を渡し、GradObj オプションと Hessian オプションを [on] に設定し、次の 3 つの出力のリストを与えます。[f(x),gradf(x),hessf(x)]。

matlabFunction は、3 つの入力のリストから 3 つの出力のこのリストを正確に生成します。matlabFunction はさらに、vars オプションを使用してベクトル入力を受け取ります。

fh = matlabFunction(f,gradf,hessf,'vars',{x});

点 [-1,2] から始めて、最小化問題を解きます。

options = optimoptions('fminunc','GradObj','on','Hessian','on', ...
    'Display','final');
[xfinal fval exitflag output] = fminunc(fh,[-1;2],options)
Local minimum possible.

fminunc stopped because the final change in function value relative to 
its initial value is less than the default value of the function tolerance.




xfinal =

    1.3333
    1.0370


fval =

   7.6623e-12


exitflag =

     3


output = 

         iterations: 14
          funcCount: 15
       cgiterations: 11
      firstorderopt: 3.4391e-05
          algorithm: 'large-scale: trust-region Newton'
            message: 'Local minimum possible.

fminunc stopped because the...'
    constrviolation: []

これを、勾配またはヘッセ行列の情報を使用しない場合の反復の回数と比較します。比較に必要な中規模のアルゴリズムは、Algorithm オプションを 'quasi-newton' に設定することで得られます。

options = optimoptions('fminunc','Display','final','Algorithm','quasi-newton');
fh2 = matlabFunction(f,'vars',{x});
% fh2 = objective with no gradient or Hessian
[xfinal fval exitflag output2] = fminunc(fh2,[-1;2],options)
Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the function tolerance.




xfinal =

    1.3333
    1.0371


fval =

   2.1985e-11


exitflag =

     1


output2 = 

       iterations: 18
        funcCount: 81
         stepsize: 1
    firstorderopt: 2.4587e-06
        algorithm: 'medium-scale: Quasi-Newton line search'
          message: 'Local minimum found.

Optimization completed because t...'

反復回数は、勾配とヘッセ行列を使用した場合よりも少なく、関数評価の回数が激減します。

sprintf(['There were %d iterations using gradient' ...
    ' and Hessian, but %d without them.'], ...
    output.iterations,output2.iterations)
sprintf(['There were %d function evaluations using gradient' ...
    ' and Hessian, but %d without them.'], ...
    output.funcCount,output2.funcCount)
ans =

There were 14 iterations using gradient and Hessian, but 18 without them.


ans =

There were 15 function evaluations using gradient and Hessian, but 81 without them.

第 2 の例: fmincon の内点法アルゴリズムを使用した制約付き最小化

第 1 の例と同じ目的関数と開始値を検討しますが、第 2 の例には 2 つの非線形制約があります。

$$5/sinh(x_2/5) /ge x_1^4$$

$$5/tanh(x_1/5) /ge x_2^2 - 1.$$

この制約により、最適化が大域的最小点 [1.333,1.037] から離されます。この 2 つの制約を視覚化します。

[X,Y] = meshgrid(-2:.01:3);
Z = (5*sinh(Y./5) >= X.^4);
% Z=1 where the first constraint is satisfied, Z=0 otherwise
Z = Z+ 2*(5*tanh(X./5) >= Y.^2 - 1);
% Z=2 where the second is satisfied, Z=3 where both are
surf(X,Y,Z,'LineStyle','none');
set(gcf,'Color','w') % white background
view(0,90)
hold on
plot3(.4396, .0373, 4,'o','MarkerEdgeColor','r','MarkerSize',8);
% best point
xlabel('x');ylabel('y');
hold off

最適点の周囲に小さい赤色の円をプロットしました。

実行可能領域全体にわたって目的関数のプロットがあります。これは両方の制約を満たす領域であり、最適点の周囲の小さい赤色の円と共に濃い赤色で示されています。

W = log(1 + 3*(Y - (X.^3 - X)).^2 + (X - 4/3).^2);
% W = the objective function
W(Z < 3) = nan; % plot only where the constraints are satisfied
surf(X,Y,W,'LineStyle','none');
view(68,20)
hold on
plot3(.4396, .0373, .8152,'o','MarkerEdgeColor','r', ...
    'MarkerSize',8); % best point
xlabel('x');ylabel('y');zlabel('z');
hold off

非線形制約は、c(x) <= 0 という形式で記述しなければなりません。matlabFunction を使用して、すべてのシンボリック制約とその導関数を計算し、関数ハンドルに配置します。

これらの制約の勾配は、列ベクトルであることが望ましいです。これらは、各列が 1 つの制約関数の勾配を表すような行列として目的関数に配置しなければなりません。これは jacobian によって生成される形式の転置行列であるため、以下では転置行列を求めます。

非線形制約を関数ハンドルに配置します。fmincon は、非線形制約と勾配が [c ceq gradc gradceq] という順序で出力されると想定します。非線形等式制約が存在しないため、ceqgradceq として [ ] を出力します。

c1 = x1^4 - 5*sinh(x2/5);
c2 = x2^2 - 5*tanh(x1/5) - 1;
c = [c1 c2];
gradc = jacobian(c,x).'; % transpose to put in correct form
constraint = matlabFunction(c,[],gradc,[],'vars',{x});

内点法アルゴリズムでは、ヘッセ行列関数を、目的関数の一部ではなく別個の関数として記述しなければなりません。これは、非線形制約関数でこれらの制約をヘッセ行列に含める必要があるからです。内点法アルゴリズムのヘッセ行列は、ラグランジアンのヘッセ行列です。詳細は、ユーザー ガイドを参照してください。

ヘッセ行列関数には、2 つの入力引数を指定します。正のベクトル x とラグランジュ乗数構造体 lambda です。非線形制約に使用する lambda 構造体の部分は、lambda.ineqnonlinlambda.eqnonlin です。現在の制約の場合は、非線形等式が存在しないため、lambda.ineqnonlin(1)lambda.ineqnonlin(2) という乗数を使用します。

第 1 の例では、目的関数のヘッセ行列を計算しました。ここでは、2 つの制約関数のヘッセ行列を計算し、matlabFunction で関数ハンドル バージョンを作成します。

hessc1 = jacobian(gradc(:,1),x); % constraint = first c column
hessc2 = jacobian(gradc(:,2),x);

hessfh = matlabFunction(hessf,'vars',{x});
hessc1h = matlabFunction(hessc1,'vars',{x});
hessc2h = matlabFunction(hessc2,'vars',{x});

最終的なヘッセ行列を作成するために、3 つのヘッセ行列を統合し、適切なラグランジュ乗数を制約関数に追加します。

myhess = @(x,lambda)(hessfh(x) + ...
    lambda.ineqnonlin(1)*hessc1h(x) + ...
    lambda.ineqnonlin(2)*hessc2h(x));

内点法アルゴリズム、勾配、およびヘッセ行列を使用するためのオプションを設定し、目的関数が目的と勾配の両方を返すようにし、ソルバーを実行します。

options = optimoptions('fmincon','Algorithm','interior-point','GradObj',...
    'on','GradConstr','on','Hessian','user-supplied',...
    'HessFcn',myhess,'Display','final');
% fh2 = objective without Hessian
fh2 = matlabFunction(f,gradf,'vars',{x});
[xfinal fval exitflag output] = fmincon(fh2,[-1;2],...
    [],[],[],[],[],[],constraint,options)
Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.




xfinal =

    0.4396
    0.0373


fval =

    0.8152


exitflag =

     1


output = 

         iterations: 10
          funcCount: 13
    constrviolation: 0
           stepsize: 1.9160e-06
          algorithm: 'interior-point'
      firstorderopt: 1.9217e-08
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.
...'

この例でも、勾配とヘッセ行列が与えられている場合の方がそうでない場合よりも、ソルバーによって反復と関数評価の回数が少なくなっています。

options = optimoptions('fmincon','Algorithm','interior-point',...
    'Display','final');
% fh3 = objective without gradient or Hessian
fh3 = matlabFunction(f,'vars',{x});
% constraint without gradient:
constraint = matlabFunction(c,[],'vars',{x});
[xfinal fval exitflag output2] = fmincon(fh3,[-1;2],...
    [],[],[],[],[],[],constraint,options)

sprintf(['There were %d iterations using gradient' ...
    ' and Hessian, but %d without them.'],...
    output.iterations,output2.iterations)
sprintf(['There were %d function evaluations using gradient' ...
    ' and Hessian, but %d without them.'], ...
    output.funcCount,output2.funcCount)
Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.




xfinal =

    0.4396
    0.0373


fval =

    0.8152


exitflag =

     1


output2 = 

         iterations: 17
          funcCount: 54
    constrviolation: 0
           stepsize: 4.2417e-06
          algorithm: 'interior-point'
      firstorderopt: 3.8433e-07
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.
...'


ans =

There were 10 iterations using gradient and Hessian, but 17 without them.


ans =

There were 13 function evaluations using gradient and Hessian, but 54 without them.

シンボリック変数のクリーンアップ

この例で使用したシンボリック変数は、実数であると想定されていました。この想定をシンボリック エンジン ワークスペースから取り除くには、これらの変数を削除するだけでは不十分です。以下の構文を使用して、変数を消去するか

syms x1 x2 clear

以下のコマンドを使用して、シンボリック エンジンをリセットします。

% reset(symengine) % uncomment this line to reset the engine

シンボリック エンジンをリセットした後、すべてのシンボリック変数を MATLAB® ワークスペースから消去してください。

% clear % uncomment this line to clear the variables
この情報は役に立ちましたか?