微分動的計画法

微分動的計画法 (DDP)（びぶんどうてきけいかくほう、英: Differential Dynamic Programming）は軌道最適化のために用いられる最適制御アルゴリズムの一つである。 本アルゴリズムは1966年に Mayne によって紹介され、 その後JacobsonとMayneによるその名の由来となった著作の中で分析された。 このアルゴリズムはシステムのダイナミクスとコスト関数を局所的な二次形式によってモデル化し、2次収束を示す。 また、Pantojaのstep-wise Newton法とも密接に関連している。

有限horizon離散時間問題

次式のダイナミクスを考える。

$\mathbf {x} _{i+1}=\mathbf {f} (\mathbf {x} _{i},\mathbf {u} _{i})$

(1)

これは制御入力 $\mathbf {u}$ のもとでの状態 $\textstyle \mathbf {x}$ の時刻 $i$ から $i+1$までの時間発展を表している。  総コスト $J_{0}$ は、ランニングコスト $\textstyle \ell$ と最終コスト $\ell _{f}$　の和として与えられ、状態 $\mathbf {x}$ と終端時刻までに適用する制御系列 $\mathbf {U} \equiv \{\mathbf {u} _{0},\mathbf {u} _{1}\dots ,\mathbf {u} _{N-1}\}$ によって定まる:

$J_{0}(\mathbf {x} ,\mathbf {U} )=\sum _{i=0}^{N-1}\ell (\mathbf {x} _{i},\mathbf {u} _{i})+\ell _{f}(\mathbf {x} _{N}),$

ここで $\mathbf {x} _{0}\equiv \mathbf {x}$であり、時刻 $i>0$ における状態　$\mathbf {x} _{i}$は式(1)によって与えられる。この最適制御問題の解は、上式を最小化する制御系列 $\mathbf {U} ^{*}(\mathbf {x} )\equiv \operatorname {argmin} _{\mathbf {U} }J_{0}(\mathbf {x} ,\mathbf {U} ).$ である。 これに対し、ありえる初期状態すべてを考えるのではなく、特定の初期状態 $\mathbf {x} _{0}$ に対して最適な入力系列 $\mathbf {U} ^{*}(\mathbf {x} )$ を求めるのが 軌道最適化である。

動的計画法

今、 $\mathbf {U} _{i}$ を制御入力の部分系列 $\mathbf {U} _{i}\equiv \{\mathbf {u} _{i},\mathbf {u} _{i+1}\dots ,\mathbf {u} _{N-1}\}$　とし、 残余コスト(cost-to-go) $J_{i}$　を $i$ から $N$までの部分和として定義する:

$J_{i}(\mathbf {x} ,\mathbf {U} _{i})=\sum _{j=i}^{N-1}\ell (\mathbf {x} _{j},\mathbf {u} _{j})+\ell _{f}(\mathbf {x} _{N}).$

最適な残余コスト、または時刻 $i$ での価値関数 は、以下のような制御系列による最小化で与えられる :

$V(\mathbf {x} ,i)\equiv \min _{\mathbf {U} _{i}}J_{i}(\mathbf {x} ,\mathbf {U} _{i}).$

 $V(\mathbf {x} ,N)\equiv \ell _{f}(\mathbf {x} _{N})$を出発点として、 時間を遡る方向に制御入力を一ステップ毎に最小化することにより、動的計画法の原理は全体の最小化に必要な計算量を削減する。

$V(\mathbf {x} ,i)=\min _{\mathbf {u} }[\ell (\mathbf {x} ,\mathbf {u} )+V(\mathbf {f} (\mathbf {x} ,\mathbf {u} ),i+1)].$

(2)

これがベルマン方程式である。

微分動的計画法 (DDP)

微分動的計画法 (DDP) では、新しい制御系列を計算するために、ノミナルな軌道に沿って逆方向パスの計算を行い、次に（得られた制御系列のもとでの）新しい軌道とその評価値を得るために順方向パスの計算を行うことを繰り返す。逆方向パスの計算から始めよう。

式(2)の $\min[~]$ 演算子の引数

$\ell (\mathbf {x} ,\mathbf {u} )+V(\mathbf {f} (\mathbf {x} ,\mathbf {u} ),i+1)$

に関して、その値の変分を時刻$i$の状態-入力ペア $(\mathbf {x} ,\mathbf {u} )$ 周辺でとったものを$Q$としよう:

${\begin{aligned}Q(\delta \mathbf {x} ,\delta \mathbf {u} )\equiv &\ell (\mathbf {x} +\delta \mathbf {x} ,\mathbf {u} +\delta \mathbf {u} )&&{}+V(\mathbf {f} (\mathbf {x} +\delta \mathbf {x} ,\mathbf {u} +\delta \mathbf {u} ),i+1)\\-&\ell (\mathbf {x} ,\mathbf {u} )&&{}-V(\mathbf {f} (\mathbf {x} ,\mathbf {u} ),i+1)\end{aligned}}$

これを、二次形式に展開する。

$Q(\delta \mathbf {x} ,\delta \mathbf {u} )\approx {\frac {1}{2}}{\begin{bmatrix}1\\\delta \mathbf {x} \\\delta \mathbf {u} \end{bmatrix}}^{\mathsf {T}}{\begin{bmatrix}0&Q_{\mathbf {x} }^{\mathsf {T}}&Q_{\mathbf {u} }^{\mathsf {T}}\\Q_{\mathbf {x} }&Q_{\mathbf {x} \mathbf {x} }&Q_{\mathbf {x} \mathbf {u} }\\Q_{\mathbf {u} }&Q_{\mathbf {u} \mathbf {x} }&Q_{\mathbf {u} \mathbf {u} }\end{bmatrix}}{\begin{bmatrix}1\\\delta \mathbf {x} \\\delta \mathbf {u} \end{bmatrix}}$

(3)

ここで用いる $Q$の表記は、森本淳 (ATR脳情報研究所) の記法の変形版であり、添え字は偏微分の分母を示している。 見やすくするため、添え字の$i$を省略し、次の時間ステップの変数は $V'\equiv V(i+1)$のようにプライムをつけて示そう。テイラー展開により係数は以下となる。

${\begin{alignedat}{2}Q_{\mathbf {x} }&=\ell _{\mathbf {x} }+\mathbf {f} _{\mathbf {x} }^{\mathsf {T}}V'_{\mathbf {x} }\\Q_{\mathbf {u} }&=\ell _{\mathbf {u} }+\mathbf {f} _{\mathbf {u} }^{\mathsf {T}}V'_{\mathbf {x} }\\Q_{\mathbf {x} \mathbf {x} }&=\ell _{\mathbf {x} \mathbf {x} }+\mathbf {f} _{\mathbf {x} }^{\mathsf {T}}V'_{\mathbf {x} \mathbf {x} }\mathbf {f} _{\mathbf {x} }+V_{\mathbf {x} }'\cdot \mathbf {f} _{\mathbf {x} \mathbf {x} }\\Q_{\mathbf {u} \mathbf {u} }&=\ell _{\mathbf {u} \mathbf {u} }+\mathbf {f} _{\mathbf {u} }^{\mathsf {T}}V'_{\mathbf {x} \mathbf {x} }\mathbf {f} _{\mathbf {u} }+{V'_{\mathbf {x} }}\cdot \mathbf {f} _{\mathbf {u} \mathbf {u} }\\Q_{\mathbf {u} \mathbf {x} }&=\ell _{\mathbf {u} \mathbf {x} }+\mathbf {f} _{\mathbf {u} }^{\mathsf {T}}V'_{\mathbf {x} \mathbf {x} }\mathbf {f} _{\mathbf {x} }+{V'_{\mathbf {x} }}\cdot \mathbf {f} _{\mathbf {u} \mathbf {x} }.\end{alignedat}}$

最後の三つの方程式に現れる最終項はベクトルとテンソルの縮約(contraction)を表す。 二次形式で近似された式(3) を入力$\delta \mathbf {u}$に関して最小化することで次式を得る。

${\delta \mathbf {u} }^{*}=\operatorname {argmin} \limits _{\delta \mathbf {u} }Q(\delta \mathbf {x} ,\delta \mathbf {u} )=-Q_{\mathbf {u} \mathbf {u} }^{-1}(Q_{\mathbf {u} }+Q_{\mathbf {u} \mathbf {x} }\delta \mathbf {x} ),$

(4)

これにより、オープンループ項 $\mathbf {k} =-Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} }$ とフィードバック項$\mathbf {K} =-Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} \mathbf {x} }$が与えられる。この結果を式(3) に代入することにより 時間$i$ における価値関数が得られる。

${\begin{alignedat}{2}\Delta V(i)&=&{}-{\tfrac {1}{2}}Q_{\mathbf {u} }Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} }\\V_{\mathbf {x} }(i)&=Q_{\mathbf {x} }&{}-Q_{\mathbf {u} }Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} \mathbf {x} }\\V_{\mathbf {x} \mathbf {x} }(i)&=Q_{\mathbf {x} \mathbf {x} }&{}-Q_{\mathbf {x} \mathbf {u} }Q_{\mathbf {u} \mathbf {u} }^{-1}Q_{\mathbf {u} \mathbf {x} }.\end{alignedat}}$

$i=N-1$ から$i=1$に向かって、再帰的に価値関数 $V(i)$ の局所的な二次形式のモデルと、制御の修正量 $\{\mathbf {k} (i),\mathbf {K} (i)\}$を求めるのが逆方向パスの計算である。 先に述べたように、価値関数の初期値は $V(\mathbf {x} ,N)\equiv \ell _{f}(\mathbf {x} _{N})$である。逆方向パスの完了後、新しい軌道にそって順方向パスの計算を行う。 

${\begin{aligned}{\hat {\mathbf {x} }}(1)&=\mathbf {x} (1)\\{\hat {\mathbf {u} }}(i)&=\mathbf {u} (i)+\mathbf {k} (i)+\mathbf {K} (i)({\hat {\mathbf {x} }}(i)-\mathbf {x} (i))\\{\hat {\mathbf {x} }}(i+1)&=\mathbf {f} ({\hat {\mathbf {x} }}(i),{\hat {\mathbf {u} }}(i))\end{aligned}}$

以上の逆方向パスと順方向パスの計算を収束するまで繰り返す。

正則化と直線探索

微分動的計画法は、ニュートン法と同様に、二次収束を示すアルゴリズムである。 そのため、最小値に向けて大きな修正ステップが発生するため、収束を保証するために、正則化 および/または直線探索が必要となる。 DDPの文脈における正則化とは、式(4)の行列 $Q_{\mathbf {u} \mathbf {u} }$ の正定性を保証することである。 DDPにおける直線検索は、オープンループ制御修正量 $\mathbf {k}$ をスカラー $0<\alpha <1$でスケーリングすることに対応する。

参照

• ベルマン方程式
• ハミルトン-ヤコビ-ベルマン方程式
• 動的計画法
• 最適制御

外部リンク

• PythonによるDDPの実装例
• MATLABによるDDPの実装例