離散時間システムの最適制御
記述をシンプルにするため, 変数の属する次元や添字の範囲は省略する. 与えられる関数は適宜滑らかさを課す. ときどき思い出したように状態ベクトル x, 制御ベクトル u の次元はそれぞれ n,m であることを用いる.
非線形離散時間システムの制御問題
問題設定
終端時刻を自然数 N と書き, 時刻を表す変数 k=0,\cdots, N に対して, 以下の離散制御システムを考える;
x(k+1) = f(x(k),u(k),k)
ここで x(k), u(k) はそれぞれ時刻 k=0,\cdots, N におけるシステムの状態ベクトル, 制御入力ベクトルを表す.
このシステムに対して, 初期状態 x(0) = x_0 が固定されたときの評価関数
J(x, u) = \varphi(x(N)) + \sum_{k=0}^{N-1} L(x(k), u(k), k)
を最小化するような制御入力 \{u(k)\}_{k=0}^{N-1} を構成したい(u(N) の値は評価関数には影響しないので, k=N-1 までで十分である).
オイラーラグランジュ方程式の導出
差分方程式を等式拘束条件とみなすと, 最小化問題のラグランジュ関数はラグランジュ乗数 \{\lambda(k)\}_{k=1}^N を用いて
\begin{aligned} \bar{J}(x,u) &= \varphi(x(N)) - \lambda^T(N)x(N) \\ &\quad + \sum_{k=1}^{N-1}(H(x(k), u(k), \lambda(k+1), k)-\lambda^T(k)x(k))\\ &\quad + H(x(0), u(0),\lambda(1),0) \end{aligned}
と表せる. ただしここに現れた関数 H(x,u,\lambda,k) := L(x,u,k) +\lambda^Tf(x,u,k) をハミルトニアンといい, 以下では H(k) = H(x(k), u(k), \lambda(k+1),k) と略記する.
初期条件 x(0)=x_0 が拘束されているため, \bar{J} は (x, u)\in ({x_0}\times (\mathbf{R}^{n})^{(N-1)})\times (\mathbf{R}^{m})^N を定義域に持つ関数である.
最小化問題の必要条件となる \bar{J} の停留条件を求めよう. 方向ベクトル (y, v)\in ({0}\times (\mathbf{R}^{n})^{(N-1)})\times (\mathbf{R}^{m})^N とパラメータ h>0 を用いると
\begin{aligned} &\bar{J}(x+hy,u+hv) - \bar{J}(x,u)\\ &\quad= \varphi(x(N)+hy(N)) - \varphi(x(N)) - h\lambda^T(N)y(N)\\ &\qquad + \sum_{k=1}^{N-1}(H(x(k)+hy(k), u(k)+hv(k), \lambda(k+1), k) - H(x(k), u(k), \lambda(k+1), k)-h\lambda^T(k)y(k))\\ &\qquad + H(x(0), u(0)+hv(0),\lambda(1),0) -H(x(0), u(0),\lambda(1),0) \end{aligned}
と表せる. ハミルトニアン H やコスト関数 \varphi の微分を用いると, ラグランジュ関数の方向微分は
\begin{aligned} d\bar{J}((x,u);(y,v)) &= \left(\dfrac{\partial \varphi}{\partial x}(x(N))-\lambda^T(N)\right) y(N)\\ &\quad + \sum_{k=1}^{N-1}\left(\dfrac{\partial H}{\partial x}(k) - \lambda^T(k)\right)y(k) + \sum_{k=1}^{N-1}\dfrac{\partial H}{\partial u} v(k)\\ &\quad + \dfrac{\partial H}{\partial u}v(0) \end{aligned}
と表せる. これにより, 停留条件として以下のオイラー・ラグランジュ方程式が得られる
\begin{aligned} &x(k+1) = f(x(k), u(k), k), \quad &x(0)=x_0, \\ &\lambda(k) = \left(\dfrac{\partial H}{\partial x}\right)^T (x(k), u(k), \lambda(k+1), k), \quad &\lambda(N) = \left(\dfrac{\partial \varphi}{\partial x}\right)^T (x(N)), \\ &\dfrac{\partial H}{\partial u}(x(k), u(k), \lambda(k+1), k)=0. \end{aligned}
ここで, ハミルトニアンの定義から状態方程式は
x(k+1) = \left(\dfrac{\partial H}{\partial \lambda}\right)^T (x(k), u(k), \lambda(k+1), k)
とも表せるため, 上で得られた状態方程式は解析力学での正準方程式に似ている. 随伴変数 \lambda の式(随伴方程式という)にマイナスがつかないのは, 時間逆向きの差分方程式となっているため.
上記の連立差分方程式は状態方程式で初期条件, 随伴方程式で終端条件が与えられているため, 2点境界値問題と呼ばれる. 拘束条件 \frac{\partial H}{\partial u}(x(k), u(k), \lambda(k+1), k)=0 は H のヘッセ行列 \frac{\partial^2 H}{\partial u^2} が正則であれば u(k) について解け, これを用いて x(k), \lambda(k) について反復法により解くことができる.
ベルマン方程式の導出
まず, x_k\in \mathbf{R}^n, k=0,\cdots, N について, x(k)=x_k, x(k+1)=f(x(k), u(k),k) となるようにする. このとき, 値関数 V(x,k) を次のように定める:
V(x_k,k) = \min_{\{u(\ell)\}_{\ell=k}^{N-1}} \left( \varphi(x(N)) + \sum_{\ell=k}^{N-1}L(x(\ell), u(\ell), \ell) \right).
これは x(k)=x_k 以降の評価関数の最小値と, その最小値を達成する \{u(\ell)\}_{\ell=k}^{N-1} を求める最小化問題である. ここでは簡単のため, 最小値は達成されるものとする. 以下のことは直ちに分かる; \begin{aligned} V(x_0,0) &= \min_{x,u} J(x,u),\\ V(x_N, N) &= \varphi(x_N), \quad & x_N\in \mathbf{R}^n \end{aligned}
さらに, 値関数の右辺のうち L(x(k), u(k), k) はそれ以降の入力には依存しないことに注意すると \begin{aligned} V(x_k,k) &= \min_{\{u(\ell)\}_{\ell=k}^{N-1}} \left( L(x(k), u(k), k) + \varphi(x(N)) + \sum_{\ell=k+1}^{N-1}L(x(\ell), u(\ell), \ell) \right)\\ &= \min_{u(k)} \left( L(x(k), u(k), k) +\min_{\{u(\ell)\}_{\ell=k+1}^{N-1}} \left( \varphi(x(N)) + \sum_{\ell=k+1}^{N-1}L(x(\ell), u(\ell), \ell) \right) \right)\\ &= \min_{u(k)} \left( L(x_k, u(k), k) + V(f(x_k,u(k),k),k+1) \right) \end{aligned}
となる. これより, 値関数は
\begin{aligned} V(x,k) &= \min_{u} \left( L(x, u, k) + V(f(x,u,k),k+1) \right), \quad & x\in \mathbf{R}^n, k=0,\cdots, N-1,\\ V(x, N) &= \varphi(x), \quad & x\in \mathbf{R}^n \end{aligned}
という後ろ向きの漸化式の解となる. この漸化式をベルマン方程式と呼び, このように評価関数の最小値を再帰的に表す手法を動的計画法と言う.
以上により, 値関数が存在すれば, ベルマン方程式が成り立つことが分かる.
逆に, ベルマン方程式が解けた場合, つまり値関数 V(x,k) と optimizer u_{\mathrm{opt}}(x,k) が得られた場合, \bar{x}(k+1) = f(\bar{x}(k), u_{\mathrm{opt}}(x(k),k), k) という閉ループフィードバック制御によって最適制御を構成できることが値関数の定義からわかる.
ベルマン方程式からのオイラー・ラグランジュ方程式の導出
ベルマン方程式 \begin{aligned} V(x,k) &= \min_{u}\left(L(x,u,k)+ V(f(x,u,k), k+1)\right), \quad & x\in \mathbf{R}^n, k=0,\cdots N-1,\\ V(x,N) &= \varphi(x), \quad & x\in \mathbf{R}^n \end{aligned}
が任意の x, k に対して解けたとし, u_{\mathrm{opt}}(x,k) を optimizer とする. さらに, \bar{x} を, 制御入力 \bar{u}(k) = u_{\mathrm{opt}}(\bar{x(k)},k) としたときの閉ループシステムの状態ベクトルとする. すなわち, 次のように構成する;
\begin{aligned} & \bar{x}(0) = x_0\\ & \bar{u}(k) = u_{\mathrm{opt}}(\bar{x}(k), k), \quad & k=0,\cdots N,\\ & \bar{x}(k+1) = f(\bar{x}(k), \bar{u}(k), k), \quad & k=0,\cdots N-1. \end{aligned}
このとき, (\bar{x}, \bar{u}) はオイラー・ラグランジュ方程式の解となることを見ていこう.
そのためには, (\bar{x}, \bar{u}) はオイラー・ラグランジュ方程式の状態方程式を満たすことは明らかなので, 随伴方程式を満たす \bar{\lambda} を構成し, 拘束条件を満たすことを確かめれば良い.
まず, 任意の x, k=0,\cdots, N-1 について,
\begin{aligned} \alpha(x, k) & := -V(x,k) + L(x,\bar{u}(k), k) + (f(x,\bar{u}(k),k), k+1)\\ & \ge -V(x,k) + \min_{u}\left(L(x,u, k) + (f(x,u, k+1), k )\right)\\ &= 0 \end{aligned}
が値関数であることから従う. さらに, ベルマン方程式から \alpha(\bar{x}(k),k)=0 を満たすので,
\begin{aligned} 0 &= \left.\dfrac{\partial \alpha}{\partial x}(x,k)\right|_{x=\bar{x}(k)}\\ &= - \dfrac{\partial V}{\partial x}(\bar{x}(k), k) + \dfrac{\partial L}{\partial x}(\bar{x}(k), \bar{u}(k), k) + \dfrac{\partial V}{\partial x}(f(\bar{x}(k), \bar{u}(k),k), k+1) \dfrac{\partial f}{\partial x}(\bar{x}(k), \bar{u}(k), k)\\ &= - \dfrac{\partial V}{\partial x}(\bar{x}(k), k) + \dfrac{\partial L}{\partial x}(\bar{x}(k), \bar{u}(k), k) + \dfrac{\partial V}{\partial x}(\bar{x}(k+1), k+1) \dfrac{\partial f}{\partial x}(\bar{x}(k), \bar{u}(k), k)\\ &= - \dfrac{\partial V}{\partial x}(\bar{x}(k), k) + \dfrac{\partial H}{\partial x} \left( \bar{x}(k), \bar{u}(k), \left(\dfrac{\partial V}{\partial x}\right)^T(\bar{x}(k+1), k+1), k \right) \end{aligned}
が成り立つ. さらに, ベルマン方程式の終端条件により,
\dfrac{\partial V}{\partial x}(\bar{x}(N), N) = \dfrac{\partial \varphi}{\partial x}(\bar{x}(N))
も成り立つ. 以上を2式を合わせると,
\bar{\lambda}(k) = \dfrac{\partial V}{\partial x}(\bar{x}(k), k), \quad k= 0, \cdots N
と定めれば, オイラー・ラグランジュ方程式の随伴方程式を満たす.
さらに, ベルマン方程式を再び用いることで, 任意の u, k=0,\cdots, N について,
\begin{aligned} \beta(u, k) & := L(\bar{x}(k),u, k) + (f(\bar{x}(k),u,k), k+1)\\ & \ge L(\bar{x}(k), \bar{u}(k), k) + (f(\bar{x}(k),\bar{u}(k),k), k+1)\\ &= 0 \end{aligned}
を定めると, \beta(\cdot, k) は \bar{u}(k) で最小値を取ることから
\begin{aligned} 0 &= \left.\dfrac{\partial \beta}{\partial u}(x,k)\right|_{u=\bar{u}(k)}\\ &= \dfrac{\partial L}{\partial u}(\bar{x}(k), \bar{u}(k), k) + \dfrac{\partial V}{\partial u}(f(\bar{x}(k), \bar{u}(k),k), k+1) \dfrac{\partial f}{\partial u}(\bar{x}(k), \bar{u}(k), k)\\ &= \dfrac{\partial L}{\partial u}(\bar{x}(k), \bar{u}(k), k) + \dfrac{\partial V}{\partial u}(\bar{x}(k+1), k+1) \dfrac{\partial f}{\partial u}(\bar{x}(k), \bar{u}(k), k)\\ &= \dfrac{\partial H}{\partial u} \left( \bar{x}(k), \bar{u}(k), \left(\dfrac{\partial V}{\partial x}\right)^T(\bar{x}(k+1), k+1), k \right)\\ &= \dfrac{\partial H}{\partial u} \left( \bar{x}(k), \bar{u}(k), \bar{\lambda}(k+1), k \right) \end{aligned}
が成り立つ. 以上により, (\bar{x},\bar{\lambda},\bar{u}) はオイラー・ラグランジュ方程式の解である.
数値解法
離散LQ制御問題
システムは次のような線形システム
x(k+1) = Ax(k) + B u(k)
を仮定し, 評価関数は次の二次形式で与えられるものとする
J=\dfrac12 x^T(N) S_fx(N) + \sum_{k=1}^{N-1}\dfrac12\left(x^T(k)Qx(k) + u^T(k)Ru(k)\right).
ただし簡単のため S_f, Q, Rは対称行列とし, 特にRは正定値とする.
オイラー・ラグランジュ方程式
ハミルトニアンは時間には依存せずH(x,u,\lambda) = \frac12(x^TQx + u^TRu)+ \lambda^T(Ax+Bu)となるのでオイラー・ラグランジュ方程式は
\begin{aligned} &x(k+1) = Ax(k) + Bu(k), \quad x(0)=x_0, \\ &\lambda(k) =Qx(k)+A^T \lambda(k+1) , \quad \lambda(N) = S_f (x(N)), \\ &u^T(k)R + \lambda^T(k+1)B=0. \end{aligned}
となる. ここで, 関係式\lambda(k)=S(k)x(k)を代入してS(k)の条件を求めると
\begin{aligned} S(k) &=Q+A^TS(k+1)A\\ &\quad - A^TS(k+1)B(R+B^TS(k+1)B)^{-1}B^TS(k+1)A,\quad k=1,\cdots, N-1,\\ S(N)&=S_f \end{aligned}
を得る. これをリッカチ方程式と言い, 条件R>0からS(N), S(N-1), \cdots, S(1)と順次解くことができる.
また, オイラー・ラグランジュ方程式の第1式と第3式から
u(k) = -(R+B^TS(k+1)B)^{-1}B^TS(k+1)Ax(k)
を得る. 以上により, 以下のステップで最適制御と状態遷移を求めることができる
- S(N), S(N-1), \cdots, S(1)をリッカチ方程式を解いて求める
- u(k) = -(R+B^TS(k+1)B)^{-1}B^TS(k+1)Ax(k)から, x(0), S(1) を用いてu(0)を定める
- 状態方程式x(k+1)=Ax(k)+Bu(k)からx(0), u(0)を用いてx(1)を求める
- 以上を繰り返して, u(1), x(1), u(2),\cdotsを順次求める
ベルマン方程式
ベルマン方程式は
\begin{aligned} V(x,k) &= \min_{u}\left(\dfrac12(x^TQx+u^TRu)+ V(Ax ＋Bu, k＋１)\right), \quad & x\in \mathbf{R}^n, k=0,\cdots N-1,\\ V(x,N) &= \dfrac12 x^TS_fx, \quad & x\in \mathbf{R}^n \end{aligned}
となる. 発見的に対称行列S(k)とおいてV(x,k)=\dfrac12 x^TS(k)xを代入すると, ベルマン方程式の右辺の min の内部をuについて平方完成することにより,
\begin{aligned} &\dfrac12(x^TQx+u^TRu)+ \dfrac{1}{2}(Ax+Bu)^TS(k+1)(Ax+Bu)\\ &\quad = \dfrac12( u^T(R+B^TS(k+1)B)u + 2u^TB^TS(k+1)Ax + x^T(Q+A^TS(k+1)A)x)\\ &\quad = \dfrac12\bigg( (u+B^TS(k+1)Ax)(R+B^TS(k+1)B)(u+B^TS(k+1)Ax)\\ &\qquad \quad +x^T(Q+A^TS(k+1)A-A^TS(k+1)B(R+B^TS(k+1)B)^{-1}B^TS(k+1)A)x \bigg). \end{aligned}
ただしここで平方完成は
u^TAu+2u^Tv = (u+v)^TA(u+v) - v^TA^{-1}v
を用いた. 右辺より, 値関数の min を実現する optimizer u は u = -(R+B^TS(k+1)B)^{-1}B^TS(k+1)Ax
であり, このとき
\begin{aligned} \dfrac12 x^TS(k)x &= \dfrac12 x^T\left\{ Q+A^TS(k+1)A\right.\\ &\left.\qquad \qquad - A^TS(k+1)B(R+B^TS(k+1)B)^{-1}B^TS(k+1)A \right\}x \end{aligned}
を得る. 任意のxについてこれが成り立つには,
\begin{aligned} S(k) &=Q+A^TS(k+1)A\\ &\quad - A^TS(k+1)B(R+B^TS(k+1)B)^{-1}B^TS(k+1)A\\ \end{aligned}
を得る. これはリッカチ方程式である.
以上により, 離散LQ制御問題においては, オイラー・ラグランジュ方程式及びベルマン方程式はリッカチ方程式を経由して解くことができ, どちらも同じ最適制御を与える.
他の問題設定
終端条件固定(x(N)が与えられた場合)
オイラー・ラグランジュ方程式の導出の際に, 変分について条件y(N)=0が追加される. これにより,
以下のオイラー・ラグランジュ方程式が得られる
\begin{aligned} &x(k+1) = f(x(k), u(k), k), \quad &x(0)=x_0,x(N)=x_N \\ &\lambda(k) = \left(\dfrac{\partial H}{\partial x}\right)^T (x(k), u(k), \lambda(k+1), k),\\ &\dfrac{\partial H}{\partial u}(x(k), u(k), \lambda(k+1), k)=0. \end{aligned}
終端条件が固定されているので, \varphiは現れず, \lambda(N)の代わりにx(N)についての条件が与えられる.