$\def\bm{\boldsymbol}$
概要
- ここ１０ヶ月ぐらいに渡ってオンラインで開催していた『モデルベース深層学習と深層展開』読み会が無事終わったので、これまで中身の理解がてらJuliaサンプルコードをPythonに書き直してたコードから、前回の記事までに取り扱えなかったものを晒していく
- 自動微分ライブラリにはJAXを使用する
- 今回は、最適制御問題に深層展開を使った6.3.5の内容を取り上げる
- 完全に私情にまみれた記事で恐縮だが、せっかくなのでアドベントカレンダーに応募してみた
本の内容のおさらい
- 最適制御問題については、他にたくさん解説があるので、ここでは簡本で触れられてるレベルで簡単に復習しておく
- 制御対象の動特性や入力エネルギーに関する指標を設計し、その指標を最小化（もしくは最大化）するように入力を決める問題が最適制御問題
- この指標は状態量を$\bm{x}(t)$、入力量を$\bm{u}(t)$とすれば、例えば以下のように書ける
- $J[\bm{u}(t)]\equiv \phi(\bm{x}(t)+\int^T_0F(\bm{x}(t),\bm{u}(t))dt$
- この$J$は$\phi$と$F$で特徴付けられた汎関数である。この汎関数の停留値を与えるような関数$\bm{u}(t)$を見つけるという意味で、この問題はひとつ前の節で扱った変分問題の１種である
- モデルベースの制御手法では$\bm{x}$のダイナミクス$\dot{\bm{x}}=h(\bm{x}(t),\bm{u}(t))$は何らかの方法でモデリングし、これを$J$の評価に陽に用いる
- $J[\bm{u}(t)]\equiv \phi(\bm{x}(t)+\int^T_0F(\bm{x}(t),\bm{u}(t))dt$
- この指標は状態量を$\bm{x}(t)$、入力量を$\bm{u}(t)$とすれば、例えば以下のように書ける
- 制御対象の動特性や入力エネルギーに関する指標を設計し、その指標を最小化（もしくは最大化）するように入力を決める問題が最適制御問題
プログラムでの理解
問題設定
本で出てきた簡単な例は以下である
- ダイナミクスが以下となるような１次元の系を制御対象とする
- $\dot{x}=u(t)$
- $J$を以下と仮定して$J$を最小化する最適な$u(t)$を導きたい
- $J[u(t)]=\int^T_0\{x(t)^2+u(t)^2\}dt$
- ただし、$x(0)=1$という境界条件が与えられている
- これを解く方法は解析的、近似的合わせて種々ある
- 一番真っ当（？）なのは最適性の原理に従って、ハミルトン・ヤコビ・ベルマン方程式と呼ばれる偏微分方程式を解く問題に帰着して、動的計画法などの解法で解く方法であるが、この教科書では、$u$を決定する関数を万能近似器で近似し、$J$の数値解を起点に誤差逆伝搬で万能近似のパラメータを更新するという大胆（？）な方法で$u(t)$を導く
- このとき、$J$の計算は時間方向に展開された数値積分で行うので、深層展開の文脈で扱われている
- 本にはこの問題の数値解が以下のようになることが知られていると書いてある
- $u(t)=^\frac{\sinh(1-t)}{\cosh(1)}$
- 本中には導出過程は書かれておらず、引用論文が載っている
- だた、この論文中にも導出方法は書かれていない…？
- 時間があれば導出したい。。※要更新箇所
- 本中には導出過程は書かれておらず、引用論文が載っている
- $u(t)=^\frac{\sinh(1-t)}{\cosh(1)}$
- 一番真っ当（？）なのは最適性の原理に従って、ハミルトン・ヤコビ・ベルマン方程式と呼ばれる偏微分方程式を解く問題に帰着して、動的計画法などの解法で解く方法であるが、この教科書では、$u$を決定する関数を万能近似器で近似し、$J$の数値解を起点に誤差逆伝搬で万能近似のパラメータを更新するという大胆（？）な方法で$u(t)$を導く
Pythonで実装
必要ライブラリインポート
import jax
import jax.numpy as jnp
from jax.example_libraries import optimizers
from tqdm.notebook import trange
from functools import partial
準備
delta = 0.05
beta = 20.0
N = 100
T = 1.0
放射基底関数の定義
def rbf(x, beta):
return jnp.exp(-beta*x**2)
@jax.jit
def rbf_func(x, theta):
s = 0.
for i in range(len(theta)):
s +=theta[i] * rbf(x - delta*i+0.5, beta)
return s
- ガウス放射基底関数$\phi(x)=\exp(-\beta x^2)$を使って、最適制御則の近似器を
- $u_{\bm{\theta}}(t)=\sum_{i}\theta_i\phi(t-\delta\cdot i + \frac{1}{2})$とおく
- $\beta$と$\delta$はハイパラ
- $u_{\bm{\theta}}(t)=\sum_{i}\theta_i\phi(t-\delta\cdot i + \frac{1}{2})$とおく
微分可能数値積分モジュールの定義
def Integrate(T, N, theta):
dt = T/N
J = 0.
x = 1.
for i in range(N):
u = rbf_func(i*dt, theta)
x += u*dt
F = x**2 + u**2
J += F*dt
return J
- 区分求積法を用いて$J$を計算
学習
train_itr = 200
adam_lr = 0.1
opt_init, opt_update, get_params = optimizers.adam(adam_lr)
def step(T, N, step_num, opt_state):
value, grads = jax.value_and_grad(Integrate, argnums=-1)(T, N, get_params(opt_state))
new_opt_state = opt_update(step_num, grads, opt_state)
return value, new_opt_state
def train(T, N, theta):
opt_state = opt_init(theta)
for itr in trange(train_itr, leave=False):
value, opt_state = step(T, N, itr, opt_state)
print("\r"+"\rloss:{}".format(value), end=" ")
return get_params(opt_state)
theta_init = jnp.ones(50)
theta_trained = train(T, N, theta_init)
- これまでの章で出てきた深層展開の計算とやってることは大差ないが、今回は数値積分の繰り返し計算が、展開対象になっていることに注意
- 学習対象は$u_{\bm{\theta}}$のパラメタ$\bm{\theta}$
結果（厳密解との比較）
- 教科書と同様、厳密解をよく近似できている
その他
- 今回の手法は、冷静に考えてみると$u(t)$を大胆に近似している結構強引な手法に思える。
- もっと複雑な系になればモデル化誤差の影響を大きく受けるだろうし、簡単な系においても、線形システムに近似する制御手法のように強力な安定性解析が使えるわけでもない。
- 実際のトコどんくらい利用価値があるのかちょっとわからん…
- もっと複雑な系になればモデル化誤差の影響を大きく受けるだろうし、簡単な系においても、線形システムに近似する制御手法のように強力な安定性解析が使えるわけでもない。
バックナンバー
- 『モデルベース深層学習と深層展開』読み会 レポート（開催前準備編）
- 『モデルベース深層学習と深層展開』読み会レポート#0
- 『モデルベース深層学習と深層展開』読み会レポート#1
- 『モデルベース深層学習と深層展開』読み会レポート#2
- 『モデルベース深層学習と深層展開』読み会レポート#3
- 『モデルベース深層学習と深層展開』読み会レポート#4
- 『モデルベース深層学習と深層展開』読み会レポート#5
- 『モデルベース深層学習と深層展開』読み会レポート#6
- 『モデルベース深層学習と深層展開』読み会レポート#7
- 『モデルベース深層学習と深層展開』読み会レポート#8
- 『モデルベース深層学習と深層展開』読み会レポート#9
- 『モデルベース深層学習と深層展開』読み会レポート#10
- 『モデルベース深層学習と深層展開』読み会レポート#11