Parameter-Varying な系に対するKoopman作用素を用いたモデル化と制御：モデル化編

当記事は東京工業大学の藤井成美（https://www.linkedin.com/in/narumi-fujii-32304330b/ ）によって執筆されました

はじめに

本記事では，IEEE CDC 2023で発表されたParameter-Varying Koopman Operator for Nonlinear System Modeling and Controlという論文を紹介します．Koopman作用素論は，ある非線形系が在る時にそれに等価な時間発展を示す高次元空間内の線形系の存在を保証しています．このKoopman作用素は，元の非線形系から得られる時系列データを用いてeDMDという方法で推定され，非線形系をData-Driven型で線形系としてモデル化する場合に用いられます．Koopman作用素を用いたモデル化では，空間的に大域的な線形系を得られることが局所線形化と比較して利点となっています．とはいえ，今までKoopman作用素を用いたモデル化の研究では，主に時不変な非線形系を対象としていました(Korda and Mezić (2016))．この論文では，実現象の多くのダイナミクスが何らかのパラメータの値に依存することに着目し，時間依存パラメータ$p=p(t)p=p(t)$を含むNonlinear Parameter-Varying system (NPVs)に対するKoopman作用素を用いたモデル化が提案されています．時間依存パラメータを含む非線形系とは，例えば以下のような系のことです．

$$x˙=10(y−x)y˙=px−y−xzz˙=xy−z \begin{align} & \dot{x}=10(y-x) \\ & \dot{y}= px-y-x z \\ & \dot{z}=x y-z \end{align}$$

モデル化の具体的な方法としては，異なる$pp$の値に対して局所的に時不変であることを仮定してKoopman作用素を同定し，得られたKoopman作用素の間を補間することでNPVsに対するKoopman作用素を推定します．得られる線形系はLinear Parameter-Varying system (LPVs)となり，パラメータ方向には大域的でないものの空間方向の大域性を持ちます．さらに，得られたLPVsに対して最近発展しているLPV-MPCを統合することで元のNPVsの制御を行うことが提案されています．本記事は，「LPVsとしてのNPVsのモデル化」と「LPV-MPC適用による制御」の二部構成とし，ここでは「LPVsとしてのNPVsのモデル化」に焦点を当てて紹介します．また，Koopman作用素やeDMDの知識をある程度仮定しているため，気になる方はこちらの記事も参考にしてください(Koopman作用素，eDMD)．

Parameter-Varying Koopman作用素

先に述べたように，Koopman作用素は空間的に大域的な線形系を提供するので，NPVsに適用した場合はParameter-Varyingな性質が残りKoopman作用素もParameter-Varying Koopman作用素 (PV Koopman作用素)となります．以下では，PV Koopman作用素の推定方法について述べます．
まず，以下の離散時間的なNPVsを考えます^[1]．

$$xk+1=f(xk, uk, pk) \begin{align} {\bm x}_{k+1} = {f} ({\bm x}_k,~{\bm u}_k,~{p}_k) \tag{1} \end{align}$$

ここで，$xk∈X⊂Rn{\bm x}_k\in \mathbb{X}\subset\mathbb{R}^n$，$uk∈U⊂Rm{\bm u}_k\in \mathbb{U}\subset\mathbb{R}^m$は状態ベクトルと入力ベクトル，$pk∈P⊂R{ p}_k\in \mathbb{P}\subset\mathbb{R}$はパラメータで時間に依存します．$X, U, P\mathbb{X},~\mathbb{U},~\mathbb{P}$などは，入力はこの大きさまでと言うようなそれぞれの拘束条件で決まります．そして，$pkp_k$に依存するPV Koopman作用素は，

$$Kpk(Ψ(xk, uk))=Ψ(xk+1, uk+1) \begin{align} \mathcal{K}_{p_k}\left({\bm \Psi}\left({\bm x}_k,~{\bm u}_k\right)\right)={\bm \Psi}\left({\bm x}_{k+1},~{\bm u}_{k+1}\right) \tag{2} \end{align}$$

のようにして高次元空間における線形系を与えます．ここで，$Ψ∈G{\bm \Psi} \in \mathbb{G}$は観測関数，$Kpk:G→G\mathcal{K}_{p_k}: \mathbb{G}\rightarrow \mathbb{G}$はPV Koopman作用素です．添え字の$pk{p_k}$でパラメータに依存することを表しています．$G\mathbb{G}$は無限次元の関数空間ですが数値計算上扱えないので，設計者が有限個の適当な関数を基底として選びます．一般に状態に関する基底の数$qq$は元の系の状態の次元$nn$よりも大きく，$Ψ∈G:Rn+m→Rq+m{\bm \Psi}\in \mathbb{G}: \mathbb{R}^{n+m}\rightarrow\mathbb{R}^{q+m}$ ($q>nq> n$)のように次元を拡張し，

$$Ψ(xk, uk)=[ψ1(xk), ψ2(xk),⋯ , ψq(xk), uk⊤]⊤(3) {\bm \Psi}\left({\bm x}_k,~{\bm u}_k\right)=\left[\psi_1({\bm x}_k),~ \psi_2({\bm x}_k), \cdots,~ \psi_q({\bm x}_k),~{\bm u}_k^{\top}\right]^{\top} \tag{3}$$

ここで，$ψi:Rn→R\psi_i:\mathbb{R}^n \rightarrow\mathbb{R}$です．特に，高次元空間にリフトされた状態は，

$$yk=Ψx(xk)=[ψ1(xk), ψ2(xk),⋯ , ψq(xk)]⊤(4) {\bm y}_k={\bm \Psi}_{{\bm x}}\left({\bm x}_k\right)=\left[\psi_1\left({\bm x}_k\right),~\psi_2\left({\bm x}_k\right), \cdots, ~\psi_q\left({\bm x}_k\right)\right]^{\top} \tag{4}$$

Koopman作用素の推定には，$xk, uk, yk{\bm x}_k,~{\bm u}_k,~{\bm y}_k$の時系列データを用いてeDMDを行うというのが一般的ですが，今回はPV Koopman作用素を推定したいので少し修正が必要です．ここでは，パラメータは観測可能であり，さらに設計者が適当な$p1,⋯ , pl∈Pp^1,\cdots,~p^l\in\mathbb{P}$の$ll$個の時不変な値に設定できることを仮定して，各パラメータ値における時系列データを集めてそれぞれeDMDを行いKoopman作用素を推定します．PV Koopman作用素はそれらのKoopman作用素の間を補間することで推定されます (Fig.1)．

Fig. 1 PV Koopman作用素 (出典：Lee et al. (2023))

例えば，時不変なパラメータの$i(i=1,⋯ , l)i(i=1,\cdots,~l)$番目を$pip^i$とすると，$pk=pip_k=p^i$におけるKoopman作用素はeDMDを用いて次のように求められます．
対象とする系の式$(1)(1)$において$pk=pip_k=p^i$とした時に得られる状態と入力の$MM$時刻における時系列データ，

$$X(i)=[x1,x2,…,xM−1]∈Rn×(M−1)X+(i)=[x2,x3,…,xM]∈Rn×(M−1)U(i)=[u1,u2,…,uM−1]∈Rm×(M−1) \begin{align} {\bm X}(i) & =\left[{\bm x}_1, {\bm x}_2, \ldots, {\bm x}_{M-1}\right] \in \mathbb{R}^{n \times(M-1)} \tag{5a}\\ {\bm X}^{+}(i) & =\left[{\bm x}_2, {\bm x}_3, \ldots, {\bm x}_M\right] \in \mathbb{R}^{n \times(M-1)} \tag{5b}\\ {\bm U}(i) & =\left[{\bm u}_1, {\bm u}_2, \ldots, {\bm u}_{M-1}\right] \in \mathbb{R}^{m \times(M-1)} \tag{5c} \end{align}$$

を用意します．$X(i), X+(i){\bm X}(i),~{\bm X}^{+}(i)$についてはリフトされた状態に変換し，

$$Y(i)=[y1,y2,…,yM−1]∈Rq×(M−1)Y+(i)=[y2,y3,…,yM]∈Rq×(M−1) \begin{align} {\bm Y}(i) & =\left[{\bm y}_1, {\bm y}_2, \ldots, {\bm y}_{M-1}\right] \in \mathbb{R}^{q \times(M-1)} \tag{6a}\\ {\bm Y}^{+}(i) & =\left[{\bm y}_2, {\bm y}_3, \ldots, {\bm y}_M\right] \in \mathbb{R}^{q \times(M-1)}\tag{6b} \end{align}$$

所望のリフトされた高次元空間内の時不変な線形系は，

$$Y+(i)=A(pi)Y(i)+B(pi)U(i)X(i)=CY(i) \begin{align} {\bm Y}^{+}(i) & ={\bm A}\left(p^i\right) {\bm Y}(i)+{\bm B}\left(p^i\right) {\bm U}(i) \tag{7a}\\ {\bm X}(i) & ={\bm C} {\bm Y}(i)\tag{7b} \end{align}$$

を満たし，$A(pi), B(pi), C{\bm A}\left(p^i\right),~{\bm B}\left(p^i\right),~{\bm C}$は，

$$min⁡A(pi),B(pi)∣Y+(i)−(A(pi)Y(i)+B(pi)U(i))∣Fmin⁡C∣X(i)−CY(i)∣F \begin{align} \min _{{\bm A}\left(p^i\right), {\bm B}\left(p^i\right)}\left|{\bm Y}^{+}(i)-\left({\bm A}\left(p^i\right) {\bm Y}(i)+{\bm B}\left(p^i\right) {\bm U}(i)\right)\right|_F \tag{8a}\\ \min _C|{\bm X}(i)-{\bm C} {\bm Y}(i)|_F\tag{8b} \end{align}$$

のようなフロべニウスノルムの最小化より求められます．これは，擬似逆行列を用いて，

$$[A(pi)B(pi)]=Y+(i)[Y(i)U(i)]†C=X(i)Y(i)† \begin{align} {\left[{\bm A}\left(p^i\right) \quad {\bm B}\left(p^i\right)\right] } & ={\bm Y}^{+}(i)\left[\begin{array}{l} {\bm Y}(i) \\ {\bm U}(i) \end{array}\right]^{\dagger}\tag{9a}\\ {\bm C} & ={\bm X}(i){\bm Y}(i)^{\dagger}\tag{9b} \end{align}$$

式$(9a)(9a)$より各$pip^i$について計算した$A(pi){\bm A}\left(p^i\right)$と$B(pi){\bm B}\left(p^i\right)$を用いて，ある$pkp_k$に対する$A(pk){\bm A}\left(p_k\right)$と$B(pk){\bm B}\left(p_k\right)$を

$$A(pk)=α1(pk)A(p1)+α2(pk)A(p2)+⋯+αl(pk)A(pl)B(pk)=α1(pk)B(p1)+α2(pk)B(p2)+⋯+αl(pk)B(pl) \begin{align} & {\bm A}\left(p_k\right)=\alpha_1\left(p_k\right) {\bm A}\left(p^1\right)+\alpha_2\left(p_k\right) {\bm A}\left(p^2\right)+\cdots+\alpha_l\left(p_k\right) {\bm A}\left(p^l\right) \tag{10a}\\ & {\bm B}\left(p_k\right)=\alpha_1\left(p_k\right) {\bm B}\left(p^1\right)+\alpha_2\left(p_k\right) {\bm B}\left(p^2\right)+\cdots+\alpha_l\left(p_k\right) {\bm B}\left(p^l\right) \tag{10b} \end{align}$$

とします．ここで，$α1(pk), α2(pk),⋯ , αl(pk)\alpha_1\left(p_k\right),~\alpha_2\left(p_k\right),\cdots,~\alpha_l\left(p_k\right)$は重み係数であり$pkp_k$に依存します．重み係数の関数同定方法としては，線形補間が用いられています．以上より，PV Koopman作用素を用いてモデル化されたLPVsとして以下を得ることができます．

$$yk+1=A(pk)yk+B(pk)uk \begin{align} {\bm y}_{k+1}={\bm A}\left(p_k\right) {\bm y}_k+{\bm B}\left(p_k\right) \mathbf{u}_k \tag{11} \end{align}$$

現在の状態と未来のパラメータの変化が既知であれば同定したLPVsを用いてNPVsの未来を予測することができ，LPV-MPCを適用することができます．

予測精度の評価

今回は，モデルの予測精度を評価することが目的なので制御入力は加えておらず，推定するのは$A(pk){\bm A}\left(p_k\right)$の部分のみです．対象とするNPVsとして，以下のローレンツ方程式を考えます．

$$x˙=10(y−x)y˙=px−y−xzz˙=xy−z \begin{align} & \dot{x}=10(y-x)\tag{12a} \\ & \dot{y}=px-y-xz\tag{12b}\\ & \dot{z}=xy-z\tag{12c} \end{align}$$

ここで，時間依存するパラメータを，

$$p(t)=25+∑i=120aisin⁡(fit) ∑i=120ai=5 (ai>0) \begin{align} p(t)=25+\sum_{i=1}^{20}a_i\sin(f_it)~~~~~~ \sum_{i=1}^{20}a_i=5~~~(a_i>0) \tag{13} \end{align}$$

とし，$fif_i$は$00$から$1010$の一様分布から決定されるとします．この系に対して，従来のKoopman作用素と今回のPV Koopman作用素を推定し，得られたモデルの予測精度を比較します．

• Koopman作用素
  - データの時刻数: $M=15000M=15000$ ($150s150s$のデータ)
  - 基底関数 (radial basis function)の数: $q=50q=50$
• PV Koopman作用素
  • $pi=20, 25, 30 (l=3)p^i=20,~25,~30~~(l=3)$
  • データの時刻数: 各パラメータ$pip^i$ごとに$M=5000M=5000$とし全部で$1500015000$時刻 (各パラメータ$pip^i$ごとに$50s50s$のデータ)
  • 基底関数 (radial basis function)の数: $q=50q=50$

Fig2は，Koopman作用素(KO)とPV Koopman作用素(PVKO)の推定から得られた，時不変な線形系とLPVsのモデルついて2秒間の予測結果を示し，元のNPVsの結果(True)と比較しています．

Fig. 2 予測した状態の比較 (出典：[Lee et al. (2023)]

PV Koopman作用素の方がより元の系の結果に合っているように見えます．定量的には，基底関数 (radial basis function)の次数$NrbfN_{rbf}$を変化させ，それぞれの$NrbfN_{rbf}$ついて2秒間の予測シミュレーションを500回行った時の，

$$RMSE=100∑k∥x^k−xk∥22∑k∥xk∥22 RMSE=100\frac{\sqrt{\sum_k\|\hat{\bm x}_k-{\bm x}_k\|_2^2}}{\sqrt{\sum_k\|{\bm x}_k\|_2^2}}$$

を比較します．$RMSERMSE$は二乗平均平方根誤差であり，$x^k\hat{\bm x}_k$は予測した状態，$xk{\bm x}_k$は式$(12a)(12a)$から式$(12c)(12c)$によって得られる状態です．Fig.3に示すように，より多くの条件においてPV Koopman作用素の方がKoopman作用素に比べて二乗平均平方根誤差が小さい傾向にあります．以上より，PV Koopman作用素がNPVsな系に対して有効なLPVsモデルを与えることが分かりました．

Fig. 3 二乗平均平方根誤差の平均$±σ\pm\sigma$の比較 (出典：[Lee et al. (2023)]

最後に

次回は，得られたLPVsのモデルに対してLPV-MPCを適用する話をしたいと思います．

参考文献

Lee, Changyu, Kiyong Park and Jinwhan Kim. “Parameter-Varying Koopman Operator for Nonlinear System Modeling and Control.” 2023 62nd IEEE Conference on Decision and Control (CDC) (2023): 3700-3705.

Korda, Milan and Igor Mezić. “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control.” Autom. 93 (2016): 149-160.

脚注

1. 連続時間系のデータサンプリングによって得られる離散時間系を想定してます．