Main Content

シミュレーションまたは常微分方程式の最適化

シミュレーションまたは常微分方程式の最適化とは

目的関数や非線形制約関数の値がシミュレーションまたは常微分方程式 (ODE) の数値解によってのみ得られる場合があります。このような最適化の問題には、潜在的な問題と解に説明するようにいくつかの共通の特徴と課題があります。

ODE を最適化する問題ベースの例については、最適化変数を使用して ODE のパラメーターを当てはめるを参照してください。ソルバーベースの例については、常微分方程式 (ODE) の適合を参照してください。

他の方法で発生する問題の多くを回避する方法については、離散化された最適軌跡 (問題ベース)を参照してください。この方法では、最適化プロセスで自動微分を使用できます。ただし、この方法は固定ステップに基づいており、低次の ODE 解のアルゴリズムを使用する可能性があるため、比較的精度が低くなる可能性があります。

Simulink® モデルを簡単に最適化するには、Simulink Design Optimization™ を使用してみてください。

潜在的な問題と解

有限差分法の問題

Optimization Toolbox™ ソルバーは、目的関数と制約関数の導関数を内部で使用します。既定では、これらの導関数は有限差分近似によって推定されます。これには以下の式を使用します。

F(x+δ)F(x)δ

または

F(x+δ)F(xδ)2δ.

これらの有限差分法による近似値は、次の理由で不正確な場合があります。

  • δ の値が大きいと、非線形性による有限差分への影響が大きくなる。

  • δ の値が小さいと、数値精度の制限により不正確になる。

特に シミュレーションや ODE の数値解では、次のことが言えます。

  • シミュレーションは通常、パラメーターの小さな変化には影響を受けません。つまり、摂動 δ が小さすぎると、シミュレーションにより導関数が間違って推定され、0 が返されることがあります。

  • シミュレーションと ODE の数値解は両方とも、関数評価に不正確さをもつ可能性があります。このような不正確さは有限差分法の近似において増幅されることがあります。

  • ODE の数値解では、機械精度よりもかなり大きい値がノイズとして導入されます。このノイズは有限差分法の近似において増幅されることがあります。

  • ODE ソルバーが可変ステップ サイズを使用する場合、F(x + δ) の評価の ODE ステップ数が F(x) の評価のステップ数と異なることがあります。この違いが原因となって戻り値に偽の差分が生じ、導関数の推測が誤ったものとなる可能性があります。

有限差分法に関するヒント

直接探索法を使用した有限差分法の回避.  Global Optimization Toolbox のライセンスをお持ちの場合は、ソルバーに patternsearch (Global Optimization Toolbox) を使用してください。patternsearch では勾配の推測は試みられないので、有限差分法の問題に説明するような制限の影響は受けません。

計算量の多い (時間のかかる) 関数評価に patternsearch を使用する場合は、次のように Cache オプションを使用します。

options = optimoptions('patternsearch','Cache','on');

patternsearch を使用できず、比較的低次元の制約なしの最小化問題を解く場合は、代わりに fminsearch を使用してください。fminsearch は有限差分を使用しません。ただし、fminsearch は高速なソルバーでも調整可能なソルバーでもありません。

大きな有限差分の設定.  有限差分法の問題の問題は、既定値よりも大きな有限差分のステップを設定することで回避できる場合があります。

  • MATLAB® R2011b 以降のリリースをお使いの場合は、有限差分のステップ サイズ オプションを既定の sqrt(eps) または eps^(1/3) よりも大きな値に設定します。たとえば、次のようにします。

    • R2011b – R2012b の場合:

      options = optimset('FinDiffRelStep',1e-3);
    • R2013a ~ R2015b およびソルバー名 'solvername' の場合:

      options = optimoptions('solvername','FinDiffRelStep',1e-3);
    • R2016a 以降およびソルバー名 'solvername' の場合:

      options = optimoptions('solvername','FiniteDifferenceStepSize',1e-3);

    さまざまな成分で異なるスケールを使用している場合は、有限差分のステップ サイズを成分のスケールに比例するベクトルに設定します。

  • MATLAB R2011a 以前のリリースをご使用の場合は、DiffMinChange オプションを既定の 1e-8 よりも大きな値に設定します。また、DiffMaxChange オプションを設定することもできます。たとえば、次のようにします。

    options = optimset('DiffMinChange',1e-3,'DiffMaxChange',1);

メモ

これらの有限差分のサイズをどのように設定するかを判断するのは困難です。

中心有限差分を設定することもできます。

options = optimoptions('solvername','FiniteDifferenceType','central');

勾配評価関数の使用.  有限差分推定の問題を回避するために、目的と非線形制約に近似勾配関数を与えることができます。optimoptions を使用して SpecifyObjectiveGradient オプションを true に設定し、該当する場合は SpecifyConstraintGradient オプションも true に設定することを忘れないでください。

  • 一部の ODE では、ODE を解くのと同時に勾配を数値的に評価できます。たとえば、目的関数 z(t,x) の微分方程式が次のようであるとします。

    ddtz(t,x)=G(z,t,x),

    ここで、x は最小化を行うパラメーターのベクトルです。x はスカラーとします。この場合、その導関数 y の微分方程式

    y(t,x)=ddxz(t,x)

    は以下になります。

    ddty(t,x)=G(z,t,x)zy(t,x)+G(z,t,x)x,

    ここで z(t,x) は目的関数 ODE の解です。y(t,x) は、z(t,x) と同じ連立微分方程式系で解くことができます。この解により、有限差分法に依らず近似された導関数が与えられます。非スカラーの x に対しては、成分ごとに 1 つの ODE を解きます。

    このメソッドの理論上および計算上の側面については、Leis と Kramer[2]を参照してください。このメソッドと有限差分メソッドによる計算上の経験については、資料 (Raue 他著) の図 7 を参照してください。[3].

  • 一部のシミュレーションでは、シミュレーション内で導関数を推定できます。たとえば、Reiman と Weiss の[4]に説明されている尤度比の手法や、Heidelberger、Cao、Zazanis および Suri の[1]で解析されている無限小摂動解析の手法は、目的関数や制約関数を評価するのと同じシミュレーションで導関数を推定します。

ODE 許容誤差を小さく設定.  odeset を使用して、AbsTol または RelTol の ODE ソルバー許容誤差を既定値よりも小さい値に設定できます。ただし、許容誤差を小さくしすぎると、解を求めるのに時間がかかったり、収束が失敗したりするなどの問題が発生する場合があります。RelTol では許容誤差は必ず 1e-9 以上にしてください。AbsTol の各成分の下限は、問題のスケールによって異なるため、特にアドバイスはありません。

確率的な関数での問題

シミュレーションで乱数を使用する場合、目的関数または制約関数を 2 回評価すると異なる結果が得られることがあります。これは、関数推定と有限差分推定の両方に影響を与えます。有限差分の値は、異なる評価点 x および x + δ による変動ではなく、無作為性による変動によって決まる可能性があります。

確率的な関数に関するヒント

制御ストリームからの乱数をシミュレーションで使用する場合は、乱数ストリームをリセットしてから目的関数または制約関数を評価します。これによって結果のばらつきを抑えることができます。たとえば、目的関数では次のようにします。

function f = mysimulation(x)
rng default % or any other resetting method
...
end

詳細については、再現可能な乱数の生成を参照してください。

目的と制約の共通の計算

多くの場合、シミュレーションでは目的関数と制約関数の両方が評価されますが、これは 1 回のシミュレーション実行中に行われます。つまり、時間のかかる計算を使用する場合でも、目的関数と非線形制約関数の両方で同じ計算を使用します。一方、fmincon などのソルバーでは、目的関数と非線形制約関数が別々に評価されます。この場合、ソルバーによって時間のかかる計算が 2 回呼び出されるため、効率性が大幅に低下することがあります。この問題を回避するには、同じ関数における目的と非線形制約の手法を使用するか、問題ベースのアプローチを使用する場合は共通の関数を持つ目的関数と制約の逐次評価または並列評価、問題ベースの手法を使用します。

目的関数または制約関数の評価の失敗

一部のパラメーター値では、シミュレーションまたは ODE が失敗することがあります。

評価の失敗に関するヒント

適切な範囲の設定.  パラメーター空間のあらゆる制限が既知とは限りませんが、すべてのパラメーターに適切な範囲 (上限と下限の両方) を設定するようにします。これによって最適化にかかる時間を短縮でき、問題のあるパラメーター値をソルバーで使わないようにできます。

範囲制約に従うソルバーの使用.  反復は制約に違反する可能性ありに説明するように、一部のアルゴリズムでは中間の反復で範囲制約に対する違反が発生してもかまいません。一方、シミュレーションと ODE の最適化では、範囲制約に常に従うアルゴリズムを使用します。詳細については、範囲制約を満たすアルゴリズムを参照してください。

NaN を返す.  シミュレーションまたは ODE ソルバーが点 x で目的関数または非線形制約関数を正常に評価できなかった場合は、関数で NaN を返すようにします。Optimization Toolbox と Global Optimization Toolbox のほとんどのソルバーは、NaN 値が発生した場合に別の反復ステップを試みるロバスト性をもっています。これらのロバストなソルバーには次のものがあります。

  • fmincon interior-pointsqp、および trust-region-reflective アルゴリズム

  • fminunc

  • lsqcurvefit

  • lsqnonlin

  • patternsearch

失敗に終わった点や実行不可能な点などの問題が発生するような点で、任意の大きな目的関数値を返すように試みる人もいますが、ソルバーでは戻り値が任意の値であることを認識できないため、ソルバーに混乱をきたす場合があります。ただし、NaN を返すと、ソルバーは別の点での評価を試みます。

参考文献

[1] Heidelberger, P., X.-R. Cao, M. A. Zazanis, and R. Suri. Convergence properties of infinitesimal perturbation analysis estimates. Management Science 34, No. 11, pp. 1281–1302, 1988.

[2] Leis, J. R. and Kramer, M.A. The Simultaneous Solution and Sensitivity Analysis of Systems Described by Ordinary Differential Equations. ACM Trans. Mathematical Software, Vol. 14, No. 1, pp. 45–60, 1988.

[3] Raue, A. et al. Lessons Learned from Quantitative Dynamical Modeling in Systems Biology. Available at https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0074335, 2013.

[4] Reiman, M. I. and A. Weiss. Sensitivity analysis via likelihood ratios. Proc. 18th Winter Simulation Conference, ACM, New York, pp. 285–289, 1986.