制御器と学習エージェントが通信ネットワークを介してプラントと接続しているシステムを考える。つまり, 制御器は学習エージェントよって作られた制御方策と通信則をもとにプラントを操作するために制御入力を送信するという働きをする。ただし, パケットロス, 量子化誤差, 通信遅延は本記事では考慮しない。

制御対象とするプラントの設定

制御対象であるプラントのダイナミクスを次のように与える。

\begin{align}
x_{k+1} &= f(x_k, u_k) + \varepsilon,\ \ u_k\in U \\
(x_{i, k+1} &= f_i(x_k, u_k) + \varepsilon, \ \ i\in\mathbb{N}_{1:n_x})
\end{align}
\tag{1}

ただし, $k \in \mathbb{N}_{\ge 0}$, $x_k \in \mathbb{R}^{n_x}$ は状態, $u_k\in\mathbb{R}^{n_u}$は制御入力, $U \subset R^{n_u}$ は制御入力の集合, $f:\mathbb{R}^{n_x}\times \mathbb{R}^{n_u}\to \mathbb{R}^{n_x}$は未知の遷移関数, $\varepsilon \sim \mathcal{N}(0, \sigma_{\varepsilon}^2)$とする。

関数$f$を構成する要素$f_i, i \in \mathbb{N}_{1:n_x}$ がガウス過程(GP)回帰によって学習し, モデリングする。ただし, カーネル関数として以下のようにガウスカーネルを選択する。　

k_i(x, x') = \alpha_i^2 {\rm exp}(-\frac{1}{2}(x - x')^\top \Lambda_i^{-1} (x - x'))
\tag{2}

ただし, $x, x'\in \mathbb{R}^{n_x + n_u}$の時, $\Lambda_i= {\rm diag}(\lambda_{i, 1}^2, \ldots, \lambda_{i, n_x + n_u}^2)$で$\alpha_i$, $\lambda_{i, n}$ $i\in\mathbb{N}_{1:n_x}$, $n\in\mathbb{N}_{1:n_x + n_u}$は$f_i$のGPモデルに対するハイパーパラメータとする。また, 入出力学習データ $\mathcal{D}_i = \bigl\{X, y_i\bigr\}$について, $X, y_i, i \in \mathbb{N}_{1:n_x}$は次のように与える。

X =  \left[
    \begin{array}{rrrr}
      x_0^* & x_1^* & \ldots & x_{N - 1}^* \\
      u_0^* & u_1^* & \ldots & u_{N - 1}^* 
    \end{array}
  \right] = [\tilde{x}_0^*, \tilde{x}_1^*, \ldots, \tilde{x}_{N - 1}^*]\\
y_i^* = [x_{i,1}^*, x_{i, 2}^*,\ldots, x_{i, N}^*]^\top
\tag{3}

ただし, $N > 0$は学習データの数。ここで, 学習データ$\mathcal{D}_i$のもと, テストデータ$x\in\mathbb{R}^{n_x}, u\in\mathbb{R}^{n_u}$に対する出力データ$\hat{f}_i(x, u)$の予測分布$p(\hat{f}_i(x, u) |x, u, \mathcal{D}_i)$はベイズの定理を用いるとガウス分布として求まる。つまり,

\hat{f}_i(x, u) \sim p(\hat{f}_i(x, u) |x, u, \mathcal{D}_i) = \mathcal{N}(\hat{\mu}_i(x, u, \mathcal{D}_i), \hat{\sigma}_i(x, u, \mathcal{D}_i))
\tag{4}

\hat{\mu}_i(x, u, \mathcal{D}_i) = \int \hat{f}_i(x, u) p(\hat{f}_i(x, u) |x, u) {\rm d}\hat{f}_i(x, u) = k_{*, i}^\top(\tilde{x})(K_i + \sigma_{\varepsilon}^2 I)^{-1} y_i^* = k_{*, i}^\top(\tilde{x})\beta_i 
\tag{5}

\hat{\sigma}_i(x, u, \mathcal{D}_i)^2 = \int (\hat{f}_i(x, u) - \mu_i(x, u))^2 p(\hat{f}_i(x, u) |x, u) {\rm d}\hat{f}_i(x, u) = k_i(\tilde{x}, \tilde{x}) - k_{*, i}^\top(\tilde{x})(K_i + \sigma_{\varepsilon}^2 I)^{-1}k_{*, i}(\tilde{x})
\tag{6}

ただし, $\tilde{x} = [x^\top, u^\top]^\top$, $k_{*, i}(\tilde{x}) = [k_i(\tilde{x}, \tilde{x}_1^{\ast}), \ldots, k_i(\tilde{x}, \tilde{x}_N^{\ast})]^\top$, $K_i$の$(p, q)-$成分は $K_{i, pq} = k_i(\tilde{x}_{p}^{\ast},\tilde{x}_{q}^{\ast})$。 $\hat{f}_i$, $i\in\mathbb{N}_{1: n_x}$全体の予測分布は異なる要素間の共分散を$0$として，出力の各要素に独立にGP回帰を適用すれば良いので

\hat{f}(x, u) \sim p(\hat{f}(x, u)|x, u, \mathcal{D}) = \mathcal{N}(\hat{\mu}(x, u, \mathcal{D}), \hat{\Sigma}(x, u, \mathcal{D}))
\tag{7}

\hat{\mu}(x, u, \mathcal{D}) = [\hat{\mu}_1(x, u), \ldots, \hat{\mu}_{n_x}(x, u)]^\top
\tag{8}

\hat{\Sigma}(x, u, \mathcal{D}) = {\rm diag}(\hat{\sigma}_1(x, u)^2, \ldots, \hat{\sigma}_{n_x}(x, u)^2)
\tag{9}

ただし, $\hat{f}(x, u) = [\hat{f}_1(x, u), \ldots, \hat{f}_{x_x}(x, u)]^\top$, $\mathcal{D} = \bigl\{\mathcal{D}_i \bigr\}_{i\in\mathbb{N}_{1:n_x}}$。また, 関数$f$の平衡点は$0 = f (0, 0) $と仮定し, システムを原点で安定化させることを制御目的とする。

自己駆動型制御器の設定

自己駆動型制御器を定義する。通信時刻を$k_i\in\mathbb{N}_{\ge 0}$, $i = 0, 1, 2, \ldots, $とおく。ただし, $k_0 = 0$, $k_{i + 1} > k_i$, $i\in \mathbb{N}_{\ge 0}$。また, $i\in\mathbb{N}_{\ge 0}$ に対応するinter-communication time step を

m_{i} = k_{i + 1} - k_i > 0 
\tag{10}

とおく。そして, 自己駆動型制御器に関する次の方策$\pi = \bigl\{\pi_{inp}, \pi_{com}\bigr\}$を設計することでプラントと制御器の通信回数を削減すること考える。

$\pi_{\rm inp}$: $\mathbb{R}^{n_x} \to \mathbb{R}^{n_u}$は制御方策。
$\pi_{\rm com}$: $\mathbb{R}^{n_x} \to \mathbb{N}_{1: M}$は通信方策。ただし, $M\in\mathbb{N}_{> 0}$ は inter-communication time stepの最大値。

自己駆動制御器の手続きは各$k_i$, $i \in\mathbb{N}_{\ge 0}$に対して次の3ステップで書くことができる。

プラントが状態$x_{k_i}$を観測し, $x_{k_i}$を制御器に送信する。
制御器は制御入力とinter-communication time stepを次のように計算する。$u_{k_i} = \pi_{\rm inp}(x_{k_i})$, $m_i = \pi_{\rm com}(x_{k_i})$。
制御器は$\bigl\{u_{k_i}, m_i\bigr\}$をプラントに送信し, 次の送信時刻$k_{i + 1} = k_i + m_i$になる直前まで一定で$u_{k_i}$を適用し続ける。つまり, $u_k = u_{k_i}$, $\forall k\in \mathbb{N}_{k_il:k_{i + 1} - 1}$

ここで, 方策$\pi$に関するコスト関数(方策コスト関数)を次のように設定する。

J^{\pi}(x_{k_i}) = \sum_{\ell = i + 1}^{\infty} \mathbb{E}_{x_{k_\ell}^{\pi}}[C_1(x_{k_\ell}) + \gamma C_2(m_{\ell})]
\tag{11}

ただし, $C_1$は状態に関するコスト関数 (状態コスト関数), $C_2$はinter-communication time stepに関するコスト関数 (通信コスト関数), $\gamma > 0$は通信コストに関する重み。
例として, 状態コスト関数$C_1$と通信コスト関数$C_2$を次のように設定する。

C_1(x_{k_\ell}) = 1 - {\rm exp} \bigl\{-\frac{1}{2}x^{\top}_{k_{\ell}}Q^{-1} x_{k_\ell}\bigr\} 
\tag{12}

C_2(m_{\ell}) = M - m_\ell
\tag{13}

ただし, $Q>0$は与えられた正定値行列。つまり, 方策コスト関数$J^{\pi}$を最小化することで最適な制御方策$\pi_{\rm inp}$と通信則$\pi_{\rm com}$を設計することができる。また, $f$は未知の関数であるがGP回帰で学習することができるので, GPモデル$\hat{f}$にを使うことで最適な方策を設計する。

状態空間の近似

プラントのGPモデルが式(7)で得られる時, 現在の時刻$k$の状態$x_k$から一定入力$u$のもとで$m \ge 2$ステップ後の状態空間$x_{k + m}$の予測分布$p(x_{k + m}| x_k, u, \mathcal{D})$を解析モーメントマッチングを用いてガウス近似する。まず, 状態$x_k$から$\ell \in \mathbb{N}_{1:m - 1}$ 時刻後まで一定値入力$u$を加えた時の状態$x_{k + \ell}$ の予測分布が次のように与えられると仮定する。　　

p(x_{k + \ell} | x_k, u, \mathcal{D})  = \mathcal{N} (\mu(x_{k + \ell - 1}, u, \mathcal{D}), \Sigma(x_{k + \ell - 1}, u, \mathcal{D}))
\tag{14}

ただし, この過程のもとでは状態$x_{k + 1}$の予測分布の平均と分散はあらかじめ以下のように定まることに注意する。　

\begin{align}
\mu(x_k, u, \mathcal{D}) &= \hat{\mu}(x_k, u, \mathcal{D})\\
\Sigma(x_k, u, \mathcal{D})&= \hat{\Sigma}(x_k, u, \mathcal{D})
\end{align}
\tag{15}

この時, $\hat{f}(x_{k + \ell}, u)$の予測分布$p(\hat{f}(x_{k + \ell}, u)| x_k, u, \mathcal{D})$は次の予測ステップにより推定される。　

\begin{align}
p(\hat{f}(x_{k + \ell}, u)| x_k, u, \mathcal{D}) &= \int p(\hat{f}(x_{k + \ell}, u)| x_{k + \ell}, u_{k + \ell}, \mathcal{D})p(x_{k + \ell}, u_{k + \ell} |x_k, u, \mathcal{D}) {\rm d}[x_{k + \ell}^\top, u_{k + \ell}^\top]^\top \\
& = \int p(\hat{f}(x_{k + \ell}, u)| x_{k + \ell}, u, \mathcal{D})p(x_{k + \ell}|x_k, u, \mathcal{D})p(u_{k + \ell} |x_k, u, \mathcal{D}) {\rm d}[x_{k + \ell}^\top, u_{k + \ell}^\top]^\top \\
& = \int p(\hat{f}(x_{k + \ell}, u)| x_{k + \ell}, u, \mathcal{D})p(x_{k + \ell}|x_k, u, \mathcal{D}){\rm Dirac}(u_{k + \ell} - u ) {\rm d}[x_{k + \ell}^\top, u_{k + \ell}^\top]^\top \\
& = \int p(\hat{f}(x_{k + \ell}, u)| x_{k + \ell}, u, \mathcal{D})p(x_{k + \ell}|x_k, u, \mathcal{D}){\rm d}x_{k + \ell}
\end{align}
\tag{16}

だだし, ${\rm Dirac}(\cdot)$はディラックデルタ関数。一般に, 式(15)の厳密な計算は難しいが, ガウスカーネルを用いたGP回帰モデルに対しては, $p(\hat{f}(x_{k + \ell}, u)| x_k, u)$の平均$\mu (x_{k + \ell}, u, \mathcal{D})$と共分散行列$\Sigma(x_{k + \ell}, u, \mathcal{D})$はガウス近似により解析的に計算ができる。まず, $p(\hat{f}(x_{k + \ell}, u)| x_k, u)$の平均$\mu (x_{k + \ell}, u, \mathcal{D})$は

\begin{align}
\mu (x_{k + \ell}, u, \mathcal{D}) &= \int \int \hat{f}(x_{k + \ell}, u) p(\hat{f}(x_{k + \ell})| x_{k + \ell}, u, \mathcal{D})p(x_{k + \ell}|x_k, u, \mathcal{D}){\rm d}\hat{f}(x_{k + \ell}, u) {\rm d}x_{k + \ell}\\ 
&=  \int \hat{\mu}(x_{k + \ell}, u, \mathcal{D})\mathcal{N} (\mu(x_{k + \ell - 1}, u, \mathcal{D}), \Sigma(x_{k + \ell - 1}, u, \mathcal{D})){\rm d}x_{k + \ell} \\
&= [\beta_i^\top \int k_{*, i}(\tilde{x}_{k + \ell})  \mathcal{N} (\mu(x_{k + \ell - 1}, u, \mathcal{D}), \Sigma(x_{k + \ell - 1}, u, \mathcal{D})){\rm d}x_{k+ \ell}]_{i\in\mathbb{N}_{1:n_x}} \\
&= [\beta_i^\top \eta_i]_{i\in\mathbb{N}_{1:n_x}}
\end{align}
\tag{17}

ただし, $\eta_i = [\eta_{i, 1}, \ldots, \eta_{i, N}]$は次のように与えられる。

\eta_{i, n} = \alpha_i^2 |\tilde{\Sigma}_{k + \ell}\Lambda_i^{-1} + I|^{-\frac{1}{2}} \times {\rm exp} (-\frac{1}{2}(\tilde{x}_n^* - \tilde{\mu}_{k + \ell})^\top(\Lambda_i + \tilde{\Sigma}_{k + \ell})^{-1}(\tilde{x}_n^* - \tilde{\mu}_{k + \ell}))
\tag{18}

ここで, $n\in\mathbb{N}_{1: N}$, $\tilde{x}_n^* = [{x}_n^{* \top}, u_{n}^{* \top}]^\top$, $\tilde{\mu}_{k + \ell} = [\mu(x_{k + \ell - 1}, u, \mathcal{D})^\top, u^\top]^\top$, $\tilde{\Sigma}_{k + \ell} = {\rm diag}(\Sigma(x_{k + \ell - 1}, u, \mathcal{D}), 0_{n_u \times n_u})$とおいた。

次に, $p(\hat{f}(x_{k + \ell}, u)| x_k, u)$の共分散行列$\Sigma(x_{k + \ell}, u, \mathcal{D})$もガウス近似により解析的に計算する。共分散行列$\Sigma(x_{k + \ell}, u, \mathcal{D})$の各成分を$\sigma_{ij}(x_{k + \ell}, u, \mathcal{D})^2$, $i, j \in \mathbb{N}_{1: n_x}$とおく。ただし,

\begin{align}
\sigma_{ij}(x_{k + \ell}, u, \mathcal{D})^2 &= \int\int \int (\hat{f}_i(x_{k + \ell}, u) - \mu_i(x_{k + \ell}, u, \mathcal{D})) (\hat{f}_j(x_{k + \ell}, u) - \mu_j(x_{k + \ell}, u, \mathcal{D})) p(\hat{f}_i(x_{k + \ell}, u), \hat{f}_j(x_{k + \ell}, u) |x_{k + \ell}, u, \mathcal{D}) {\rm d}\hat{f}_i(x_{k + \ell}, u) {\rm d}\hat{f}_j(x_{k + \ell}, u) {\rm d}x_{k + \ell} \\
&= \begin{cases}
    \mathbb{V}_{x_{k + \ell}}[\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]] + \mathbb{E}_{x_{k + \ell}}[\mathbb{V}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]] &, (i = j) \\
    \mathbb{V}_{x_{k + \ell}}[\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u), \hat{f}_j(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u), \hat{f}_j(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]] &, (i \neq j)
  \end{cases}\\
&= \begin{cases}
    \mathbb{E}_{x_{k + \ell}}[\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]^2] - \mathbb{E}_{x_{k + \ell}}[\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]]^2 + \mathbb{E}_{x_{k + \ell}}[\mathbb{V}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]] &, (i = j)\\
    \mathbb{E}_{x_{k + \ell}}[\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]\mathbb{E}_{\hat{f}_j(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]] - \mathbb{E}_{x_{k + \ell}}[\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_i(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]\mathbb{E}_{\hat{f}_i(x_{k + \ell}, u)}[\hat{f}_j(x_{k + \ell}, u)|x_{k + \ell}, u, \mathcal{D}]] &, (i \neq j)
  \end{cases}\\
&= \begin{cases}
\beta_i^\top L_{ii} \beta_i - \mu_i(x_{k + \ell}, u, \mathcal{D})^2 + \alpha_i^2 - {\rm tr}((K_i + \sigma_{\varepsilon}^2 I)^{-1} L_{ii}) + \sigma_{\varepsilon}^2&, (i = j) \\
\beta_i^\top L_{ij} \beta_j - \mu_i(x_{k + \ell}, u, \mathcal{D}) \mu_j(x_{k + \ell}, u, \mathcal{D})&, (i \neq j)\\
\end{cases}
\end{align}
\tag{19}

ここで, 式(19)の式変形について, $p(\hat{f}_i(x_{k + \ell}, u)| x_k, u, \mathcal{D})$と$p(\hat{f}_j(x_{k + \ell}, u)| x_k, u, \mathcal{D})$, $i\neq j$の要素間の共分散は$0$であることに注意する。また, $\tilde{x}_{pq, ij}^* = \Lambda_i^{-1}(\tilde{x}_{p}^* - \tilde{\mu}_{k + \ell}) + \Lambda_j^{-1}(\tilde{x}_{q}^* - \tilde{\mu}_{k + \ell})$と$R_{ij} = (\Lambda_i^{-1} + \Lambda_j^{-1})\tilde{\Sigma}_{k + \ell} + I$とすると, 行列$L_{ij}\in\mathbb{R}^{N\times N}$の$(p, q)-$成分は次のように与えられる。

L_{ij, pq} = |R_{ij}|^{-\frac{1}{2}}k_i(\tilde{x}_p^*,\tilde{\mu}_{k + \ell})k_j(\tilde{x}_q^*,\tilde{\mu}_{k + \ell}) \times {\rm exp} (\frac{1}{2}\tilde{x}_{pq, ij}^{* \top} (\tilde{\Sigma}_{k + \ell}^{-1} + \Lambda_i^{-1} + \Lambda_j^{-1})^{-1}\tilde{x}_{pq, ij}^*)
\tag{20}

よって, $\hat{f}(x_{k + \ell}, u)$の予測分布はガウス近似が可能であり,

p(\hat{f}(x_{k + \ell}, u)| x_k, u, \mathcal{D}) = \mathcal{N} (\mu (x_{k + \ell}, u, \mathcal{D}), \Sigma (x_{k + \ell}, u, \mathcal{D}))
\tag{21}

以上より, $\ell = 1,\ldots, m - 1$から順に上記の操作を繰り返すことで, $p(x_{k + m} | x_k, u) \approx p(\hat{f}(x_{k + m - 1}, u)| x_k, u, \mathcal{D})$を近似することができる。

近似価値反復

プラントのGPモデルが式(7)で得られる時, 方策コスト関数(11)を最小化する自己駆動制御器を設計する。方策コスト関数(11)の最適ベルマン方程式は次のように与えられる。

J^*(x_{k_i}) = \min_{u_{k_i}, m_i}\bigl\{\mathbb{E}_{x_{k_{i + 1}}}[C_1(x_{k_{i + 1}}) + \gamma C_2(m_{i + 1}) + J^*(x_{k_{i + 1}})] \bigr\} 
\tag{22}

ただし, 状態空間$\mathbb{R}^{n_x}$と入力空間$U$が無限であるので, 式(22)の解析解を求めることは不可能なため, 離散化すること考える。まず, $N_X$と$N_U$を状態空間と入力空間を表現する数として, 有限の状態空間と入力空間を次のように与える。

X_R = \bigl\{x_{R, 1},\ldots, x_{R, N_X}\bigr\}
\tag{23}

U_R = \bigl\{u_{R, 1},\ldots, u_{R, N_U}\bigr\}
\tag{24}

ただし, $x_{R, i}\in\mathbb{R}^{n_x}$, $i\in\mathbb{N}_{1:N_X}$と$u_{R, j}\in\mathbb{R}^{n_u}$, $j\in\mathbb{N}_{1:N_U}$。この時, 状態コスト関数を$C_1$を式(12)で, また通信コスト関数を式(13)で定義していることから, 放射基底関数(RBFs)を用いることで, 最適方策コスト関数$\hat{J}^{*}$は

\hat{J}^* (x) = \sum_{n = 1}^{N_X} w_{J, n} {\rm exp}(-\frac{1}{2\sigma_{J}^2}\|x - x_{R, n}\|^2) 
\tag{25}

最適制御方策$\hat{\pi}_{\rm inp}^* $は

\hat{\pi}_{\rm inp}^* (x) = \sum_{n = 1}^{N_X} w_{u, n} {\rm exp}(-\frac{1}{2\sigma_{u}^2}\|x - x_{R, n}\|^2) 
\tag{26}

最適通信方策$\hat{\pi}_{\rm com}^* $は

\hat{\pi}_{\rm com}^* (x) = {\rm argmin}_{j\in\mathbb{N}_{>0}}|\sum_{n = 1}^{N_X} w_{c, n} {\rm exp}(-\frac{1}{2\sigma_{c}^2}\|x - x_{R, n}\|^2) - j|
\tag{27}

で近似することができる。ただし, $\bigl\{w_{J, n}\bigr\}_{n=1}^{N_X}$, $\bigl\{w_{u, n}\bigr\}_{n=1}^{N_X}$, $\bigl\{w_{c, n}\bigr\}_{n=1}^{N_X}$は重み, $\sigma_{J}$,$\sigma_{u}$, $\sigma_{c}$はRBFsのハイパーパラメータ。

ここで, 各$x\in X_R$に対して, 全ての$u\in U_R$と$m\in\mathbb{N}_{1: M}$に対応する $D(x, u, m)$を次のように与える。

\begin{align}
D(x, u, m) &= \mathbb{E}_{x_m} [C_1(x_m) + \gamma C_2(m') + \hat{J}^*(x_m)] \\
&= \int p(x_m | x, u, \mathcal{D}) C_1 (x_m) {\rm d}x_m + \gamma\int p(x_m | x, u, \mathcal{D}) C_2 (m') {\rm d}x_m + \int p(x_m | x, u, \mathcal{D}) \hat{J}^* (x_m) {\rm d}x_m
\end{align}
\tag{28}

ただし, $x_m$は状態$x$から$m$時間ステップ制御入力$u$を適用することによって到達する状態, $m'$は状態$x_m$に対応するinter-communicaiton time stepであり, $m' = \hat{\pi}^*_{\rm com}(x_m)$。この$D(x, u, m)$を価値反復することでハイパーパラメータ$\bigl\{w_{J, n}\bigr\}_{n=1}^{N_X}$, $\bigl\{w_{u, n}\bigr\}_{n=1}^{N_X}$, $\bigl\{w_{c, n}\bigr\}_{n=1}^{N_X}$, $\sigma_{J}$,$\sigma_{u}$, $\sigma_{c}$を更新する。ここで, 式(28)について解析的に計算する。

まず, $\int p(x_m | x, u) C_1 (x_m) {\rm d}x_m$ の項について　

\begin{align}
\int p(x_m | x, u) C_1 (x_m) {\rm d}x_m & \approx \int \mathcal{N}(\mu_m, \Sigma_m)C_1(x_m)dx_m \\
&= 1 - \int \mathcal{N}(\mu_m, \Sigma_m) {\rm exp}(-\frac{1}{2} x_m^\top Q^{-1} x_m) dx_m \\
&=1 - |I + \Sigma_m Q|^{-\frac{1}{2}} {\rm exp} \bigl\{-\frac{1}{2}\mu_m (Q + \Sigma_m)^{-1} \mu_m\bigr\}
\end{align}
\tag{29}

ただし, 最初の近似については式(15)と同じ操作を行なった。　

次に, $\int p(x_m | x, u) C_2 (m') {\rm d}x_m$ の項について　

\begin{align}
\int p(x_m | x, u) C_2 (m') {\rm d}x_m &= M - \int p(x_m | x, u) \hat{\pi}_{\rm com}^* (x_m) {\rm d}x_m \\
& \approx M - \sum_{n = 1}^{N_X} w_{c, n} \int \mathcal{N}(\mu_m, \Sigma_m) {\rm exp}(-\frac{1}{2\sigma_{c}^2}\|x - x_{R, n}\|^2) {\rm d}x_m \\
&= M -  \sum_{n = 1}^{N_X} w_{c, n}\Bigl\{|I + \sigma_c^{-2}\Sigma_m|^{-\frac{1}{2}} {\rm exp} \bigl\{-\frac{1}{2\sigma_c^2}(\mu_m - x_{R, n})^\top (I + \Sigma_m \sigma_c^{-2})^{-1} (\mu_m - x_{R, n})\bigr\} \Bigr\}
\end{align}
\tag{30}

ただし, $\hat{\pi}_{\rm com}^* (x) \approx \sum_{n = 1}^{N_X} w_{c, n} {\rm exp}(-\frac{1}{2\sigma_{c}^2}|x - x_{R, n}|^2)$として近似した。

最後に, $\int p(x_m | x, u) \hat{J}^* (x_m) {\rm d}x_m$ の項について　

\begin{align}
\int p(x_m | x, u) \hat{J}^* (x_m) {\rm d}x_m &\approx \sum_{n = 1}^{N_X} w_{J, n} \int \mathcal{N}(\mu_m, \Sigma_m) {\rm exp}(-\frac{1}{2\sigma_{J}^2}\|x - x_{R, n}\|^2) {\rm d}x_m\\
&= \sum_{n = 1}^{N_X} w_{c, n}\Bigl\{|I + \sigma_J^{-2}\Sigma_m|^{-\frac{1}{2}} {\rm exp} \bigl\{-\frac{1}{2\sigma_J^2}(\mu_m - x_{R, n})^\top (I + \Sigma_m \sigma_J^{-2})^{-1} (\mu_m - x_{R, n})\bigr\} \Bigr\}
\end{align}
\tag{31}

以上の準備により, 近似価値反復アルゴリズムは以下のように与えることができる。

アルゴリズム1
1. $\hat{\pi}^{\ast}_{\rm com}$, $\hat{\pi}^{ \ast}_{\rm com}$, $\hat{\pi}^{\ast}_{\rm com}$のハイパーパラメータを初期化する。　
2. $N_{\rm ite}$回, step3からstep8の操作を繰り返す。
3. $x\in X_{R}$に対して, step4からstep7の操作を繰り返す。
4. $D_{\rm min} \leftarrow \infty$
5. $u\in U_R$, $m\in\mathbb{N}_{1: M}$に対して, step6からstep7の操作を繰り返す。
6. D(x, u, m)を計算する。
7. もし, $D(x, u, m) < D_{\rm min}$ ならば, $D_{\rm min} \leftarrow D(x, u, m)$, $u^{\ast}(x) \leftarrow u$, $m^{\ast}(x) \leftarrow m$ とする。
8. $\hat{\pi}^{\ast}_{\rm com}$, $\hat{\pi}^{ \ast}_{\rm com}$, $\hat{\pi}^{\ast}_{\rm com}$のハイパーパラメータを新しい学習データ$x_{R, n}$, $D^{\ast}(x_{R, n})$, $u^{\ast}(x)$, $m^{\ast}(x_{R, n})$で更新する。ただし, $n\in\mathbb{N}_{1:N_X}$。

強化学習に基づいたプラントと自己駆動型制御器のjoint学習

強化学習に基づいてプラントと自己駆動型制御器のダイナミクスをjoint学習する。まず, 学習エージェントは初期状態で知識を持っていないと仮定しているので, 効率的に学習するために, step1で初期化をする。次に, $\varepsilon$の確率でランダムで制御方策と通信則を取り出し, 学習データを更新する。それ以外の場合は,　学習データを望にプラントのGPモデルを更新し, アルゴリズム1に従って最適制御方策と最適通信方策を更新する。　

上記のアルゴリズムの概要を以下にに示す。

アルゴリズム2
1. $\hat{\pi}^{\ast}_{\rm inp}$を初期化し, $x\in\mathbb{R}^{n_x}$に対して, $\hat{\pi}^{\ast}_{\rm com} \leftarrow 1$とする。　
2. $i\in\mathbb{N}_{1:n_x}$に対して, $X \leftarrow \bigl\{\bigr\}$, $y_i \leftarrow \bigl\{\bigr\}$, $\mathcal{D}_i \leftarrow \bigl\{X, y_i\bigr\}$で学習データを初期化する。
3. $N_{\rm epi}$回, step4からstep11のエピソードを繰り返す。
4. $\ell \leftarrow 0$, $k_{\ell} \leftarrow 0$, $x_{k_{\ell}} \leftarrow x_{\rm init}$とし, プラントに$x_{k_{\ell}}$を送信する。　
5. $\ell \in \mathbb{N}_{0: N_{\rm max} - 1}$に対して, step6からstep10の操作を行う。
6. $r \sim {\rm Uniform}[0, 1]$から取り出す。
7. もし, $r < \varepsilon$ ならば, $m_{\ell} \leftarrow 1$として, $U$の中からランダムで$u_{k_\ell}$を取り出す。それ以外の場合, $u_{k_{\ell}} \leftarrow \hat{\pi}^{\ast}_{\rm inp} (x_{k_\ell})$と $m_{k_{\ell}} \leftarrow \hat{\pi}^{\ast}_{\rm com} (x_{k_\ell})$ を更新する。　　
8. $k_{\ell + 1} \leftarrow k_{\ell} + m_{\ell}$で通信時刻を更新し, プラントに制御方策$u_{k_\ell}$とinter-communication time step $m_\ell$ を送信する。
9. プラントは$m_\ell$ステップ制御方策$u_{k_\ell}$を適用した後, 制御器に状態$x_{k_{\ell} + 1}$を送信する。
10. もし, $m_{\ell} == 1$ならば, $i\in\mathbb{N}_{1:n_x}$に対して, $X \leftarrow \bigl\{X \cup [x_{k_\ell}^{\top}, u_{k_\ell}^{\top}]^\top \bigr\}$, $y_i \leftarrow \bigl\{y_i \cup x_{k_{\ell + 1}, i} \bigr\}$, $\mathcal{D}_i \leftarrow \Bigl\{\mathcal{D}_i \cup \bigl\{X, y_i\bigr\}\Bigr\}$で学習データを更新し, 新しい学習データ$\mathcal{D} = \bigl\{\mathcal{D}_i \bigr\}_{i ~ 1}^{n_x}$を用いて, 学習エージェントはプラントのGPモデルを学習する。
11. アルゴリズム1を実行し, 最適制御方策$\hat{\pi}^{\ast}$と最適通信方策$\hat{\pi}^{\ast}_{\rm com}$を更新する。　

アルゴリズム2の主な改善点

学習データが$\varepsilon$の確率でしか更新されない。
安定・安全性が保証されていない。
価値反復に勾配法を適用する。

参考文献

[1] : K. Hashimoto, Y. Yoshimura, and T.Ushio, (2020). Learning self-triggered controllers with Gaussian processes. IEEE Transactions on Cybernetics.
[2] : K. G. Vamvoudakis and H. Ferraz, “Model-free event-triggered control algorithm for continuous-time linear systems with optimal performance,” Automatica, vol. 87, pp. 412–420, Jan. 2018.
[3] : X. Zhong, Z. Ni, H. He, X. Xu, and D. Zhao, “Event-triggered reinforce- ment learning approach for unknown nonlinear continuous-time system,” in Proc. Int. Joint Conf. Neural Netw., 2014, pp. 3677–3684
[4] : X. Yang and H. He, “Adaptive critic designs for event-triggered robust control of nonlinear systems with unknown dynamics,” IEEE Trans. Cybern., vol. 49, no. 6, pp. 2255–2267, Jun. 2019
[5] : Y. Yang, K. G. Vamvoudakis,H. Ferraz, and H. Modares, “Dynamic intermittent Q-learning-based model-free suboptimal co-design of L2- stabilization,” Int. J. Robust Nonlinear Control, vol. 29, no. 9, pp. 2673–2694, 2019
[6] : Y. Yang, K. G. Vamvoudakis, H. Ferraz, and H. Modares, “Dynamic intermittent Q-learning for systems with reduced bandwidth,” in Proc. IEEE Conf. Decis. Control (IEEE CDC), 2018, pp. 924–931.
[7] : Y. Yang, H. Modares, K. G. Vamvoudakis, Y. Yin, and D. C. Wunsch, “Dynamic intermittent feedback design for H∞ containment control on a directed graph,” IEEE Trans. Cybern., early access.
[8] : M. P. Deisenroth, D. Fox, and C.E. Rasmussen (2013). Gaussian processes for data-efficient learning in robotics and control. IEEE transactions on pattern analysis and machine intelligence, 37(2), 408-423.

未知のプラントに対する自己駆動制御器の学習