MPPIのお気持ち

モデル予測経路積分（Model Predictive Path Integral, MPPI）は実装も簡単で動作も速い便利なMPC/最適制御ソルバとして注目されている．しかしその裏側にある理論を見てみると分かりにくいところがちらほらある．ここでは，MPPIの理論について分かりやすく説明することを試みる．
* **対象読者**：制御分野の人で，MPPIの理論を調べようとしたけどよく分からなかった方．
* 参考文献としては，[Information Theoretic Model Predictive Control: Theory and Applications to Autonomous Driving](https://arxiv.org/abs/1707.02342)などがよい．
* ここでは簡単のため，ホライズンを$H=1$，入力次元数を1として最適な入力$u\in\mathbb R$を求める問題を扱う．

# 登場する確率分布
シンボル|名前|形|役割
---|---|---|---
$P$|ノミナル分布|正規分布$\mathcal N(u_\mathrm{nom},\sigma^2)$に限る|$Q^\ast$を作るための種．$u_\mathrm{nom}$はノミナル入力で，これから離れる入力$u$ほど入力コストが高くなる．
$Q^*$|最適分布|複雑|もしもこの分布の元で制御できれば総コストが最小になる．
$Q(u)$|調整用分布|正規分布$\mathcal N(u,\sigma^2)$に限る|$u$を調整することで複雑な$Q^\ast$を近似する．我々が一番知りたいのは，$u$の最適値である．
$G$|サンプル生成器|何でもいい|$Q(u)$の最適な$u$をモンテカルロ法で求める時に使うサンプルを生成する．

> 各分布の確率・確率密度関数は分布シンボルの小文字とする．例えば$P$の確率・確率密度関数は$p(v)$．


# 最適制御問題
入力$u$によって更新された状態に対する状態コストを$S(u)$と表すこととする．また，入力コスト$S_u(u)$を次のように定める：
```math
S_u(u) = \frac{\lambda}{2\sigma^2}(u-u_\mathrm{nom})^2.
```

:::note warn
**伏線：変な係数**
入力コスト$S_u(u)$には$\frac{\lambda}{2\sigma^2}$という変な係数が入っているがこれは伏線である．
:::

また，総コストを
```math
J(u) := \mathbb E_{v \sim Q(u)}[S(v)] + S_u(u)
```
と定める．これを最小化することを今回の主問題とする：
```math
\min_u J(u).
```
:::note warn
何と名付けよう．最適化問題？うーん，分布を最適化するから分布最適化問題とか？
:::


# 自由エネルギーと最適分布の導出
ノミナル分布$P$から$v$が生成されている状況を考える．この時，次のような自由エネルギーを天下り的に定義する：
```math
F := -\lambda \log \mathbb E_{v\sim P} \left[ 
    e^{-S(v)/\lambda} 
\right] .
```
> ここで$\lambda$は入力コスト$S_u$にも出てきた定数である．

ここでどっかからやってきた別の分布$Q$を想定し，自由エネルギーを記述すると：
```math
F = -\lambda \log \mathbb E_{v\sim Q} \left[
    e^{-S(v)/\lambda} \ \frac{p(v)}{q(v)}
\right] .
```
ところで自由エネルギー$F$にはイェンセン不等式が成り立つということで，
```math
\begin{align}
F 
& \le -\lambda \mathbb E_{v\sim Q} \left[
\log\left(
    e^{-S(v)/\lambda} \ \frac{p(v)}{q(v)}
\right) \right] \\
& = \mathbb E_{v\sim Q(u)} [S(v)] + \lambda \mathbb D_\mathrm{KL}(Q(u)||P) .
\end{align}
```

> 最後に出てきた"$\mathbb E_{v\sim Q(u)} [S(v)] + \lambda \mathbb D_\mathrm{KL}(Q(u)||P)$"のことは**正則化項付き状態コスト期待値**とでも呼ぼう．

この自由エネルギーイェンセン不等式に関して面白い事実が2つある：

:::note info
正則化項付き状態コスト期待値は，実は**総コスト$J(u)$そのものである**．何故なら，入力コスト部分が
```math
\lambda \mathbb D_\mathrm{KL}(Q(u)||P) = \frac{\lambda}{2} \ \frac{(u-u_\mathrm{nom})^2}{\sigma^2}
```
となるからである（伏線回収）．
:::

:::note info
```math
Q^\ast: \ q^\ast(v) := \eta^{-1} e^{-S(v)/\lambda} p(v)
```
という天下りな分布$Q^\ast$を$Q=Q^\ast$に代入すると，等式が成り立つ．証明してみてね⭐️
:::

$Q=Q^\ast$とすると，自由エネルギーイェンセン不等式は等式となり，正則化項付き状態コスト期待値，改め，総コストは最小を取る．その通り，$Q^\ast$こそ**最適分布**である．

# 最適分布の近似
ここからは，なるべく最適分布$Q^\ast$に近い調整用分布$Q(u)$を作るための最適な$u$を求める問題を考える．具体的には次のような最適問題となる：
```math
\min_u \mathbb D_\mathrm{KL}
(Q^\ast || Q(u)) .
```

:::note warn
この問題の最適解が$\min_u J(u)$の最適解と同じという保証はない（よね？）
:::

この問題の最適解は
```math
\begin{align}
u^* &= \mathbb E_{v\sim Q^\ast}[v]
\end{align}
```
となる．証明してみてね⭐️

:::note info
**最適解がこの形に収まるのはかなり都合がいい**．なぜかというと，モンテカルロ法などで$v\sim Q^\ast$を生成できれば，その平均がすなわち最適な$u^\ast$となるからである．
最適入力を求めるのに，よく分からない最適化ソルバを引っ張り出してきて，いつ終わるか分からない反復求解アルゴリズムを実行しなくていいのだ．
:::

なーんだ，じゃあ$Q^\ast$からモンテカルロサンプリングすればいいじゃん！と思うかもしれないが待ってほしい．今のところ，$Q^\ast$は不明な分布であり，生成器を実装できないのだ．

# 別分布からのモンテカルロサンプリング
ここからの問題は
```math
\begin{align}
u^* &= \mathbb E_{v\sim Q^\ast}[v]
\end{align}
```
をモンテカルロ法で求めることであるが，$Q^\ast$の生成器は実装不可能なので，別の分布から行なわねばならない．そこで使われるのがサンプル生成器$G$である．
$G$は何でもよいのだが，$G$の種類によって$u^*$を求める数式がちょっとずつ違っている．実際，いくつかのMPPIに関する論文をあたると最後の$u^\ast$を求める数式がちょっとずつ違っており，混乱のもとになっている．

### 一般的な議論
モンテカルロ法において，別の分布のサンプルを用いて所望の期待値を取り出すには，重みづけというのが必要になる：
```math
\begin{cases}
u^* &= \mathbb E_{v\sim Q^\ast}[v] = \mathbb E_{v\sim G}\left[ w(v) v \right] \\
w(v) &= \frac{q^\ast(v)}{g(v)}
\end{cases}
```
となるため，十分大きな$N$を用いて
```math
u^\ast \approxeq \sum_{i=1}^N \frac{w(v_i)}{\sum_{j=1}^N w(v_j)} v_i ,\quad \forall i\in\{1,\ldots,N\}, v_i\sim G
```
とする必要がある．サンプル生成器$G$によってこの重み$w(v)$の表式が変わってくる．

### サンプル生成器＝任意の正規分布
$G = \mathcal N(\hat u, \sigma^2)$とした時（ただし$\hat u\ne u_\mathrm{nom}$），
```math
w(v)=\frac{q^\ast(v)}{Z^{-1}\exp((u-\hat u)^2/\sigma^2)}
\propto \exp \left(
-\frac{1}{\lambda}S(v)-(\hat u-u_\mathrm{nom})v/\sigma^2
\right)
```
となる（導出しみてね🌟）．

:::note info
$\hat u$としては，予想される最適解$u^\ast$に近いほうが精度が高い（モンテカルロ法の結果の分散が小さい）．例えばMPC制御であれば，前ステップの入力$u^\ast$を入れたりするのがよいだろう．
:::

このサンプル生成器として話が進んでいる文献は例えば
* [Information Theoretic Model Predictive Control: Theory and Applications to Autonomous Driving](https://arxiv.org/abs/1707.02342) 式22
* [モデル予測経路積分制御によるクワッドロータの制御](https://www.jstage.jst.go.jp/article/isciesci/64/8/64_316/_article/-char/ja) 式33
    * ノミナル入力$u_\mathrm{nom}=0$としている．
    * $w(v)$の中に$\hat u^2/2\sigma^2$という定数項があるが，これを外の規格化定数に追い出してしまえば，私が示したものと同じになる．
* [Contingency Constrained Planning with MPPI within MPPI](https://arxiv.org/abs/2412.09777) 式7
* [モデル予測経路積分制御 (Model Predictive Path Integral Control; MPPI) 入門](https://zenn.dev/takuya_fukatsu/articles/36c0d6911d18b7) 2章

### サンプル生成器＝ノミナル分布と同じ
$G=P=\mathcal N(u_\mathrm{nom},\sigma^2)$とするものである．この時は単純になるので導出するけど，
```math
\begin{align}
w(v)&=\frac{q^\ast(v)}{p(v)} = \frac{\eta^{-1} e^{-S(v)/\lambda} p(v)}{p(v)} \\
&=\eta^{-1} e^{-S(v)/\lambda} .
\end{align}
```
このサンプル生成器として話が進んでいる文献は例えば
* [Real-Time Whole-Body Control of Legged Robots with Model-Predictive Path Integral Control](https://arxiv.org/abs/2409.10469) 式1

# さいごに
数学の理論で分かんなくなったりした場合，LLMが結構役に立つのでどしどし議論するとよい．

▼これはCopilotに別分布からのモンテカルロサンプリングについての理論を質問した例．

https://copilot.microsoft.com/shares/xoHSnwAP5AzeNZJF5oSrN

あと，何か間違えている箇所があれば指摘して頂ければと思います．