MPPI（Model Predictive Path Integral）の数理：クロスエントロピー法との統一的理解

MPPIとは

クロスエントロピー法（CEM） では、重点サンプリングにおけるサンプリング分布のパラメータを反復的に最適化する手法を紹介しました。

Model Predictive Path Integral（MPPI） は、確率最適制御に基づくサンプリングベースのモデル予測制御（MPC）アルゴリズムです。CEMとMPPIは、一見異なるアプローチに見えますが、実は重点サンプリングと変分推論という共通の数理的枠組みで統一的に理解できます。

• CEM: エリートサンプルの「ハード」な選択（上位 \(P\) %）
• MPPI: コストに基づく指数関数的な「ソフト」な重み付け

参考文献：Williams, G., et al. (2017). “Information Theoretic MPC for Model-Based Reinforcement Learning.” ICRA 2017.

問題設定

離散時間の確率的力学系を考えます：

ここで \(\mathbf{x}_t\) は状態、\(\mathbf{u}_t\) は制御入力、\(\boldsymbol{\epsilon}_t\) はシステムノイズです。

制御入力列 \(\mathbf{U} = (\mathbf{u}_0, \mathbf{u}_1, \ldots, \mathbf{u}_{T-1})\) に対するコスト関数を：

とします。\(\phi\) は終端コスト、\(q\) はステージコストです。

MPPIの導出

なぜ最適な制御分布がコストの指数関数 \(\exp(-J/\lambda)\) という特別な形をとるのでしょうか。これを天下り的に受け入れるのではなく、確率最適制御の基礎方程式から導出します。

経路積分制御理論の背景：HJB方程式の線形化（Kappenの理論）

MPPIの"Path Integral"（経路積分）という名前は、Kappen (2005) による**線形可解な確率制御（linearly solvable stochastic control）**の理論に由来します。連続時間・制御アフィンな確率微分方程式

を考え、コスト汎関数を終端コストと状態コスト \(q(\mathbf{x})\) ・2次形式の制御コスト \(\frac12 \mathbf{u}^T R \mathbf{u}\) の和とします：

価値関数 \(V(\mathbf{x}, t) = \min_{\mathbf{u}(\cdot)} C(\mathbf{x}, t, \mathbf{u}(\cdot))\) が満たすハミルトン・ヤコビ・ベルマン（HJB）方程式は、伊藤の公式から次のようになります：

右辺は \(\mathbf{u}\) について凸な2次関数なので、微分してゼロと置くことで解析的に最小化できます：

これを (P3) に代入すると：

このPDEは \(\nabla V\) の2次項を含むため非線形です。これが一般の確率最適制御問題が解析的に解けない理由です。

Kappenはここで、ノイズ共分散と制御コストの重みが次の関係を満たすという仮定を置きました（「システムが受けるノイズの大きさは、制御の安さに比例する」）：

この仮定のもとで「望ましさ (desirability)」関数を

と定義します。\(V = -\lambda \log \psi\) なので、連鎖律により

が成り立ちます。これらを (P5) に代入すると、\(\nabla\psi(\nabla\psi)^T/\psi^2\) に比例する項が2箇所（\((\nabla V)^T\mathbf{G}R^{-1}\mathbf{G}^T\nabla V\) の展開と \(\mathrm{tr}(\mathbf{G}\boldsymbol{\Sigma}\mathbf{G}^T\nabla^2 V)\) の展開）に現れ、その係数は \(\mathbf{G}\boldsymbol{\Sigma}\mathbf{G}^T - \lambda\mathbf{G}R^{-1}\mathbf{G}^T\) という形にまとまります。これは仮定(P6)より 恒等的にゼロ です。つまり非線形項は打ち消し合うのではなく、Kappenの仮定によって係数そのものが消えます。残った項を整理すると：

境界条件 \(\psi(\mathbf{x}, T) = \exp(-\phi(\mathbf{x})/\lambda)\) のもとで、式(P9)は \(\psi\) について線形です。非線形だったHJB方程式(P5)が、指数変換(P7)とノイズ・コストの整合条件(P6)によって線形の2階偏微分方程式に変換されました。これが「線形可解な確率制御」と呼ばれる理由です。

式(P9)の右辺第1・2項は、まさに無制御（\(\mathbf{u}=\mathbf{0}\) 、以下「受動的」と呼ぶ）な拡散過程 \(d\mathbf{x}=\mathbf{f}\,dt+\mathbf{G}\,d\boldsymbol{\omega}\) の生成作用素 \(\mathcal{L}\psi = \mathbf{f}^T\nabla\psi + \frac12\mathrm{tr}(\mathbf{G}\boldsymbol{\Sigma}\mathbf{G}^T\nabla^2\psi)\) です。したがって(P9)は \(\partial_t\psi + \mathcal{L}\psi - \frac{q(\mathbf{x})}{\lambda}\psi = 0\) という標準形の線形PDEであり、Feynman-Kacの定理により、受動的な拡散過程の経路に関する期待値として表現できます：

この期待値は無数のランダムな経路上の積分（=経路積分）であり、これがPath Integral Controlという名前の由来です。最適制御則は \(\psi\) の勾配から直接復元できます：

MPPIとの接続：MPPIは、(P10)の期待値をPDEを解析的に解く代わりにモンテカルロ近似する手法です。受動的な力学（ノイズのみを乗せた軌道）から \(N\) 本のロールアウトをサンプリングし、各軌道のコスト \(J(\mathbf{V}^{(i)})\) を評価して \(\frac1N\sum_i \exp(-J^{(i)}/\lambda)\) を計算することは、まさに(P10)の \(\psi\) の経験的推定にほかなりません。次節で導出する重み付き平均則は、この \(\psi\) の推定値から最適制御(P11)を復元する近似手続きに対応します。すなわち、離散時間・有限ホライズンのMPPIアルゴリズムは、Kappenの連続時間理論を後退ホライズンで数値的に解くレシピとして位置づけられます。

最適制御分布

確率最適制御の枠組みでは、制御列を確率変数とみなし、最適制御分布を以下のように定義します：

ここで \(\mathbf{V}\) はノイズ付き制御列、\(\lambda > 0\) は温度パラメータ、\(p(\mathbf{V})\) は事前分布（現在のサンプリング分布）です。前節の記法では、離散時間近似のもとで \(J(\mathbf{V})/\lambda\) が連続時間の \(-\log\psi\) に対応し、式(3)はまさに望ましさ関数 \(\psi\) による経路の再重み付けです。

温度 \(\lambda\) の役割：

• \(\lambda \to 0\) : 最小コストの制御列のみを選択（ハードな選択）
• \(\lambda \to \infty\) : すべての制御列を均等に重み付け

変分推論としての導出：モーメント射影（M-projection）

式(3)の \(p^*(\mathbf{V})\) から直接サンプリングすることはできません（それができるなら最適制御問題はすでに解けています）。そこで、現在のサンプリング方策と同じガウス族 \(q(\mathbf{V}; \mathbf{U}) = \mathcal{N}(\mathbf{U}, \boldsymbol{\Sigma})\) の中から、\(p^*\) に最も近い分布を探すという近似問題に置き換えます。「近さ」の尺度としてKLダイバージェンスを用い、次のモーメント射影 (M-projection) を解きます：

（モードを追いかける逆向きのKL、\(D_{\mathrm{KL}}(q \| p^*)\) を最小化するI-projectionとは方向が逆であることに注意してください。I-projectionは一般に閉形式を持ちませんが、M-projectionは指数型分布族に対して閉形式のモーメント整合則を持ちます。）

\(\boldsymbol{\Sigma}\) を固定して \(\mathbf{U}\) についてのみ展開すると、

となり、\(\mathbf{U}\) に関する勾配をゼロと置くと \(\boldsymbol{\Sigma}^{-1}\mathbb{E}_{p^*}[\mathbf{V}-\mathbf{U}] = \mathbf{0}\) 、すなわち

というモーメント一致が最適解であることが厳密に示せます（ガウス分布へのKL射影が常に一次モーメントを一致させることに帰着するという、指数型分布族の一般的性質の一例です）。式(3)の \(p^*(\mathbf{V}) \propto \exp(-J(\mathbf{V})/\lambda)p(\mathbf{V})\) を代入すると、

この積分は \(J\) が非線形な力学のロールアウトを含む汎関数であるため閉形式では計算できませんが、これは分布 \(p(\mathbf{V})\) （現在のサンプリング方策）のもとでの期待値そのものです。\(N\) 個のサンプル \(\mathbf{V}^{(i)} \sim p(\mathbf{V})\) を生成し、自己正規化重点サンプリング（self-normalized importance sampling）でこの期待値をモンテカルロ推定すると：

が得られます。これは次節の式(4)・式(5)そのものです。つまり MPPIの重み付き平均更新則は、「最適経路分布 \(p^*\) へのガウス族の中でのモーメント射影を、自己正規化重点サンプリングでモンテカルロ推定したもの」 として厳密に導出できます。「変分推論としてのMPC」という呼び方は、この特定の射影方向（M-projection）とガウス族という制約のもとで初めて意味を持つ主張です。

重み付き制御更新

\(N\) 個のサンプル \(\mathbf{V}^{(1)}, \ldots, \mathbf{V}^{(N)}\) を事前分布 \(p(\mathbf{V})\) から生成し、各サンプルのコスト \(J(\mathbf{V}^{(i)})\) を計算します。

重点サンプリングの重みは以下で計算されます（式(3d)と同一）：

最適制御入力は重み付き平均として求められます：

温度パラメータ \(\lambda\) の役割：数値実験

式(4)の重み \(w^{(i)}\) は理論的には次のように振る舞います：

• \(\lambda \to 0\) ：コスト差 \(J(\mathbf{V}^{(i)}) - J_{\min}\) がわずかでも指数関数的に増幅され、重みは最小コストのサンプル1つに退化的に集中する（ハードな選択に近づく）
• \(\lambda \to \infty\) ：すべてのサンプルがほぼ等しい重みを持ち、更新は事前分布からのランダムな平均に近づく（探索は広いが情報を活用できない）

この集中度合いを定量化する指標として、重点サンプリングの実効サンプルサイズ (Effective Sample Size, ESS) を用います：

ESSは重みが一様なら \(N\) に、単一サンプルへの退化的集中では \(1\) に近づきます。以下、本記事のPython実装（後述）をロールアウトについてNumPyでベクトル化した上で、\(N=500\) 固定・乱数シード42で \(\lambda \in \{0.05, 0.2, 1.0, 5.0, 20.0\}\) を比較し、各ステップのESSと閉ループでの累積コストを計測しました。

import numpy as np

def dynamics_batch(x, u, dt=0.1):
    """バッチ版力学モデル。x: (N, 4), u: (N, 2) -> (N, 4)"""
    x_next = np.empty_like(x)
    x_next[:, 0] = x[:, 0] + x[:, 2] * dt
    x_next[:, 1] = x[:, 1] + x[:, 3] * dt
    x_next[:, 2] = x[:, 2] + u[:, 0] * dt
    x_next[:, 3] = x[:, 3] + u[:, 1] * dt
    return x_next

def stage_cost_batch(x, u, goal):
    dist = np.sum((x[:, :2] - goal) ** 2, axis=1)
    ctrl = 0.01 * np.sum(u**2, axis=1)
    return dist + ctrl

def terminal_cost_batch(x, goal):
    return 10.0 * np.sum((x[:, :2] - goal) ** 2, axis=1)

class MPPIControllerVec:
    """実験用にロールアウトをベクトル化したMPPIコントローラ。
    数式上は元のMPPIControllerと同一で、Nサンプルをforループでなく
    NumPy配列演算として一括評価するだけの違い（実運用のMPPIも
    GPU上でこの形の並列化を行う）。"""
    def __init__(self, N, lam, T=20, n_ctrl=2, sigma=0.5, goal=np.array([5.0, 5.0])):
        self.N, self.lam, self.T, self.sigma, self.goal = N, lam, T, sigma, goal
        self.mu = np.zeros((T, n_ctrl))

    def compute_control(self, x0, rng):
        N, T = self.N, self.T
        noise = rng.standard_normal((N, T, self.mu.shape[1])) * self.sigma
        V = self.mu[np.newaxis, :, :] + noise

        x = np.tile(x0, (N, 1))
        costs = np.zeros(N)
        for t in range(T):
            costs += stage_cost_batch(x, V[:, t, :], self.goal)
            x = dynamics_batch(x, V[:, t, :])
        costs += terminal_cost_batch(x, self.goal)

        costs_shifted = costs - np.min(costs)
        weights = np.exp(-costs_shifted / self.lam)
        weights /= np.sum(weights)
        ess = 1.0 / np.sum(weights**2)  # 式(6)

        self.mu += np.sum(weights[:, np.newaxis, np.newaxis] * noise, axis=0)
        u_opt = self.mu[0].copy()
        self.mu = np.roll(self.mu, -1, axis=0)
        self.mu[-1] = 0.0
        return u_opt, ess


def run_episode(N, lam, n_steps=150, seed=42, goal=np.array([5.0, 5.0])):
    rng = np.random.default_rng(seed)
    x = np.array([0.0, 0.0, 0.0, 0.0])
    controller = MPPIControllerVec(N=N, lam=lam, goal=goal)
    ess_hist, total_cost = [], 0.0
    for _ in range(n_steps):
        u, ess = controller.compute_control(x, rng)
        ess_hist.append(ess)
        total_cost += stage_cost_batch(x[np.newaxis, :], u[np.newaxis, :], goal)[0]
        x = dynamics_batch(x[np.newaxis, :], u[np.newaxis, :])[0]
    return np.mean(ess_hist), total_cost


for lam in [0.05, 0.2, 1.0, 5.0, 20.0]:
    mean_ess, total_cost = run_episode(N=500, lam=lam, seed=42)
    print(f"lambda={lam:>6}: mean_ESS={mean_ess:6.1f} ({mean_ess/500*100:5.1f}% of N), "
          f"total_cost={total_cost:.2f}")

実行結果（乱数シード42、\(N=500\) 、150ステップ）：

理論通り、\(\lambda\) が大きいほどESSは単調に増加し（\(N\) の3.5%→86.8%）、重みがより多くのサンプルに分散していることが確認できます。一方、累積コストは \(\lambda\) が小さいほど良好（\(\lambda=0.05\) で709.41が最小、\(\lambda=20\) で878.48が最悪）という結果になりました。今回のようになめらかで単峰な2次コスト関数では、コストの低いサンプルを積極的に信頼する「貪欲な」重み付け（小さい \(\lambda\) ）が有利であることを示しています。逆に \(\lambda\) が大きすぎると、低コストサンプルが持つ情報が事前分布に近いほぼ一様な平均によって薄められ、性能が劣化します。

これは「探索の広さ」と「活用の精度」のトレードオフですが、その最適点はコスト関数の滑らかさに依存する点に注意が必要です。後述の「コスト関数設計の落とし穴」では、コスト関数側の性質（勾配情報の乏しさ）が同じ数式的メカニズム（重みの一様化）を引き起こし、\(\lambda\) を下げるだけでは解決できないケースを示します。

CEMとの比較

CEMとMPPIは共に、現在の分布から制御列をサンプリングし、コスト情報に基づいて分布を更新するという構造を持ちます。

CEM（ハードな選択）

CEMでは上位 \(P\) % のエリートサンプルを選択し、分布パラメータを更新します：

ここで \(\mathcal{E}\) はエリートサンプルの集合です。

MPPI（ソフトな重み付け）

MPPIでは全サンプルを指数関数的な重みで使用します：

統一的な理解：どちらもモーメント射影である

式(3a)-(3d)で導出した通り、MPPIの更新則(9)は「重み関数 \(w(\mathbf{V}) \propto \exp(-J(\mathbf{V})/\lambda)\) で再重み付けした分布へのモーメント射影」でした。CEMの更新則(7)-(8)も、実は全く同じ枠組みに収まります。すなわち、重み関数を

（\(J_{(P)}\) は下位 \(P\) % 分位点、\(\mathbb{1}[\cdot]\) は指示関数）というハードな0/1の重みに置き換えれば、式(3a)-(3d)と全く同じモーメント射影の議論がそのままCEMの更新則を再現します。つまり CEMとMPPIの違いは、「\(p^*\) への同じモーメント射影という問題を、どのような重み関数 \(w(\mathbf{V})\) で近似するか」という一点に集約されます：

CEMは \(w\) を指示関数（ハード閾値）、MPPIは \(w\) を指数関数（ソフト重み）にとった特殊ケースです。この統一的な視点から、両者の関係をより精密に述べられます。\(\lambda \to 0\) の極限では、式(4)の重みはコスト最小の1サンプルのみに退化的に集中します（\(N\) 個のサンプル中、最小コストのサンプルのみが重み1を持ち、他は0）。これは、CEMにおいてエリート数を \(|\mathcal{E}|=1\) まで絞った極限（分散更新なしの点推定）に対応する、退化的なケースです。一般の \(0 < P < 1\) に対するCEMの上位 \(P\) % 選択は、MPPIの温度 \(\lambda\) を0にする操作とは異なる重み関数（指示関数 vs 指数関数）を採用しているため、両者は「同じ枠組みの2つの異なるインスタンス」であって、片方がもう片方の厳密な極限というわけではない点に注意してください。

MPPIアルゴリズム

入力: 初期制御列 μ = (μ_0, ..., μ_{T-1}), サンプル数 N, 温度 λ, ノイズ分散 Σ
繰り返し:
  1. ノイズ付き制御列をN個サンプリング:
     ε_t^(i) ~ N(0, Σ),  V^(i) = μ + ε^(i)
  2. 各サンプルのコストを計算:
     J(V^(i)) = φ(x_T^(i)) + Σ_t q(x_t^(i), v_t^(i))
  3. 重みを計算:
     w^(i) = exp(-J(V^(i))/λ) / Σ_j exp(-J(V^(j))/λ)
  4. 制御入力を更新:
     μ_t ← μ_t + Σ_i w^(i) ε_t^(i)
  5. 最初の制御入力 μ_0 を適用し、制御列をシフト

Python実装

2次元のゴール到達問題でMPPIを実装します。

問題の定義

import numpy as np
import matplotlib.pyplot as plt

# ---- 問題設定 ----
n_state = 4     # 状態次元 [x, y, vx, vy]
n_ctrl = 2      # 制御次元 [ax, ay]
dt = 0.1        # 時間刻み
T = 20          # 予測ホライズン
N = 500         # サンプル数
lam = 1.0       # 温度パラメータ
sigma = 0.5     # 制御ノイズの標準偏差
goal = np.array([5.0, 5.0])  # ゴール位置

def dynamics(x, u):
    """線形力学モデル（質点の2次元運動）"""
    x_next = np.zeros(n_state)
    x_next[0] = x[0] + x[2] * dt  # px
    x_next[1] = x[1] + x[3] * dt  # py
    x_next[2] = x[2] + u[0] * dt  # vx
    x_next[3] = x[3] + u[1] * dt  # vy
    return x_next

def stage_cost(x, u):
    """ステージコスト"""
    pos = x[:2]
    dist = np.sum((pos - goal)**2)
    ctrl = 0.01 * np.sum(u**2)
    return dist + ctrl

def terminal_cost(x):
    """終端コスト"""
    pos = x[:2]
    return 10.0 * np.sum((pos - goal)**2)

MPPIコントローラ

class MPPIController:
    def __init__(self):
        # 制御列の初期化（ゼロ）
        self.mu = np.zeros((T, n_ctrl))

    def compute_control(self, x0):
        """MPPIによる最適制御入力の計算"""
        # ステップ1: ノイズ付き制御列のサンプリング
        noise = np.random.randn(N, T, n_ctrl) * sigma
        V = self.mu[np.newaxis, :, :] + noise  # (N, T, n_ctrl)

        # ステップ2: 各サンプルのコスト計算
        costs = np.zeros(N)
        for i in range(N):
            x = x0.copy()
            for t in range(T):
                costs[i] += stage_cost(x, V[i, t])
                x = dynamics(x, V[i, t])
            costs[i] += terminal_cost(x)

        # ステップ3: 重みの計算（式4）
        # 数値安定性のためにコストの最小値を引く
        costs_shifted = costs - np.min(costs)
        weights = np.exp(-costs_shifted / lam)
        weights /= np.sum(weights)  # 正規化

        # ステップ4: 制御入力の更新（式5）
        weighted_noise = np.sum(
            weights[:, np.newaxis, np.newaxis] * noise, axis=0
        )
        self.mu += weighted_noise

        # 最適制御入力（最初のタイムステップ）
        u_opt = self.mu[0].copy()

        # 制御列のシフト（次のタイムステップの準備）
        self.mu = np.roll(self.mu, -1, axis=0)
        self.mu[-1] = 0.0  # 最後の要素をゼロに

        return u_opt

シミュレーション

np.random.seed(42)

# 初期状態
x = np.array([0.0, 0.0, 0.0, 0.0])
controller = MPPIController()

# シミュレーション
n_steps = 100
trajectory = [x.copy()]
controls = []

for step in range(n_steps):
    u = controller.compute_control(x)
    controls.append(u.copy())
    x = dynamics(x, u)
    trajectory.append(x.copy())

    # ゴール到達判定
    if np.linalg.norm(x[:2] - goal) < 0.1:
        print(f"Goal reached at step {step + 1}")
        break

trajectory = np.array(trajectory)
controls = np.array(controls)

# ---- 結果のプロット ----
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 軌道
axes[0].plot(trajectory[:, 0], trajectory[:, 1], "b-o",
             markersize=3, label="MPPI trajectory")
axes[0].plot(*goal, "r*", markersize=15, label="Goal")
axes[0].plot(0, 0, "gs", markersize=10, label="Start")
axes[0].set_xlabel("x")
axes[0].set_ylabel("y")
axes[0].set_title("MPPI - Goal Reaching")
axes[0].legend()
axes[0].grid(True)
axes[0].set_aspect("equal")

# 制御入力
axes[1].plot(controls[:, 0], label="$a_x$")
axes[1].plot(controls[:, 1], label="$a_y$")
axes[1].set_xlabel("Step")
axes[1].set_ylabel("Control input")
axes[1].set_title("Control inputs over time")
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.savefig("mppi_result.png", dpi=150)
plt.show()

サンプル数 \(N\) の影響（数値実験）

式(3c)-(3d)の導出からわかる通り、重み付き平均は積分 \(\mathbb{E}_{p(\mathbf{V})}[\mathbf{V}\exp(-J(\mathbf{V})/\lambda)]\) のモンテカルロ推定です。モンテカルロ推定の誤差は一般に \(O(1/\sqrt{N})\) で減少するため、サンプル数 \(N\) を増やすほど制御性能のばらつきが小さくなり、平均的な性能も改善することが期待されます。これを確認するため、\(\lambda=1.0\) に固定し、\(N \in \{20, 50, 100, 500, 2000\}\) について乱数シードを10通り変えて閉ループシミュレーションを行いました。

import numpy as np

goal = np.array([5.0, 5.0])
dt, T, sigma = 0.1, 20, 0.5

def dynamics_batch(x, u):
    x_next = np.empty_like(x)
    x_next[:, 0] = x[:, 0] + x[:, 2] * dt
    x_next[:, 1] = x[:, 1] + x[:, 3] * dt
    x_next[:, 2] = x[:, 2] + u[:, 0] * dt
    x_next[:, 3] = x[:, 3] + u[:, 1] * dt
    return x_next

def stage_cost_batch(x, u):
    dist = np.sum((x[:, :2] - goal) ** 2, axis=1)
    ctrl = 0.01 * np.sum(u**2, axis=1)
    return dist + ctrl

def terminal_cost_batch(x):
    return 10.0 * np.sum((x[:, :2] - goal) ** 2, axis=1)

def run_episode(N, lam, seed, n_steps=150):
    rng = np.random.default_rng(seed)
    x = np.array([0.0, 0.0, 0.0, 0.0])
    mu = np.zeros((T, 2))
    total_cost = 0.0
    reached_step = None
    for step in range(n_steps):
        noise = rng.standard_normal((N, T, 2)) * sigma
        V = mu[np.newaxis, :, :] + noise
        xb = np.tile(x, (N, 1))
        costs = np.zeros(N)
        for t in range(T):
            costs += stage_cost_batch(xb, V[:, t, :])
            xb = dynamics_batch(xb, V[:, t, :])
        costs += terminal_cost_batch(xb)
        w = np.exp(-(costs - costs.min()) / lam)
        w /= w.sum()
        mu += np.sum(w[:, np.newaxis, np.newaxis] * noise, axis=0)
        u = mu[0].copy()
        total_cost += stage_cost_batch(x[np.newaxis, :], u[np.newaxis, :])[0]
        x = dynamics_batch(x[np.newaxis, :], u[np.newaxis, :])[0]
        mu = np.roll(mu, -1, axis=0)
        mu[-1] = 0.0
        if np.linalg.norm(x[:2] - goal) < 0.1 and reached_step is None:
            reached_step = step + 1
    final_dist = np.linalg.norm(x[:2] - goal)
    return final_dist, total_cost, (reached_step if reached_step is not None else n_steps)

for N in [20, 50, 100, 500, 2000]:
    dists, costs, reached = [], [], []
    for seed in range(10):
        d, c, r = run_episode(N=N, lam=1.0, seed=seed)
        dists.append(d)
        costs.append(c)
        reached.append(r)
    dists, costs, reached = np.array(dists), np.array(costs), np.array(reached)
    print(f"N={N:>5}: reached_step={reached.mean():.1f}+-{reached.std():.1f}, "
          f"final_dist={dists.mean():.4f}+-{dists.std():.4f}, "
          f"total_cost={costs.mean():.2f}+-{costs.std():.2f}")

実行結果（10シード平均±標準偏差、150ステップ、\(\lambda=1.0\) ）：

すべての条件で最終的にゴールへの到達（距離0.1未満）には成功しましたが（10/10）、\(N\) を20から2000に増やすことで、累積コストの平均は833.58から697.11へと約16%改善し、同時に標準偏差も24.91から8.58へと約3分の1に縮小しています。これはモンテカルロ推定の分散が \(N\) に反比例して減少するという理論的予測と整合的です。実用上は、計算予算（1制御周期あたりに評価できるロールアウト数）と要求される制御性能のばらつき許容度とのトレードオフとして \(N\) を選ぶことになります。

コスト関数設計の落とし穴：勾配情報が乏しい場合の発散

MPPIの重み(4)は、コストの相対的な差 \(J(\mathbf{V}^{(i)}) - J_{\min}\) に基づいて決まります。したがって、コスト関数がある領域でほぼ一定（勾配情報がほとんどない、あるいは局所的に平坦）だと、その領域にいる限りサンプル間のコスト差がほぼゼロになり、重みはほぼ一様になります。これは前節で見た「\(\lambda\) が大きすぎる場合」と数式的に全く同じ現象（重みの一様化）が、温度ではなくコスト関数側の性質によって引き起こされるケースです。

これを確認するため、通常の2次コスト（stage_cost_batch）に対し、ゴールから半径1.0より外側では一定値（勾配ゼロ）、内側でのみ2次関数になる「プラトー型」コストを用意し、比較しました：

import numpy as np

goal = np.array([5.0, 5.0])
dt, T, sigma, lam, N = 0.1, 20, 0.5, 1.0, 500

def dynamics_batch(x, u):
    x_next = np.empty_like(x)
    x_next[:, 0] = x[:, 0] + x[:, 2] * dt
    x_next[:, 1] = x[:, 1] + x[:, 3] * dt
    x_next[:, 2] = x[:, 2] + u[:, 0] * dt
    x_next[:, 3] = x[:, 3] + u[:, 1] * dt
    return x_next

def stage_cost_quadratic(x, u):
    return np.sum((x[:, :2] - goal) ** 2, axis=1) + 0.01 * np.sum(u**2, axis=1)

def terminal_cost_quadratic(x):
    return 10.0 * np.sum((x[:, :2] - goal) ** 2, axis=1)

def stage_cost_flat(x, u):
    """ゴール近傍(半径1.0)の外は一定値50（勾配ゼロ）、内側のみ2次関数。"""
    dist = np.linalg.norm(x[:, :2] - goal, axis=1)
    ctrl = 0.01 * np.sum(u**2, axis=1)
    return np.where(dist <= 1.0, 50.0 * dist**2, 50.0) + ctrl

def terminal_cost_flat(x):
    dist = np.linalg.norm(x[:, :2] - goal, axis=1)
    return np.where(dist <= 1.0, 500.0 * dist**2, 500.0)

def run_episode(stage_cost, terminal_cost, seed, n_steps=150):
    rng = np.random.default_rng(seed)
    x = np.array([0.0, 0.0, 0.0, 0.0])
    mu = np.zeros((T, 2))
    ess_hist = []
    for _ in range(n_steps):
        noise = rng.standard_normal((N, T, 2)) * sigma
        V = mu[np.newaxis, :, :] + noise
        xb = np.tile(x, (N, 1))
        costs = np.zeros(N)
        for t in range(T):
            costs += stage_cost(xb, V[:, t, :])
            xb = dynamics_batch(xb, V[:, t, :])
        costs += terminal_cost(xb)
        w = np.exp(-(costs - costs.min()) / lam)
        w /= w.sum()
        ess_hist.append(1.0 / np.sum(w**2))
        mu += np.sum(w[:, np.newaxis, np.newaxis] * noise, axis=0)
        u = mu[0].copy()
        x = dynamics_batch(x[np.newaxis, :], u[np.newaxis, :])[0]
        mu = np.roll(mu, -1, axis=0)
        mu[-1] = 0.0
    final_dist = np.linalg.norm(x[:2] - goal)
    return final_dist, np.mean(ess_hist)

for label, sc, tc in [("quadratic (baseline)", stage_cost_quadratic, terminal_cost_quadratic),
                       ("flat/plateau (adversarial)", stage_cost_flat, terminal_cost_flat)]:
    dists, ess_fracs = [], []
    for seed in range(10):
        d, ess = run_episode(sc, tc, seed=seed)
        dists.append(d)
        ess_fracs.append(ess / N * 100)
    dists = np.array(dists)
    n_success = int(np.sum(dists < 0.1))
    print(f"{label:>28}: success={n_success}/10, final_dist={dists.mean():.4f}+-{dists.std():.4f}, "
          f"mean_ESS={np.mean(ess_fracs):.2f}% of N")

実行結果（10シード、\(N=500\) 、\(\lambda=1.0\) 、150ステップ）：

結果は劇的です。通常の2次コストでは10/10のシードでゴールに到達しました（最終距離誤差0.0154）が、プラトー型コストでは10シード全てで到達に失敗し、150ステップ経過時点でゴールから平均6.98離れた位置で停滞しました。平均ESSを見ると、通常コストでは51.9%であるのに対し、プラトー型コストでは**99.95%**とほぼ \(N\) に等しく、重みがほぼ完全に一様化していることがわかります。これは前節の \(\lambda \to \infty\) の極限（式(6)のESSが\(N\) に近づく状況）と全く同じ数式的帰結です。

この現象のメカニズムは次の通りです。ロボットが半径1.0の外側にいる間、すべてのロールアウトのコストはほぼ同じ定数（プラトー値50付近、制御コストの微小な差のみ）になるため、\(J(\mathbf{V}^{(i)}) - J_{\min}\) が \(\lambda\) に対して常に小さく、重みはほぼ一様になります。式(3d)の重み付き平均はほぼ事前分布（ノイズの単純平均、期待値ゼロ）に一致し、制御入力の更新はほぼランダムウォークに退化します。したがってロボットはゴールに向かう有意な駆動力を得られず、プラトー領域を彷徨い続けます。この失敗は温度 \(\lambda\) をいくら下げても、勾配情報がそもそもコスト関数に存在しない限り解消しません。実務上の教訓として、MPPIを含むサンプリングベース最適化ではコスト関数の設計が性能を支配する最重要要素であり、目的地に向けた密な（できれば滑らかで単調な）勾配情報をコスト関数自体に埋め込む必要があります。

近年の研究動向（2023年以降）

MPPIは元来ロボティクス分野で発展してきましたが、2023年以降も応用・拡張の両面で活発に研究が進んでいます。

• 脚ロボットのリアルタイム全身制御：Alvarez-Padilla et al. (2024) は、接触モデルを含む脚ロボットの全身力学に対して、オフライン学習なしにMPPIをハードウェア上でリアルタイム実行する手法を報告しています。MPPIは勾配を必要とせず並列化が容易なため、接触が絡む非平滑な力学・コスト関数でも安定に動作する点が評価されています。
• 学習済みニューラルネットワークとの統合：自動運転分野では、Ryu et al. (2025) の IANN-MPPI が、周囲車両との相互作用を考慮したニューラルネットワークベースの事前分布・コスト設計をMPPIに組み込む手法を提案しています。本記事で見たように、コスト関数の勾配情報の乏しさがMPPIの性能を大きく左右するため、学習によってコスト関数や事前分布（サンプリング分布）自体をデータから獲得しようとする研究の流れは、この記事のエッジケースの実験結果とも整合的です。
• CEMとのハイブリッド化によるサンプル効率の改善：本記事のサンプル数実験で見た通り、MPPIの性能は \(N\) に依存してばらつきます。近年は、CEMのようなエリート選択と組み合わせて実効的なサンプル効率を高めるハイブリッド手法（脚ロボットの歩容最適化などで報告例があります）も研究されています。

これらの研究は、いずれも本記事で導出した基本原理（式(3d)の重み付きモンテカルロ近似、式(6)のESSによる重み集中度、勾配情報とコスト関数設計の重要性）の上に成り立っています。不確実性が大きい実システムでは、理論的な導出だけでなく、こうした実装上の工夫が実用性を大きく左右します。

関連記事

• 焼きなまし法（Simulated Annealing）の理論とPython実装 - MPPIと同じく確率的探索に基づく最適化手法を解説しています。
• ベイズ最適化の基礎とPython実装 - 代理モデルベースの効率的な最適化手法を解説しています。
• 遺伝的アルゴリズム（GA）の理論とPython実装 - 進化的計算に基づくメタヒューリスティクスを解説しています。
• PID制御の実装とチューニング - MPPIの比較対象となる古典的な制御手法を解説しています。
• クロスエントロピー法：モンテカルロ最適化の実践的手法 - MPPIのソフト重み付けと対をなすハードなエリート選択を行うCEM。同じ重点サンプリング枠組みからの導出を確認できます。
• 粒子群最適化（PSO）の理論とPython実装 - 集団ベース最適化のもう一つの代表例。MPPIのコスト重み付き平均と、PSOの速度・位置更新を対比できます。
• ガウス過程回帰の理論とPython実装 - MPPIが大量サンプリングで探索するのに対し、評価コストが高い問題ではGPサロゲートを使うベイズ最適化が有利です。
• 拡張カルマンフィルタ（EKF）の理論とPython実装 - MPPIで利用する状態推定の代表的手法。状態推定 + サンプリング制御の組み合わせが理解できます。
• 粒子フィルタのPython実装：リサンプリング手法の比較 - MPPI のソフト重み付けと同じ重点サンプリング枠組みを、状態推定の文脈で適用する補完的な手法です。
• モンテカルロ最適化メソッド比較（CEM/SA/GA/MPPI/PSO） - MPPI を含む5つのモンテカルロ最適化手法を「サンプリング→評価→更新」の共通骨格で整理するハブ記事。MPPI のソフト重み付けが他手法とどう接続するかを俯瞰できます。

参考文献

• Kappen, H. J. (2005). “Linear Theory for Control of Nonlinear Stochastic Systems.” Physical Review Letters, 95(20), 200201.
• Kappen, H. J. (2005). “Path Integrals and Symmetry Breaking for Optimal Control Theory.” Journal of Statistical Mechanics: Theory and Experiment, 2005(11), P11011.
• Theodorou, E., Buchli, J., & Schaal, S. (2010). “A Generalized Path Integral Control Approach to Reinforcement Learning.” Journal of Machine Learning Research, 11, 3137-3181.
• Williams, G., Aldrich, A., & Theodorou, E. A. (2017). “Model Predictive Path Integral Control: From Theory to Parallel Computation.” Journal of Guidance, Control, and Dynamics, 40(2), 344-357.
• Williams, G., et al. (2017). “Information Theoretic MPC for Model-Based Reinforcement Learning.” ICRA 2017.
• Alvarez-Padilla, J., Zhang, J. Z., Kwok, S., Dolan, J. M., & Manchester, Z. (2024). “Real-Time Whole-Body Control of Legged Robots with Model-Predictive Path Integral Control.” arXiv:2409.10469.
• Ryu, K., Sung, M., Gupta, P., D’sa, J., Tariq, F. M., Isele, D., & Bae, S. (2025). “IANN-MPPI: Interaction-Aware Neural Network-Enhanced Model Predictive Path Integral Approach for Autonomous Driving.” arXiv:2507.11940.
• Rubinstein, R. Y., & Kroese, D. P. (2013). The Cross-Entropy Method. Springer.