Sobol感度指標 — 分散分解に基づくグローバル感度解析

理論と物理

Sobol指標とは

Sobol感度指標は、1993年にロシアの数学者 I.M. Sobol が提案した分散分解（ANOVA分解）に基づくグローバル感度解析手法である。局所的な偏微分に基づく感度解析とは異なり、パラメータ空間全域にわたるモデル応答のばらつきを評価する。

ANOVA分解（分散分解）

$k$ 個の独立な入力変数 $X_1, X_2, \ldots, X_k$ を持つモデル出力 $Y = f(X_1, \ldots, X_k)$ について、ANOVA分解（高次元モデル表現: HDMR）は以下のように書ける：

高次交互作用指標

2次の交互作用指標は以下のように定義される：

ここで $S_{ij}^{\text{closed}} = \frac{V_{X_i, X_j}[E_{X_{\sim ij}}(Y|X_i, X_j)]}{V[Y]}$ は $X_i$ と $X_j$ の閉じた2次指標である。理論上は任意の次数まで計算可能だが、計算コストが指数的に増大するため、実務では $S_i$ と $S_{Ti}$ の2つで十分なことが多い。

数値解法と実装

Saltelli推定量

Saltelli推定量の手順：

Step 1：$N$ 行 $k$ 列のサンプル行列 $\mathbf{A}$ と $\mathbf{B}$ を独立に生成（準モンテカルロ列推奨）。

Step 2：各 $i = 1, \ldots, k$ に対して混合行列 $\mathbf{A}_B^{(i)}$ を構成。$j$ 列目は：

Step 3：モデル $f$ を $\mathbf{A}$, $\mathbf{B}$, $\mathbf{A}_B^{(1)}, \ldots, \mathbf{A}_B^{(k)}$ の全行に対して評価。合計 $N(k+2)$ 回のモデル評価が必要。

Step 4：推定量を計算。

ここで $V[Y]$ は $\mathbf{A}$ と $\mathbf{B}$ の全サンプルから推定する。

準モンテカルロサンプリング

サロゲートモデルとの併用

• 多項式カオス展開（PCE）：展開係数から解析的にSobol指標を直接計算できる。最も効率的
• Kriging（ガウス過程回帰）：予測の不確かさも定量化できるのが強み
• ニューラルネットワーク：高次元・強非線形に向くが、学習データ量が必要

Bootstrap信頼区間

Sobol指標の推定値には統計的な不確かさが伴う。Bootstrap法で信頼区間を評価するのが標準的だ：

1. 元の $N$ サンプルから重複を許して $N$ 個を再サンプリング
2. 再サンプリングデータでSobol指標を再計算
3. これを $B$ 回（通常1000回以上）繰り返し
4. 得られた $B$ 個の推定値からパーセンタイル信頼区間を構成

信頼区間が広い場合はサンプル数 $N$ を増やす必要がある。SALibでは sobol.analyze() の calc_second_order オプションで自動的にBootstrap信頼区間が計算される。

実践ガイド

Sobol解析のワークフロー

1. 問題定義：出力QoI（Quantity of Interest）と入力パラメータ $X_1, \ldots, X_k$ を明確にする。各入力の確率分布（一様分布・正規分布等）を決める
2. スクリーニング（任意）：パラメータ数が多い場合（20個以上）、まずMorris法で重要でないパラメータを除外する
3. サンプル行列生成：Saltelli法で $\mathbf{A}, \mathbf{B}, \mathbf{A}_B^{(i)}$ を生成。SALibの saltelli.sample() を使うのが最も楽
4. モデル評価：生成された全サンプル点でCAEモデルを実行。並列計算を推奨
5. Sobol指標の計算：評価結果を集めて推定量を計算。SALibの sobol.analyze() を使う
6. 結果の解釈と意思決定：$S_i, S_{Ti}$ からパラメータの重要度をランキングし、設計の最適化や不確かさ低減に反映する

サンプルサイズの決め方

結果の解釈とパラメータ固定

• $S_{Ti} \approx 0$（例えば $< 0.01$）→ そのパラメータは公称値で固定してよい。次元削減
• $S_i$ が大きい → そのパラメータのばらつきを制御（公差を厳しくする等）すれば出力のばらつきを効果的に低減できる
• $S_{Ti} - S_i$ が大きい → 他のパラメータとの交互作用が強い。単独で制御しても効果が薄い場合がある

• 板厚: $S_1 = 0.40,\ S_{T1} = 0.48$ — 最も支配的
• 降伏応力: $S_2 = 0.25,\ S_{T2} = 0.35$ — 2番目に重要、交互作用もあり
• 摩擦係数: $S_3 = 0.12,\ S_{T3} = 0.18$ — 中程度の寄与
• プレス速度: $S_4 = 0.02,\ S_{T4} = 0.03$ — ほぼ影響なし → 固定可
• ブランクホルダー力: $S_5 = 0.08,\ S_{T5} = 0.10$ — 小さいが無視はできない

この結果から、「板厚と降伏応力の管理が最優先」「プレス速度は気にしなくてよい」という明確な設計指針が出る。

Python SALib実装例

SALib（Sensitivity Analysis Library）を使った最小限のコード例：

import numpy as np
from SALib.sample import saltelli
from SALib.analyze import sobol

# 問題定義
problem = {
    'num_vars': 3,
    'names': ['thickness', 'yield_stress', 'friction'],
    'bounds': [[0.8, 1.2],       # mm
               [200, 350],       # MPa
               [0.05, 0.20]]     # 摩擦係数
}

# Saltelliサンプリング（Sobol列使用）
param_values = saltelli.sample(problem, 2048, calc_second_order=True)
# → shape: (2048*(3+2), 3) = (10240, 3)

# モデル評価（ここを実際のCAEソルバーに置き換える）
Y = np.zeros(param_values.shape[0])
for i, X in enumerate(param_values):
    Y[i] = your_cae_model(X[0], X[1], X[2])

# Sobol指標の計算
Si = sobol.analyze(problem, Y, calc_second_order=True,
                   conf_level=0.95, print_to_console=True)

# 結果
print("S1:", Si['S1'])       # 第1次指標
print("ST:", Si['ST'])       # 全次指標
print("S2:", Si['S2'])       # 2次交互作用指標
print("S1_conf:", Si['S1_conf'])  # 95%信頼区間

ソフトウェア比較

主要ツール比較

選定の指針

• 初めてSobol解析をやる → SALib一択。pip installして10行で結果が出る
• PCEベースで効率よくやりたい → UQLab（MATLAB）か OpenTURNS（Python）
• 大規模HPCジョブと連携 → Dakota。SLURM/PBS統合がしっかりしている
• GUIで手軽に、Ansys既存ユーザー → optiSLang。Workbenchから直接パラメトリックスタディ
• 相関入力への対応が必要 → OpenTURNSが最も柔軟

先端技術

群Sobol指標（Group Sobol Indices）

数学的には、$\mathbf{X}_u = \{X_i : i \in u\}$ としてグループ $u$ のSobol指標を $S_u = V[E(Y|\mathbf{X}_u)] / V[Y]$ と定義する。計算はSaltelli法の拡張で可能。SALibも groups パラメータで対応している。

時間依存Sobol指標

過渡解析の場合、各時刻 $t$ でのSobol指標 $S_i(t), S_{Ti}(t)$ を計算する。時刻によって支配的なパラメータが変わることがあり、例えば衝突解析では：

• 初期（$t < 5$ ms）：衝突速度 $S_1(t) \approx 0.8$ が支配的
• 中期（$5 < t < 20$ ms）：材料パラメータが $S_2(t) \approx 0.5$ に増大
• 後期（$t > 20$ ms）：摩擦や境界条件の影響が増す

時間依存解析では、各時刻を独立にSobol解析するか、主成分分析（PCA）で時系列データを縮約してからSobol解析を行う方法がある。

Shapley効果との比較

Shapley効果 $Sh_i$ は、変数 $X_i$ のすべての可能な結託（サブセット）に対する限界的な分散寄与の平均として定義される。$\sum_{i=1}^k Sh_i = V[Y]$ が常に成立し、「公平な分配」を実現する。ただし計算コストは $O(2^k)$ のサブセット評価が必要で、パラメータ数が増えると爆発する。

トラブルシューティング

指標が負になる

1. サンプル数 $N$ が不足：推定量はモンテカルロ推定なので、$N$ が小さいと統計的な揺らぎで負値が出る。$N$ を倍にして再計算してみよう
2. そのパラメータの真の $S_i$ がほぼゼロ：寄与がゼロに近い変数は推定量のノイズに埋もれて負になりやすい。信頼区間を見て、ゼロをまたいでいるなら「影響なし」と判断して問題ない

$\sum S_i$ が1を超える

• $N$ を増やす（最低でも2048以上）
• Sobol列（QMC）を使っているか確認。擬似乱数だと収束が遅い
• モデル出力に外れ値・不連続がないか確認。シミュレーションが途中で破綻しているケースがある

$N$ を増やしても収束しない

• モデルの出力分散が非常に大きい（裾の重い分布）→ 出力のロバスト化（対数変換等）を検討
• 入力に相関があるのに独立を仮定している → 相関を考慮した手法（Kucherenko法、Shapley効果）に切り替え
• モデルが非常に強い非線形性を持つ → サロゲートモデルで平滑化してからSobol解析を行う
• モデルにノイズがある（確率的シミュレーション）→ 同一入力で複数回実行して出力を平均する

まずは「$N$ を倍にしたときに推定値の変化が10%以内に収まるか」を収束判定の目安にするといい。