計算機の飛躍的な進歩により, 非線形システムに対する最適制御問題を数値的に解くのに必要な時間が短くなってきている. そこで, 最適制御問題をサンプリング時間内で解いて得られたフィードバック制御を用いるモデル予測制御を行う環境が現実味を帯びてきた. しかし, 一般にモデル予測制御の安定性を議論することは難しく, さらに対象とするシステムが不確かな場合は尚更難しい. そこで, この記事では不確かな非線形システムをGP回帰によって学習し, 得られたGPモデルをもとに安定性を保証したモデル予測制御を設計する一つの方法について書くことにした.

プラントのダイナミクス

制御対象である未知プラントのダイナミクス$\mathbf{f}$を次のように与える。　

\begin{align}
\mathbf{x}(t_{k+1}) &= \mathbf{f}(\mathbf{x}(t_{k}), \mathbf{u}(t_{k}))
\end{align}
\tag{1}

ただし, $k \in \mathbb{N}_{\ge 0}$, $\mathbf{x}(t_k) \in \mathcal{X} \subset \mathbb{R}^{n}$は状態, $u(t_k)\in U\subset\mathbb{R}^m$は入力. そして, 状態$\mathbf{x}(t_k)$は有界なノイズ$w_i$, $|{w}_i| \le \sigma_{w_i}^2$を含んで常に観測できるものとする. ここで, 表記の簡単のため, $\tilde{\mathbf{x}} = [\mathbf{x}^\top, \mathbf{u}^\top]^\top \in \mathcal{X} \cup \mathcal{U} \in \mathbb{R}^{n + m}$とおく. また, 関数$f_i \in \mathbb{R}$が平均0かつガウスカーネル$k_{f_i}(\cdot, \cdot)$のGPによって特徴付けられるものとする. ただし, 　

k_{f_i}(\tilde{\mathbf{x}}, \tilde{\mathbf{x}}') = \sigma_{f_i}^2\exp \left(-\frac{1}{2}(\tilde{\mathbf{x}} - \tilde{\mathbf{x}}')^\top\Lambda_{f_i}^{-1}(\tilde{\mathbf{x}} - \tilde{\mathbf{x}}')\right)
\tag{2}

ここで $\Lambda_{f_i} = {\rm diag}(\ell_{f_i, 1}^2, \ldots, \ell_{f_i, n + m}^2)$はスケーリングに関する対角行列, $\sigma_{f_i}^2$は潜在関数$f_i$の分散. また, ハイパーパラメータ$\theta_{f_i} = \left\{ \Lambda_{f_i}, \sigma_{f_i}^2 \right\}$, $\theta_\mathbf{f} = \left\{\theta_{f_i}\right\}_i$とおく. そして, 学習用データ$\mathcal{D}_{i, N_k} = \left\{\tilde{X}^\ast,\mathbf{y}_i^{\ast} \right\}$が$\tilde{X}^* = [\mathbf{\tilde{x}}_0^\ast, \ldots, \mathbf{\tilde{x}}_{N_k - 1}^\ast]^\top\in\mathbb{R}^{(n + m)\times N_k}$, $\mathbf{y}_i^{\ast} = [x_{i, 1}^\ast, \ldots, x_{i, N_k}^\ast]^\top\in\mathbb{R}^{N_k}$で与えられたとき, 入力$\mathbf{\tilde{x}}(t_k)$ に対する関数値$f_i(\mathbf{\tilde{x}}(t_k))$の事後分布の平均$\mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k})$, 分散$\sigma_{f_i}^2(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k})$と入力$\mathbf{\tilde{x}}(t_k), \mathbf{\tilde{x}}(t_k)'$に対する関数値$f_i(\mathbf{\tilde{x}}(t_k)), f_i(\mathbf{\tilde{x}}(t_k)')$共分散$\bar{k}_{f_i}(\mathbf{\tilde{x}}(t_k), \mathbf{\tilde{x}}(t_k)')$は次のように求まる. 　

\mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) = \mathbf{k}_{f_i}^{*\top}(\mathbf{\tilde{x}}(t_k))(K_{g_i}^* + \sigma_{w_i}^2I)^{-1}\mathbf{y}_i^* = \mathbf{k}_{g_i}^{*\top}(\mathbf{\tilde{x}}(t_k))\boldsymbol{\beta}_{g_i}
\tag{3}

\bar{k}_{f_i}(\mathbf{\tilde{x}}(t_k), \mathbf{\tilde{x}}(t_k)' | \mathcal{D}_{i, N_k}) = k_{f_i}(\mathbf{\tilde{x}}(t_k), \mathbf{\tilde{x}}(t_k)') - \mathbf{k}_{f_i}^{*\top}(\mathbf{\tilde{x}}(t_k))(K_{f_i}^* + \sigma_{w_i}^2 I)^{-1}\mathbf{k}_{f_i}^*(\mathbf{\tilde{x}}(t_k)')
\tag{4}

\sigma_{f_i}^2(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) = \bar{k}_{f_i}(\mathbf{\tilde{x}}(t_k), \mathbf{\tilde{x}}(t_k) | \mathcal{D}_{i, N_k})
\tag{5}

ただし, $\mathbf{k}_{f_i}^\ast(\mathbf{\tilde{x}}(t_k)) = [k_{f_i}(\tilde{\mathbf{x}}_1^\ast, \mathbf{\tilde{x}}(t_k)),\ldots, k_{f_{i}}(\tilde{\mathbf{x}}_{N_k}^\ast, \mathbf{\tilde{x}}(t_k))]^\top\in\mathbb{R}^{N_k}$, $K_{f_i}^* = [k_{f_i}(\tilde{\mathbf{x}}_p^\ast, \tilde{\mathbf{x}}_q^\ast)]_{p, q}\in\mathbb{R}^{N_k\times N_k}$, $\boldsymbol{\beta}_{f_i} = (K_{f_i}^\ast + \sigma_{w_i}^2I)^{-1}\mathbf{y}_i^\ast \in\mathbb{R}^{N_k}$, ${\boldsymbol \mu}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k| \mathcal{D}_{N_k})) = [\mu_{f_1} (\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{1, N_k}), \ldots, \mu_{f_n}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{n, N_k})]^\top \in \mathbb{R}^n$とおく.

以上の準備により, 未知のプラントのダイナミクス$\mathbf{f}$の近似モデリング$\hat{\mathbf{f}}$を$N_k$回目でGPにより得られた平均として次のように与える. 　

\hat{{\mathbf{x}}}(t_{k + 1})= \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}(t_k)| \mathcal{D}_{N_k}) = {\boldsymbol \mu}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k})
\tag{6}

また, ハイパーパラメータ$\theta_{f_i}$は学習用データ$\mathcal{D}_{i, N_k}$が与えられた時, 次の対数尤度関数を最大化することで与える.

\log p(\mathbf{y}_i^* | \tilde{X}^*, \theta_{f_i}) = -\frac{1}{2}\mathbf{y}_i^{* \top} K_{f_i}^{* -1}\mathbf{y}_i^* -\frac{1}{2} \log|K_{f_i}^{*}| - \frac{n}{2}\log(2\pi)
\tag{7}

モデル予測制御

未知のプラントのダイナミクス$\mathbf{f}$の近似モデリング$\hat{\mathbf{f}}$である式(6)を用いて, 各時刻$t_k$ごとに次の最適制御問題を考える. 　

\begin{align}
\min_{\hat{\mathbf{u}}(t)_{[k, k + N_H - 1]}} & \ V_{N_H} ({\mathbf{x}}(t_k), \hat{\mathbf{u}}(t)_{[k, k + N_H - 1]}) \notag \\ 
{\rm s.t.} & \ \hat{{\mathbf{x}}}(t_{k + i + 1})= \hat{\mathbf{f}}(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i}) | \mathcal{D}_{N_k}) \notag \\
& \ \hat{\mathbf{x}}(t_k) = {\mathbf{x}}(t_k) \notag \\
& \ \hat{\mathbf{x}}(t_{k + i}) \in \mathcal{X} \subset \mathbb{R}^n \notag \\
& \ \hat{\mathbf{u}}(t_{k + i}) \in \mathcal{U} \subset \mathbb{R}^m \notag \\
& \ i = 0, \ldots, N_H - 1
\end{align}
\tag{8}

ただし, ${\mathbf{x}}(t_k)$は初期状態, $\hat{\mathbf{u}}(t)_{[k, k + N_H - 1]} = \left\{\hat{\mathbf{u}}(t_k), \ldots, \hat{\mathbf{u}}(t_{k + N_H - 1}) \right\}$は近似モデル$\hat{\mathbf{f}}$に対する最適制御の列, $\mathcal{D}_{N_k} = \left\{\mathcal{D}_{i, N_k}\right\}_{i=1,\ldots, n}$. そして, 最適制御問題(8)のコスト関数は次のように与える. 　　

 V_{N_H} ({\mathbf{x}}(t_k), \hat{\mathbf{u}}(t)_{[k, k + N_H - 1]}) = \sum_{i = 0}^{N_H - 1} L(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i})) + \lambda V_{f}(\hat{\mathbf{x}}(t_{k + N_H}))
\tag{9}

ただし, $\lambda \ge 1$. ここで, 次の2つの仮定を用意する. 　

仮定1: $V_f(\mathbf{x_f})$が制御リアプノフ関数(CLF)であり, $\mathbf{x}_f \in \mathcal{X}_f = \left\{\mathbf{x}_f\in\mathbb{R}^n : V_f(\mathbf{x}_f) \le \nu\right\} \subseteq \mathcal{X}$, $\nu > 0$に対して,

\left\{
    \begin{array}{c}
      \alpha_1(\|\mathbf{x}_f\|) \le V_f(\mathbf{x}_f) \le \alpha_2(\|\mathbf{x}_f\|) \\
      V_f(\mathbf{x}_f^+) - V_f(\mathbf{x}_f) \le - L(\mathbf{x}_f, \kappa_f(\mathbf{x}_f)) \\
    \mathbf{x}_f^+ = \hat{\mathbf{f}}(\mathbf{x}_f, \kappa(\mathbf{x}_f)) 
    \end{array}
\right.
\tag{10}

を満たす正定関数$L(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i}))$と状態フィードバック$\kappa_f(\mathbf{x}_f) \in \mathcal{U}$が存在する.

仮定２: $\forall \hat{\mathbf{x}}(t_{k + i}) \notin \mathcal{X}_f$, $\forall \hat{\mathbf{u}}(t_{k + i})\in \mathcal{U}$に対して, $L(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i})) > d >0$を満たす$d > 0$が存在する.

そして, 学習用データ$\mathcal{D}_{N_k}$のもと得られた最適制御の列$\hat{\mathbf{u}}^\star(t)_{[k, k + N_H - 1]} $による定まる値関数を$V_{N_H}^\star ({\mathbf{x}}(t_k)) = V_{N_H} ({\mathbf{x}}(t_k), \hat{\mathbf{u}}^\star(t)_{[k, k + N_H - 1]}| \mathcal{D}_{N_k})$とし, 時刻$t_k$の最適制御を$\mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) = \hat{\mathbf{u}}^\star(t_k)$とおく. 　　　　　　

Nonimal Stability

未知のプラントのダイナミクス$\mathbf{f}$がノミナルであり, GPによりモデリングした$\hat{\mathbf{f}}$と等しい時, そのDomain of Attractionと安定性について議論する. まず, [4]の補題1より, 最適制御問題(8)について, $\hat{\mathbf{x}}^\star(t_{k + N_H}) \notin \mathcal{X}_f$ならば, $\hat{\mathbf{x}}^\star(t_{k + i}) \notin \mathcal{X}_f$, $i = 0, \ldots, N_H - 1$であることが示されている (証明略). よって, $\hat{\mathbf{x}}(t_{k + N_H}) \notin \mathcal{X}_f$ならば, 仮定1, 2より, 値関数$V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$は　

 V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) = \sum_{i = 0}^{N_H - 1} L(\hat{\mathbf{x}}^\star(t_{k + i}), \hat{\mathbf{u}}^\star(t_{k + i})) + \lambda V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H})) > N_H d + \lambda \nu
\tag{11}

が成り立つ. この対偶をとると, 　　　　　

 V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \le N_H d + \lambda \nu
\tag{12}

が成り立つならば, $\hat{\mathbf{x}}(t_{k + N_H}) \in \mathcal{X}_f$が成り立つ. つまり, 領域$\mathcal{X}_{N_H}^k(\lambda) = \left\{\mathbf{x}(t_k)\in \mathcal{X} | V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_{k}}) \le N_H d + \lambda \nu \right\}$がDomain of Attraction. そして, [6]より, 短い評価区間の最適制御は長い評価区間に対する許容制御入力になりうるので, $V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \le V_{N_H - 1}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$が成り立つことから, $\mathbf{x}(t_{k + 1}) = \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k})| \mathcal{D}_{N_k}))$とおくと　　

V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_k}) \le V_{N_H - 1}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_k}) = V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_k}) - L(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}(\mathbf{x}(t_k)| \mathcal{D}_{N_k}))
\tag{13}

が成り立つ. よって, $\forall \mathbf{x}(t_k) \in \mathcal{X}_{N_H}^k(\lambda)$に対して, 仮定1より$V_{N_H}^\star (\cdot| \mathcal{D}_{N_k})$は$\mathbf{x}(t_{k + 1}) = \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k})| \mathcal{D}_{N_k}))$のリアプノフ関数であるので原点は漸近安定である.　　

ここで, $V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k + 1}}) \le V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_k})$を満たすように学習用データ$\mathcal{D}_{N_k}$を更新させると, 式(13)と合わせて

\begin{align}
V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_{k + 1}}) &\le V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k + 1}}) - L(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}(\mathbf{x}(t_k)|\mathcal{D}_{N_{k + 1}})) \notag \\
&\le V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}}) - L(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}(\mathbf{x}(t_k)|\mathcal{D}_{N_k}))
\end{align}
\tag{14}

が成り立つ. よって, 式(12)と合わせて

N_H d + \lambda \nu  \ge V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \ge V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_{k + 1}}) \ge \ldots, 
\tag{15}

であることから, $\mathbf{x}(t_{k}) \in \mathcal{X}_{N_H}^k(\lambda)$ならば, 再帰的に$\mathbf{x}(t_{k + i}) \in \mathcal{X}_{N_H}^{k + i}(\lambda)$, $i= 0, 1, \ldots$に対しても, 最適制御問題(8)は実行可能である. 以上より, $\forall \mathbf{x}(t_k) \in \mathcal{X}_{N_H}^k(\lambda)$に対して, 仮定1より$V_{N_H}^\star (\cdot| \cdot)$は$\mathbf{x}(t_{k + 1}) = \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k})| \mathcal{D}_{N_k}))$のリアプノフ関数であるので原点は漸近安定である. 　　

最後に, $\lambda \le \lambda'$の時, $\mathbf{x}(t_{k}) \in \mathcal{X}_{N_H}^k(\lambda)$に対して, 　

\begin{align} 
 V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) &= \sum_{i = 0}^{N_H - 1} L(\hat{\mathbf{x}}^\star(t_{k + i}), \hat{\mathbf{u}}^\star(t_{k + i})) + \lambda V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H})) \notag \\ 
&= \sum_{i = 0}^{N_H - 1} L(\hat{\mathbf{x}}^\star(t_{k + i}), \hat{\mathbf{u}}^\star(t_{k + i})) + \lambda' V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H})) - (\lambda' - \lambda) V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H})) \notag \\
&\ge V_{N_H}^{\star '} ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) - (\lambda' - \lambda) V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H}))
\end{align}
\tag{16}

と変形される($V_{N_H}^{\star '} ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \le \sum_{i = 0}^{N_H - 1} L(\hat{\mathbf{x}}^\star(t_{k + i}), \hat{\mathbf{u}}^\star(t_{k + i})) + \lambda' V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H}))$とおいた). よって, $\mathbf{x}(t_{k}) \in \mathcal{X}_{N_H}^k(\lambda)$に対して, 　

V_{N_H}^{\star '} ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \le V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) + (\lambda' - \lambda) V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H})) \le N_H d + \lambda' \nu
\tag{17}

であり, $\mathbf{x}(t_{k}) \in \mathcal{X}_{N_H}^k(\lambda’)$が成り立つ. よって, $\lambda\le \lambda'$ならば, $\mathcal{X}_{N_H}^k(\lambda) \subseteq \mathcal{X}_{N_H}^k(\lambda’)$であることから, 領域$\mathcal{X}_{N_H}^k(\lambda)$は$\lambda$に比例する関数であることも確認される.

以上の議論により, 以下の定理が導かれる. 　　

定理1(Nominal Case): 仮定1のもと, $\mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k} )$が最適制御問題(8)のモデル予測制御であるとする. そして, $V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_{k + 1}}) \le V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_{k}})$で学習用データ$\mathcal{D}_{N_{k}}$を更新させる時, $\forall \mathbf{x}(t_k) \in \mathcal{X}_{N_H}^{k}(\lambda) = \left\{\mathbf{x}(t_k)\in \mathcal{X} | V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_{k}}) \le N_H d + \lambda \nu \right\} \subseteq \mathcal{X}$に対して, ${\mathbf{x}}(t_{k + 1}) = \hat{\mathbf{f}}({\mathbf{x}}(t_{k}), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k} ) | \mathcal{D}_{N_k})$は再帰的に解$\mathbf{x}(t_k)$が存在し, 原点が漸近安定である. 加えて, $\mathcal{X}_{N_H}^{k}(\lambda)$は$\lambda$に比例して大きくなる. 　　

RKHSによる状態空間の解析

上記の安定性の議論はプラントのダイナミクスがノミナルな場合で考えているが, 実際の状態空間をRKHSにより解析する. 学習用データ$\mathcal{D}_{i, N_k}$によって生成されるカーネル$\bar{k}_{f_i}$によって導入されるヒルベルト空間を$H_{\bar{k}_{f_i}}$と書くことにする. そして, 次の仮定を用意する.

仮定3: $f_i \in H_{\bar{k}_{f_i}} (= H_{k_{f_i}})$に対して$\langle f_i , \bar{k}_{f_i}(\cdot, \tilde{\mathbf{x}}(t_k)) \rangle_{\bar{k}_{f_i}} = f_i(\tilde{\mathbf{x}}(t_k))$を満たし, かつ, $||f_i||_{k_{f_i}}^2 = \langle {f}_i , {f}_i \rangle_{k_{f_i}} \le b_{f_i} < \infty$を満たす$b\_{f\_i} > 0$が存在する. 　

ただし, $|f_i|_{{k}_{f_i}} \le |f_i|_{\bar{k}_{f_i}}$に注意する. このとき, 　

\begin{align}
|f_i(\tilde{\mathbf{x}}(t_k)) - \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k})| &= |\langle f_i - \mu_{f_i}, \bar{k}_{f_i}(\cdot, \tilde{\mathbf{x}}(t_k))\rangle_{\bar{k}_{f_i}}| \le \|\bar{k}_{f_i}(\cdot, \tilde{\mathbf{x}}(t_k))\|_{\bar{k}_{f_i}} \cdot \| f_i - 
 \mu_{f_i}\|_{\bar{k}_{f_i}} \notag \\
&= \sqrt{\bar{k}_{f_i}(\tilde{\mathbf{x}}(t_k), \tilde{\mathbf{x}}(t_k) | \mathcal{D}_{i, N_k})}\cdot \|f_i - \mu_{f_i}\|_{\bar{k}_{f_i}} = \sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) \| f_i - {\mu}_{f_i}\|_{\bar{k}_{f_i}}
\end{align}
\tag{18}

が成り立つ. ここで, $|| f_i - {\mu}_{f_i}||_{\bar{k}_{f_i}}$について

\| f_i - {\mu}_{g_i}\|_{\bar{k}_{g_i}}^2 = \| f_i - {\mu}_{f_i}\|_{k_{f_i}}^2 + \sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} \{f_i(\tilde{\mathbf{x}}_p^*) - \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, p})\}^2  
\tag{19}

が成り立つことに注意する. ただし, $| f_i - {\mu}_{f_i}|_{k_{f_i}}^2$は

\begin{align}
\|f_i - {\mu}_{f_i}\|_{k_{f_i}}^2 &= \|f_i\|_{k_{f_i}}^2 - 2 \langle f_i, {\mu}_{f_i} \rangle_{k_{f_i}} + \|{\mu}_{f_i}\|_{k_{f_i}}^2 \\ 
& = \|f_i\|_{k_{f_i}}^2 - 2 \langle f_i, \mathbf{k}_{f_i}^{*\top}(\cdot)\boldsymbol{\beta}_{f_i} \rangle_{k_{f_i}} + \|\mathbf{k}_{g_i}^{*\top}(\cdot)\boldsymbol{\beta}_{f_i} \|_{k_{f_i}} \notag \\ 
&= \|f_i\|_{k_{f_i}}^2 - 2 {\boldsymbol \beta}_{f_i}^\top \mathbf{f}_i + {\boldsymbol \beta}_{f_i}^\top \mathbf{y}_i^* - \sigma_{w_i}^2 \| {\boldsymbol \beta}_{f_i}\|_{k_{f_i}}^2
\end{align}
\tag{20}

であり, $\sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} \left\{f_i(\tilde{\mathbf{x}}_p^*) - \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, p})\right\}^2$は

\begin{align}
\sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} \{f_i(\tilde{\mathbf{x}}_p^*) - \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, p})\}^2 &= \sigma_{w_i}^{-2} \sum_{p = 0}^{N_k - 1} \{\sigma_{w_i}^2 {\beta}_{g_i, p} - w_{i, p} \}^2 \notag \\
&= \sigma_{w_i}^2 \|{\boldsymbol \beta}_{f_i} \|_{k_{f_i}} - 2 \sum_{p = 0}^{N_k - 1} {\beta}_{f_i, p} w_{f_i, p} + \sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} w_{f_i, p}^2
\end{align}
\tag{21}

である. よって, $|| f_i - {\mu}_{f_i}||_{\bar{k}_{f_i}}$は

\begin{align}
\| f_i - {\mu}_{f_i}\|_{\bar{k}_{g_i}}^2 &= \|f_i\|_{k_{f_i}}^2 - 2 {\boldsymbol \beta}_{f_i}^\top \mathbf{f}_i + {\boldsymbol \beta}_{f_i}^\top \mathbf{y}_i^* - 2 \sum_{p = 0}^{N_k - 1} {\beta}_{f_i, p} w_{f_i, p} + \sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} w_{f_i, p}^2 \notag \\
&= \|f_i\|_{k_{f_i}}^2 - {\boldsymbol \beta}_{f_i}^\top \mathbf{y}_i^* + \sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} w_{f_i, p}^2
\end{align}
\tag{22}

と計算される. 以上より,

\begin{align}
|f_i(\tilde{\mathbf{x}}(t_k)) - \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k})| &= \sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) \sqrt{\left(\|f_i\|_{k_{f_i}}^2 - {\boldsymbol \beta}_{f_i}^\top \mathbf{y}_i^* + \sigma_{w_i}^{-2}\sum_{p = 0}^{N_k - 1} w_{f_i, p}^2  \right)} \notag \\
&\le \sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) \sqrt{ b_i^2 - {\boldsymbol \beta}_{f_i}^\top \mathbf{y}_i^* + N_k }
\end{align}
\tag{23}

ここで, $q_{f_i, N_k} = \sqrt{b_i^2 - {\boldsymbol \beta}_{f_i}^\top \mathbf{y}_i^* + N_k}$とおくと

f_i(\tilde{\mathbf{x}}(t_k)) \in [ \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) - q_{f_i, N_k}\sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}), \mu_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k}) + q_{f_i, N_k}\sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k})]
\tag{24}

であることから, $|\Delta{q_{f_i, N_k}} | \le |q_{f_i, N_k}|$, $\Delta{Q_{\mathbf{f}, N_k}} = {\rm diag}(\Delta{q_{f_1, N_k}}, \ldots, \Delta{q_{f_n, N_k}})$, $Q_{\mathbf{f}, N_k} = {\rm diag}(|q_{f_1, N_k}|, \ldots, |q_{f_n, N_k}|)$, ${\boldsymbol \sigma}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) = [\sigma_{f_1}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{1, N_k}), \ldots, \sigma_{f_n}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{n, N_k})]^\top$, $\Delta{\boldsymbol \sigma}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) = \Delta{Q_{\mathbf{f}, N_k}} {\boldsymbol \sigma}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$とすると

\begin{align}
\mathbf{f}(\mathbf{x}(t_k), \mathbf{u}(t_k)) &= {\boldsymbol \mu}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) + \Delta{Q_{\mathbf{f}, N_k}} {\boldsymbol \sigma}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \notag\\
&= \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}(t_k)| \mathcal{D}_{N_k}) + \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{N_k})\notag \\
&:= F(\mathbf{x}(t_k), \mathbf{u}(t_k), \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})
\end{align}
\tag{25}

と得られる. 最後に, 以下での安定性の議論のため, $f_i$が有界であることから, モデルの不確かさ$|\sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{i, N_k})| \le \xi_{i, N_k} < \infty$, $\forall i = 1,\ldots, n$を満たす${\boldsymbol \xi}\_{N\_k} = [\xi\_{1, N\_k}, \ldots, \xi\_{n, N\_k}]^\top > 0$を定義しておく.　

Robust Stability

定理1をもとに未知のプラントがモデルの不確かさ$\sigma_{f_i}(\tilde{\mathbf{x}}(t_k)| \mathcal{D}_{n, N_k})$に対して, 定理1で得られたモデル予測制御$\mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k})$がInput to State Stability (ISS)を満たすことを示す. (ISSの定義については省略する.) 　

まず, 定理1より, $V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$は$\forall \mathbf{x}(t_k)\in\mathcal{X}_{N_H}^{k}(\lambda)$で$\mathbf{x}(t_{k + 1}) = \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k})| \mathcal{D}_{N_k}))$のリアプノフ関数であるので,

\begin{align}
\alpha_1(\|\mathbf{x}(t_k)\|) \le V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \le \alpha_2(\|\mathbf{x}(t_{k})\|)\notag \\
V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_{k + 1}}) - V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) \le - \alpha_3(\|\mathbf{x}(t_k)\|)
\end{align}
\tag{26}

を満たす$\mathcal{K}_\infty$関数$\alpha_1(\cdot), \alpha_2(\cdot), \alpha_3(\cdot)$が存在する. よって, $\alpha_2^{-1}(V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})) \le \left\|\mathbf{x}(t_k)\right\|$であることと$\mathcal{K}$関数$\gamma(\cdot)$を用いて　

\begin{align}
V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_{k + 1}}) &\le V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) - \alpha_3(\|\mathbf{x}(t_k)\|) + \gamma(\Delta{\boldsymbol \sigma}_{\mathbf{f}}(t_k)) \notag\\ 
&\le V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) - \alpha_3 \circ \alpha_2^{-1} \circ V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) + \gamma(\|Q_{\mathbf{f}, N_k}{\boldsymbol \xi}_{N_k}\|) \notag\\
&= (I_{\rm d} - \alpha_3 \circ \alpha_2^{-1}) \circ V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}) + \gamma(\|Q_{\mathbf{f}, N_k}{\boldsymbol \xi}_{N_k}\|)
\end{align}
\tag{27}

が成り立つ. ただし, $I_{\rm d}$は恒等関数. ここで, [7]の補題B.1より, $\forall \eta > 0$に対して, $\gamma(\left\|Q_{\mathbf{f}, N_k}{\boldsymbol \xi}_{N_k}\right\|)\le \gamma'(\eta)$を満たす$\mathcal{K}$関数$\gamma' = \alpha_3 \circ \alpha_2^{-1}$が存在する. よって, 領域$\Omega_{\eta}^{k}(\lambda)$を

\Omega_\eta^{k}(\lambda) = \{\mathbf{x}(t_k)\in \mathcal{X} | V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_{k}}) \le \eta \} \subseteq \mathcal{X}_{N_H}^{k}(\lambda)
\tag{28}

と定めると, $\forall \mathbf{x}(t_k) \in \Omega_{\eta}^{k}(\lambda)$に対して, 　

 \begin{align}
V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_{k + 1}}) &\le (I_{\rm d} - \alpha_3 \circ \alpha_2^{-1})(\eta) + \gamma'(\eta) = \eta
\end{align}
\tag{29}

が成り立つ. ただし, $0 < \eta \le N_H d + \lambda \nu$. よって,$\forall \mathbf{x}(t_k) \in \Omega_\eta^{k}(\lambda)$に対して, $\mathbf{x}(t_{k + 1}) = \hat{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)|\mathcal{D}_{N_k})| \mathcal{D}_{N_k}))$は実行可能であり, 領域$\Omega_\eta^{k}(\lambda)$は未知のプラント$\mathbf{f}$のinvariant set. 　　　　

ここで, 次の仮定を用意する. 　

仮定4: ノミナルモデル$\hat{\mathbf{x}}(t_{k + i + 1}) = \hat{\mathbf{f}}(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i})|\mathcal{D}_{N_k} )$が$i = 0, \ldots, N_H - 1$で, $\forall \hat{\mathbf{x}}(t_{k + i}) \in \mathcal{X}_{N_H}^{k + i}(\lambda)$, $\forall \hat{\mathbf{u}}(t_{k + i}) \in\mathcal{U}$と全ての$\mathcal{D}_{N_k}$において, $\mathbf{x}(t_k)$に関して一様連続. また, ステージコスト$L(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i}))$, $V_{f}(\hat{\mathbf{x}}(t_{k + N_H}))$が$\forall \hat{\mathbf{x}}(t_{k + i}) \in \mathcal{X}_{N_H}^{k + i}(\lambda)$, $\forall \hat{\mathbf{u}}(t_{k + i}) \in\mathcal{U}$, $i = 0, \ldots, N_H - 1$において, $\mathbf{x}(t_k)$に関して一様連続.

仮定4のもと, 学習用データ$\mathcal{D}_k$が与えられた時, $\forall \mathbf{x}(t_k) \in \Omega_\eta^{k}(\lambda)$において, $V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}})$が$\mathbf{x}(t_k)$に関して一様連続であることを示し, その性質を用いることで, $V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}})$がモデルの不確かさに関して未知のプラント$\mathbf{x}(t_{k+1}) = F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}), \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})$のISSリアプノフ関数となることを示す.

まず, 仮定4の各種関数の一様連続条件から, $\left\|\hat{\mathbf{f}}(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i})|\mathcal{D}_{N_k}) - \hat{\mathbf{f}}(\hat{\mathbf{z}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i})|\mathcal{D}_{N_k})\right\| \le \rho_x( \left\| \hat{\mathbf{x}}(t_{k + i}) - \hat{\mathbf{z}}(t_{k + i})\right\|)$, $\left\|L(\hat{\mathbf{x}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i})) - L(\hat{\mathbf{z}}(t_{k + i}), \hat{\mathbf{u}}(t_{k + i}))\right\| \le \rho_\ell( \left\| \hat{\mathbf{x}}(t_{k + i}) - \hat{\mathbf{z}}(t_{k + i})\right\|)$, $\left\|V_{f}(\hat{\mathbf{x}}(t_{k + N_H})) - V_{f}(\hat{\mathbf{z}}(t_{k + N_H}))\right\| \le \rho_f( \left\| \hat{\mathbf{x}}(t_{k + i}) - \hat{\mathbf{z}}(t_{k + i})\right\|)$を満たす$\mathcal{K}$関数$\rho_x(\cdot)$, $\rho_\ell(\cdot)$, $\rho_f(\cdot)$が存在する. ここで, [8]の補題2より, $\rho_x(\cdot)$, $i= 0, \ldots, N_H - 1$に関して, $\left\| \hat{\mathbf{x}}(t_{k + i}) - \hat{\mathbf{z}}(t_{k + i})\right\| \le \rho_x^i( \left\| \hat{\mathbf{x}}(t_{k}) - \hat{\mathbf{z}}(t_{k})\right\|)$が成り立つ. よって, これらの性質を組み合わせると,

\begin{align}
\|V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}}) - V_{N_H}^\star ({\mathbf{z}}(t_{k})| \mathcal{D}_{N_{k}})\| &\le \sum_{i = 0}^{N_H - 1} \|L(\hat{\mathbf{x}}^\star(t_{k + i}), \hat{\mathbf{u}}^\star(t_{k + i})) - L(\hat{\mathbf{z}}^\star(t_{k + i}), \hat{\mathbf{u}}^\star(t_{k + i})) \| + \| \lambda V_{f}(\hat{\mathbf{x}}^\star(t_{k + N_H})) - \lambda V_{f}(\hat{\mathbf{z}}^\star(t_{k + N_H})) \|\notag \\
& \le \sum_{i = 0}^{N_H - 1} \rho_\ell \circ \rho_x^i (\|\mathbf{x}(t_k) - \mathbf{z}(t_k)\|)  + \rho_f \circ \rho_x^i (\|\mathbf{x}(t_k) - \mathbf{z}(t_k)\|)  \notag \\
&:= \rho_v (\|\mathbf{x}(t_k) - \mathbf{z}(t_k)\|)
\end{align}
\tag{30}

ここで, 式(25)より, $F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}), \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})$は$\Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})$についてアフィンであるので, $\Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})$において一様連続であることから,

\|F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),  \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k}) -  F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),  \Delta{\boldsymbol \sigma}_{\mathbf{f}}'(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k}) \| \le \rho_e (\|\Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})) - \Delta{\boldsymbol \sigma}_{\mathbf{f}}'(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})) \|)
\tag{31}

を満たす$\mathcal{K}$関数$\rho_e(\cdot)$が存在する. よって,　

\begin{align}
V_{N_H}^\star ({\mathbf{x}}(t_{k + 1})| \mathcal{D}_{N_{k + 1}}) - V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}}) =& V_{N_H}^\star (F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),  \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})| \mathcal{D}_{N_{k + 1}}) - V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}})\notag\\
=& V_{N_H}^\star (F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),  \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})| \mathcal{D}_{N_{k + 1}}) - V_{N_H}^\star (F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}), 0 | \mathcal{D}_{N_{k + 1}}) \notag\\
& + V_{N_H}^\star (F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}), 0 | \mathcal{D}_{N_{k + 1}}) - V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}}) \notag \\
\le& \|V_{N_H}^\star (F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),  \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k})| \mathcal{D}_{N_{k + 1}})-V_{N_H}^\star (F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}), 0 | \mathcal{D}_{N_{k + 1}})\| - \alpha_3(\|\mathbf{x}(t_k)\|) \notag \\
\le& \rho_v(\|F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),  \Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))| \mathcal{D}_{N_k}) - F(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}),0| \mathcal{D}_{N_k})\|) - \alpha_3(\|\mathbf{x}(t_k)\|) \notag \\
\le& \rho_v \circ \rho_e(\|\Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))\|) - \alpha_3(\|\mathbf{x}(t_k)\|)
\end{align}
\tag{32}

が導かれる. 以上の議論により, 次の定理が導かれる.

定理２: 仮定1~4と定理1を満たすモデル予測制御$ \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$が存在すると仮定する. このとき, $\forall \mathbf{x}(t_k) \in \Omega_\eta^{k}(\lambda) = \left\{\mathbf{x}(t_k)\in \mathcal{X} | V_{N_H}^\star ({\mathbf{x}}(t_k)| \mathcal{D}_{N_{k}}) \le \eta \right\} \subseteq\mathcal{X}_{N_H}^k(\lambda)$において, 未知のプラント$\mathbf{x}(t_{k + 1}) = \mathbf{f}(\mathbf{x}(t_k), \mathbf{u}(t_k))$にモデル予測制御$\mathbf{u}(t_k) = \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$を適用した時, モデルの誤差$\Delta{\boldsymbol \sigma}_{\mathbf{f}}(\mathbf{x}(t_k), \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k}))$に対して$V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}})$がISSリアプノフ関数となり, その閉ループシステムはISSである. 加えて, $\left\|Q_{\mathbf{f}, N_k}{\boldsymbol \xi}_{N_k}\right\|$が小さくなるほど$\Omega_\eta^{k}(\lambda)$は大きくなる.

安定性を保証した学習ベースモデル予測制御のアルゴリズム

上記の定理1, 定理２によりガウス過程により未知のプラントを学習し, 安定性を保証したモデル予測制御を生成するアルゴリズムを以下に示す.

アルゴリズム1
初期化事項1: 学習用データ$\mathcal{D}_k$を初期化し, GPモデルのハイパーパラメータ$\theta_\mathbf{f}$を式(7)に従い更新する.
初期化事項2: [9]に基づき, 原点でGPモデルを線形近似し, 終端コスト$V_f(\cdot)$を決定する.

各ステップ$k$ごとに2 ~ 7を繰り返す.
最適制御問題(8)を解き (MATLAB では fmincon を用いれば良い), モデル予測制御$\mathbf{u}(t_k) = \mathbf{u}_{\rm MPC}({\mathbf{x}}(t_k)| \mathcal{D}_{N_k})$を得る.
プラントにモデル予測制御を適用し, $\mathbf{x}(t_{k + 1}) = \mathbf{f}(\mathbf{x}(t_k),\mathbf{u}(t_k))$を得る.
GPモデルにモデル予測制御を適用し, $\hat{\mathbf{x}}(t_{k + 1}) = \hat{\mathbf{f}}(\mathbf{x}(t_k),\mathbf{u}(t_k))$を得る.
式(5)により, $\forall i = 1,\ldots, n$に対して, $\sigma_{f_i}^2(\mathbf{x}(t_k), \mathbf{u}(t_k)|\mathcal{D}_{i, N_k}) \le \bar{\sigma}_{\mathbf{f}}^2$, かつ, $\left\|\mathbf{x}(t_{k + 1}) - \hat{\mathbf{x}}(t_{k + 1})\right\| \le \bar{e}$ならばステップ2へ進む, それ以外の場合はステップ6に進む. (モデルが不確かな場合のみ学習させる.)
$\mathcal{D}_{N_{k + 1}}' \leftarrow \mathcal{D}_{N_k} \cup \left\{(\mathbf{x}(t_k), \mathbf{u}(t_k)), \mathbf{x}(t_{k + 1})\right\}$とおく.
定理1, 2より, $V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k + 1}})' \le V_{N_H}^\star ({\mathbf{x}}(t_{k})| \mathcal{D}_{N_{k}})$ならば, $\mathcal{D}_{N_{k + 1}} \leftarrow \mathcal{D}_{N_{k + 1}}'$とする. それ以外の場合, $\mathcal{D}_{N_{k + 1}} \leftarrow \mathcal{D}_{N_{k}}$とする.

Appendix

表記

ベクトルは太字で表現 (e.g., $\mathbf{x} = [x_1, \ldots, x_n]^\top \in \mathbb{R}^n$).

RKHSの特性について

ヒルベルト空間 : ベクトル空間で内積$\langle \ , \ \rangle$が与えられており, 完備性(コーシー列が常に収束する)を満たす.

$k(x, y)$が正定値カーネル(例えば, ガウスカーネル) $\implies$ $\Omega$上の関数で形成されるヒルベルト空間$H_k$で次の3つの条件を満たすものが一意に存在する.
- $x\in\Omega$を任意に固定したとき, $k(\cdot, x)\in H_k$.
- $f = \sum_{i = 1}^n c_i k(\cdot, x_i)$の元は$H_k$で稠密.
- $\forall f \in H_k$, $x\in\Omega$に対して, $f(x) = \langle f, k(\cdot, x)\rangle$.

References

[1] : K. Hashimoto, Y. Yoshimura, and T.Ushio, (2020). Learning self-triggered controllers with Gaussian processes. IEEE Transactions on Cybernetics.
[2]: K. Hashimoto, A. Saoud, M. Kishida, T. Ushio, and D. Dimarogonas, (2020). Learning-based safe symbolic abstractions for nonlinear control systems.
[3]: M. Maiworm, D. Limon, and R. Findeisen, (2021). Online learning‐based model predictive control with Gaussian process models and stability guarantees. International Journal of Robust and Nonlinear Control.
[4]: D. Limón, T. Alamo, F. Salas, and E. F. Camacho, (2006). On the stability of constrained MPC without terminal constraint. IEEE transactions on automatic control, 51(5), 832-836.
[5]: I. Steinwart, and A. Christmann, (2008). Support vector machines. Springer Science & Business Media.
[6]: D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert, Constrained modelpredictive control: Stability and optimality, Automatica, vol. 36, pp. 789–814, 2000.
[7]: Z. P. Jiang and Y. Wang, (2001). Input-to-state stability for discrete-time nonlinear systems. Automatica, 37(6), 857-869.
[8]: Limón D, Alamo T, Raimondo DM, et al. Input-to-state stability: a unifying framework for robust model predictive control. Nonlinear Model Predictive Control. New York, NY: Springer; 2009:1-26.
[9]: Maiworm M, Bäthge T, Findeisen R. Scenario-based model predictive control: recursive feasibility and stability. Paper presented at: Proceedings of the 9th International Symposium on Advanced Control of Chemical Processes (ADCHEM). Whistler, Canada: IFAC; 2015:50-56
[10]: M. Maiworm, D. Limon, and R. Findeisen (2021). Online learning‐based model predictive control with Gaussian process models and stability guarantees. International Journal of Robust and Nonlinear Control.