微分動的計画法のソースを表示

'''微分動的計画法'''（びぶんどうてきけいかくほう、{{lang-en-short|Differential Dynamic Programming}}、略称：'''DDP'''）は軌道最適化のために用いられる最適制御アルゴリズムの一つである。本アルゴリズムは1966年に{{仮リンク|デイビッド・メイン|en|David Mayne}}によって紹介され<ref>{{Cite journal |volume=3 |pages=85–95 |last=Mayne |first=D. Q. |title=A second-order gradient method of optimizing non-linear discrete time systems |journal=Int J Control |year=1966 |doi=10.1080/00207176608921369 }}</ref>、その後JacobsonとMayneによるその名の由来となった著作の中で分析された<ref>{{cite book |last=Mayne |first= David H. and Jacobson, David Q. |title=Differential dynamic programming |year=1970 |publisher=American Elsevier Pub. Co. |location=New York|isbn=0-444-00070-4 |url=https://books.google.com/books?id=tA-oAAAAIAAJ }}</ref>。このアルゴリズムはシステムのダイナミクスとコスト関数を局所的な二次形式によってモデル化し、2次収束を示す。また、Pantojaのstep-wise Newton法とも密接に関連している<ref>{{Cite journal |doi=10.1080/00207178808906114 |issn=0020-7179 |volume=47 |issue=5 |pages=1539–1553 |last=de O. Pantoja |first=J. F. A. |title=Differential dynamic programming and Newton's method |journal=International Journal of Control |year=1988 }}</ref><ref>{{Cite journal |last=Liao |first=L. Z. |author2=C. A Shoemaker |title=Advantages of differential dynamic programming over Newton's method for discrete-time optimal control problems |publisher=Cornell University, Ithaca, NY |year=1992 |url=https://hdl.handle.net/1813/5474 }}</ref>。

== 有限horizon離散時間問題 ==
次式のダイナミクスを考える。

{{NumBlk|:|<math>\mathbf{x}_{i+1} = \mathbf{f}(\mathbf{x}_i,\mathbf{u}_i)</math>|{{EquationRef|1}}}}

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

:<math>J_0(\mathbf{x},\mathbf{U})=\sum_{i=0}^{N-1}\ell(\mathbf{x}_i,\mathbf{u}_i) + \ell_f(\mathbf{x}_N),</math>

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

== 動的計画法 ==
今、<math>\mathbf{U}_i</math> を制御入力の部分系列 <math>\mathbf{U}_i \equiv \{\mathbf{u}_i,\mathbf{u}_{i+1}\dots,\mathbf{u}_{N-1}\}</math> とし、残余コスト（cost-to-go）<math>J_i</math> を <math>i</math> から <math>N</math> までの部分和として定義する:

:<math>J_i(\mathbf{x},\mathbf{U}_i)=\sum_{j=i}^{N-1}\ell(\mathbf{x}_j,\mathbf{u}_j) + \ell_f(\mathbf{x}_N).</math>

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

:<math>V(\mathbf{x},i)\equiv \min_{\mathbf{U}_i}J_i(\mathbf{x},\mathbf{U}_i).</math>

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

{{NumBlk|:|<math>V(\mathbf{x},i)= \min_{\mathbf{u}}[\ell(\mathbf{x},\mathbf{u}) + V(\mathbf{f}(\mathbf{x},\mathbf{u}),i+1)].</math>|{{EquationRef|2}}}}

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

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

式{{EquationNote|2|(2)}}の <math>\min[~]</math> 演算子の引数

:<math>\ell(\mathbf{x},\mathbf{u}) + V(\mathbf{f}(\mathbf{x},\mathbf{u}),i+1)</math>

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

:<math>\begin{align}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{align}
</math>

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

{{NumBlk|:|<math>
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}
</math>|{{EquationRef|3}}}}

ここで用いる <math>Q</math> の表記は、森本淳（ATR脳情報研究所）の記法の変形版であり、添字は偏微分の分母を示している<ref>{{Cite conference |volume=2 |pages=1927–1932 |last=Morimoto |first=J. |author2=G. Zeglin |author3=C.G. Atkeson |title=Minimax differential dynamic programming: Application to a biped walking robot |book-title=Intelligent Robots and Systems, 2003.(IROS 2003). Proceedings. 2003 IEEE/RSJ International Conference on |date=2003
}}</ref>。

見やすくするため、添字の <math>i</math> を省略し、次の時間ステップの変数は <math>V'\equiv V(i+1)</math>のようにプライムをつけて示そう。テイラー展開により係数は以下となる。

: <math>
\begin{alignat}{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{alignat}
</math>

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

{{NumBlk|:|<math>
{\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}),
</math>|{{EquationRef|4}}}}

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

:<math>
\begin{alignat}{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{alignat}
</math>

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

:<math>
\begin{align}
\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{align}
</math>

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

== 正則化と直線探索 ==
微分動的計画法は、[[ニュートン法]]と同様に、二次収束を示すアルゴリズムである。そのため、最小値に向けて大きな修正ステップが発生するため、収束を保証するために、[[正則化]]および/または[[直線探索]]が必要となる<ref>{{Cite journal |volume=36 |issue=6 |pages=692 |last=Liao |first=L. Z |author2=C. A Shoemaker |title=Convergence in unconstrained discrete-time differential dynamic programming |journal=IEEE Transactions on Automatic Control |year=1991 |doi=10.1109/9.86943 }}</ref><ref>{{Cite thesis |publisher=Hebrew University |last=Tassa |first=Y. |title=Theory and implementation of bio-mimetic motor controllers |date=2011 |url=http://icnc.huji.ac.il/phd/theses/files/YuvalTassa.pdf }}</ref>。

DDPの文脈における正則化とは、式{{EquationNote|4|(4)}}の行列 <math>Q_{\mathbf{u}\mathbf{u}}</math> の正定値性を保証することである。DDPにおける直線検索は、オープンループ制御修正量 <math>\mathbf{k}</math> をスカラー <math>0<\alpha<1</math> でスケーリングすることに対応する。

== 脚注 ==
{{脚注ヘルプ}}
{{Reflist}}

== 関連項目 ==
* [[ベルマン方程式]]
* [[ハミルトン-ヤコビ-ベルマン方程式]]
* [[動的計画法]]
* [[最適制御]]

== 外部リンク ==
* [http://www.ros.org/wiki/color_DDP PythonによるDDPの実装例]
* [http://www.mathworks.com/matlabcentral/fileexchange/52069-ilqg-ddp-trajectory-optimization MATLABによるDDPの実装例]

{{DEFAULTSORT:ひふんとうてきけいかくほう}}
[[Category:オペレーションズリサーチ]]
[[Category:最適化アルゴリズム]]
[[Category:動的計画法]]