[DMD-3] DMDior

입출력이 포함된 확장 DMD인 DMDio (DMD with input/output) 알고리즘을 유도해 보았다 (https://pasus.tistory.com/225). 원래 시스템을 식별한 후에 축소 모델 (ROM, reduced order model)로 근사화 하는 순서였다.

 

 

이번에는 이와 약간 다른 접근 방법을 사용해 보고자 한다. 바로 축소 모델을 식별하는 방법이다. 이러한 방법을 DMDior (DMDio for reduced order model)라고 한다.

식별하고자 하는 미지의 이산시간 시스템이 식 (1)과 같이 표현된다고 하자.

 

$$\begin{array}{r}\text{(1)}& {\mathbf{x}}_{k+1}& =A{\mathbf{x}}_k+B{\mathbf{u}}_k\\ {\mathbf{y}}_k& =C{\mathbf{x}}_k+D{\mathbf{u}}_k\end{array}$$

 

여기서 ${\mathbf{x}}_k\in {\mathbb{R}}^n$, ${\mathbf{u}}_k\in {\mathbb{R}}^p$, ${\mathbf{y}}_k\in {\mathbb{R}}^q$, $A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times p}$, $C\in {\mathbb{R}}^{q\times n}$, $D\in {\mathbb{R}}^{q\times p}$ 이다. 입력과 출력의 차원 $p$ 와 $q$ 는 일반적으로 매우 작은 값으로 $p,q\ll n$ 이다.

DMDio와 마찬가지로 일정한 싯점에서 $m$ 개의 상태변수 스냅샷 데이터 ${\mathbf{x}}_k$ 와 $(m-1)$ 개의 입력과 출력의 스냅샷 데이터 ${\mathbf{u}}_k$ 와 ${\mathbf{y}}_k$ 를 수집하여 이를 다음과 같이 행렬 형태로 표현한다. 여기서 ${\mathbf{u}}_k$ 는 임의의 입력이다.

 

$$\begin{array}{r}\text{(2)}& \mathbf{X}& =[{\mathbf{x}}_1{\mathbf{x}}_2\cdots {\mathbf{x}}_{m-1}]\in {\mathbb{R}}^{n\times (m-1)}\\ {\mathbf{X}}^{\prime }& =[{\mathbf{x}}_2{\mathbf{x}}_3\cdots {\mathbf{x}}_m]\in {\mathbb{R}}^{n\times (m-1)}\\ \\ \mathbf{U}& =[{\mathbf{u}}_1{\mathbf{u}}_2\cdots {\mathbf{u}}_{m-1}]\in {\mathbb{R}}^{p\times (m-1)}\\ \\ \mathbf{Y}& =[{\mathbf{y}}_1{\mathbf{y}}_2\cdots {\mathbf{y}}_{m-1}]\in {\mathbb{R}}^{q\times (m-1)}\end{array}$$

 

여기서 ${\mathbf{X}}^{\prime }$ 은 $\mathbf{X}$ 를 한 시간스텝 만큼 시프트(shift)시킨 행렬이다.

 

 

이제 바로 $n$ 차 모델을 $r$ 차 축소 모델 (ROM, reduced order model)로 근사화 하자. 여기서 r≪n이다. 선형 기저 변환(basis transformation)을 통해서 $n$차 벡터 ${\mathbf{x}}_k$ 를 $r$ 차 벡터 ${\tilde{\mathbf{x}}}_k$ 로 축소 변환한다.

 

$$\begin{array}{}\text{(3)}& {\mathbf{x}}_k=P{\tilde{\mathbf{x}}}_k\end{array}$$

 

여기서 $P\in {\mathbb{R}}^{n\times r}$ 는 변환 행렬로서 $P^{T}P=I$ 를 만족한다.

 

 

표준 DMD와 마찬가지로 스냅샷 POD(proper orthogonal decomposition) 알고리즘을 이용하여 변환 행렬을 결정하도록 한다. 상태변수의 변환 행렬이므로 행렬 $\mathbf{X}$ 의 SVD의 왼쪽 특이 벡터(left singular vector)를 변환 행렬로 사용한다. 행렬 $\mathbf{X}$ 의 SVD를 계산하면 다음과 같다.

 

$$\begin{array}{r}\text{(4)}& \mathbf{X}& =U\mathrm{\Sigma }{V}^T=\left[\begin{array}{cc}{U}_r& U_{rem}\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Sigma }}_r& 0\\ 0& {\mathrm{\Sigma }}_{rem}\end{array}\right]\left[\begin{array}{c}{V}_r^T\\ V_{rem}^T\end{array}\right]\\ & \approx U_r{\mathrm{\Sigma }}_r{V}_r^T\end{array}$$

 

여기서 $U\in {\mathbb{R}}^{n\times n}$, $\mathrm{\Sigma }\in {\mathbb{R}}^{n\times (m-1)}$, $V\in {\mathbb{R}}^{(m-1)\times (m-1)}$, $U_r\in {\mathbb{R}}^{n\times r}$, ${\mathrm{\Sigma }}_r\in {\mathbb{R}}^{r\times r}$, $V_r\in R^{(m-1)\times r}$ 이다. $r$ 은 시스템의 차원을 $n$ 에서 $r$ 로 줄이기 위한 절단값(truncation value)으로서 설계 변수가 된다.

 

 

이제 $P=U_r$ 로 놓고 식 (3)을 식 (1)에 대입하면 다음과 같이 된다.

 

$$\begin{array}{r}\text{(5)}& {\tilde{\mathbf{x}}}_{k+1}& =U_r^T{\mathbf{x}}_{k+1}=U_r^{T}A{U}_r{\tilde{\mathbf{x}}}_k+U_r^{T}B{\mathbf{u}}_k\\ & =\tilde{A}{\tilde{\mathbf{x}}}_k+\tilde{B}{\mathbf{u}}_k\\ \\ {\mathbf{y}}_k& =C{U}_r{\tilde{\mathbf{x}}}_k+D{\mathbf{u}}_k\\ \\ & =\tilde{C}{\tilde{\mathbf{x}}}_k+D{\mathbf{u}}_k\end{array}$$

 

한편, 식 (1)의 시스템 운동식을 이용하면 다음과 같은 관계식을 유도할 수 있다.

 

$$\begin{array}{r}\text{(6)}& {\mathbf{X}}^{\prime }& =[A{\mathbf{x}}_1+B{\mathbf{u}}_{1}A{\mathbf{x}}_2+B{\mathbf{u}}_2\cdots A{\mathbf{x}}_{m-1}+B{\mathbf{u}}_{m-1}]\\ & =A\mathbf{X}+B\mathbf{U}\\ \\ \mathbf{Y}& =[C{\mathbf{x}}_1+D{\mathbf{u}}_{1}C{\mathbf{x}}_2+D{\mathbf{u}}_2\cdots C{\mathbf{x}}_{m-1}+D{\mathbf{u}}_{m-1}]\\ \\ & =C\mathbf{X}+D\mathbf{U}\end{array}$$

 

식 (5)와 (6)으로부터 다음과 같은 식을 유도할 수 있다.

 

$$\begin{array}{r}\text{(7)}& \left[\begin{array}{c}{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]& =\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\left[\begin{array}{c}\mathbf{X}\\ \mathbf{U}\end{array}\right]\\ & \approx \left[\begin{array}{cc}{U}_r\tilde{A}{U}_r^T& U_r\tilde{B}\\ \tilde{C}{U}_r^T& D\end{array}\right]\left[\begin{array}{c}\mathbf{X}\\ \mathbf{U}\end{array}\right]\\ \\ & =\left[\begin{array}{cc}{U}_r& 0\\ 0& I\end{array}\right]\left[\begin{array}{cc}\tilde{A}& \tilde{B}\\ \tilde{C}& D\end{array}\right]\left[\begin{array}{cc}{U}_r^T& 0\\ 0& I\end{array}\right]\left[\begin{array}{c}\mathbf{X}\\ \mathbf{U}\end{array}\right]\end{array}$$

 

또는

 

$$\begin{array}{}\text{(8)}& \left[\begin{array}{c}{U}_r^T{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]\approx \left[\begin{array}{cc}\tilde{A}& \tilde{B}\\ \tilde{C}& D\end{array}\right]\left[\begin{array}{c}{U}_r^T\mathbf{X}\\ \mathbf{U}\end{array}\right]\end{array}$$

 

이다. 여기서 $\mathbf{Z}=\left[\begin{array}{c}{U}_r^T{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]\in {\mathbb{R}}^{(r+q)\times (m-1)}$, $\mathbf{F}=\left[\begin{array}{cc}\tilde{A}& \tilde{B}\\ \tilde{C}& D\end{array}\right]\in {\mathbb{R}}^{(r+q)\times (r+p)}$, $\mathbf{\Omega }=\left[\begin{array}{c}{U}_r^T\mathbf{X}\\ \mathbf{U}\end{array}\right]\in {\mathbb{R}}^{(r+p)\times (m-1)}$ 로 놓는다.

DMDio에서와 달리 축소 모델인 $\mathbf{F}$ 를 직접 식별하는 것이 DMDior의 차이점이다. 프로베니어스(Frobenius) 놈이 최소화되도록 $\mathbf{F}$ 를 다음과 같이 계산할 수 있다.

 

$$\begin{array}{r}\text{(9)}& \mathbf{F}& =\arg \underset{\mathbf{F}}{min}\Vert \mathbf{Z}-\mathbf{F}\mathbf{\Omega }{\Vert }_F=\mathbf{Z}{\mathbf{\Omega }}^{+}\\ & =\left[\begin{array}{c}{U}_r^T{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]{\left[\begin{array}{c}{U}_r^T\mathbf{X}\\ \mathbf{U}\end{array}\right]}^{+}\end{array}$$

 

식 (4)를 (9)에 대입하면 다음과 같다.

 

$$\begin{array}{r}\text{(10)}& \mathbf{F}& \approx \left[\begin{array}{c}{U}_r^T{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]{\left[\begin{array}{c}{U}_r^T{U}_r{\mathrm{\Sigma }}_r{V}_r^T\\ \mathbf{U}\end{array}\right]}^{+}\\ & =\left[\begin{array}{c}{U}_r^T{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]{\left[\begin{array}{c}{\mathrm{\Sigma }}_r{V}_r^T\\ \mathbf{U}\end{array}\right]}^{+}\\ \\ & =\mathbf{Z}{\tilde{\mathbf{\Omega }}}^{+}\end{array}$$

 

여기서 ${\tilde{\mathbf{\Omega }}}^{+}\in {\mathbb{R}}^{(m-1)\times (r+p)}$ 는 무어-펜로즈 유사 역행렬 (Moore-Penrose pseudo inverse matrix)이다. 유사 역행렬은 특이값 분해(SVD, singular value decomposition)를 이용해서 계산할 수 있다. 행렬 $\tilde{\mathbf{\Omega }}$ 의 SVD는 다음과 같다.

 

$$\begin{array}{r}\text{(11)}& \tilde{\mathbf{\Omega }}& =\hat{U}\hat{\mathrm{\Sigma }}{\hat{V}}^T=\left[\begin{array}{cc}{\hat{U}}_s& {\hat{U}}_{rem}\end{array}\right]\left[\begin{array}{cc}{\hat{\mathrm{\Sigma }}}_s& 0\\ 0& {\hat{\mathrm{\Sigma }}}_{rem}\end{array}\right]\left[\begin{array}{c}{\hat{V}}_s^T\\ {\hat{V}}_{rem}^T\end{array}\right]\\ & \approx {\hat{U}}_s{\hat{\mathrm{\Sigma }}}_s{\hat{V}}_s^T\end{array}$$

 

여기서 $\hat{U}\in {\mathbb{R}}^{(r+p)\times (r+p)}$, $\hat{\mathrm{\Sigma }}\in {\mathbb{R}}^{(r+p)\times (m-1)}$, $\hat{V}\in {\mathbb{R}}^{(m-1)\times (m-1)}$, ${\hat{U}}_s\in {\mathbb{R}}^{(r+p)\times s}$, ${\hat{\mathrm{\Sigma }}}_s\in {\mathbb{R}}^{s\times s}$, ${\hat{V}}_s\in {\mathbb{R}}^{(m-1)\times s}$ 이다. $s$ 는 절단값(truncation value)으로서 설계 변수가 된다.

그러면 유사 역행렬 ${\tilde{\mathbf{\Omega }}}^{+}$ 은 다음과 같이 근사적으로 계산할 수 있다.

 

$$\begin{array}{}\text{(12)}& {\tilde{\mathbf{\Omega }}}^{+}\approx {\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_s^T\in {\mathbb{R}}^{(m-1)\times (r+p)}\end{array}$$

 

식 (12)를 (10)에 대입하면 $\mathbf{F}$ 는 다음과 같이 된다.

 

$$\begin{array}{r}\text{(13)}& \mathbf{F}& =\left[\begin{array}{cc}\tilde{A}& \tilde{B}\\ \tilde{C}& D\end{array}\right]\approx \left[\begin{array}{c}{U}_r^T{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]\left[\begin{array}{cc}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T& {\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\end{array}\right]\\ & =\left[\begin{array}{cc}{U}_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T& U_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\\ \mathbf{Y}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T& \mathbf{Y}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\end{array}\right]\end{array}$$

 

여기서 ${\hat{U}}_s^T=\left[\begin{array}{cc}{\hat{U}}_1^T& {\hat{U}}_2^T\end{array}\right]$ 이고 ${\hat{U}}_1\in {\mathbb{R}}^{r\times s}$, ${\hat{U}}_2\in {\mathbb{R}}^{p\times s}$ 이다.

따라서 행렬 $\tilde{A},\tilde{B},\tilde{C},D$ 는 각각 다음과 같이 식별할 수 있다.

 

$$\begin{array}{r}\text{(14)}& \tilde{A}& =U_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T\\ \tilde{B}& =U_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\\ \\ \tilde{C}& =\mathbf{Y}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T\\ \\ D& =\mathbf{Y}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\end{array}$$

 

 

시스템의 동특성을 파악하기 위해서 축소 시스템 행렬 $\tilde{A}$ 의 고유값 ${\lambda }_i$ 와 고유벡터 ${\mathbf{w}}_i$ 를 먼저 계산한다. 축소 시스템 행렬 $\tilde{A}$ 의 고유값과 고유벡터 식은 다음과 같다.

 

$$\begin{array}{r}\text{(15)}& \tilde{A}{\mathbf{w}}_i& =U_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T{\mathbf{w}}_i\\ & ={\lambda }_i{\mathbf{w}}_i,i=1,...,r\end{array}$$

 

이제 벡터 ${\varphi }_i$ 를 다음과 같이 정의하자.

 

$$\begin{array}{}\text{(16)}& {\varphi }_i={\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T{\mathbf{w}}_i\end{array}$$

 

그러면 식 (15)는 다음과 같이 된다.

 

$$\begin{array}{}\text{(17)}& U_r^T{\varphi }_i={\lambda }_i{\mathbf{w}}_i\end{array}$$

 

위 식의 양변에 $U_r{U}_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T$ 을 곱하면,

 

$$\begin{array}{}\text{(18)}& U_r{U}_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T{U}_r^T{\varphi }_i={\lambda }_i{U}_r{U}_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T{\mathbf{w}}_i\end{array}$$

 

이 된다. 식 (7)과 (16)에 의하면,

 

$$\begin{array}{}\text{(19)}& A{\varphi }_i={\lambda }_i{U}_r{U}_r^T{\varphi }_i\approx {\lambda }_i{\varphi }_i,i=1,...,r\end{array}$$

 

가 된다.

위 식은 행렬 $A$ 의 고유값과 고유벡터 식이다. 따라서 원래 시스템 $A$ 와 축소 시스템 $\tilde{A}$ 의 고유값은 근사적으로 같고 고유벡터는 식 (16)의 관계가 성립한다는 것을 알 수 있다. 원래 시스템 $A$ 의 고유벡터 ${\varphi }_i$ 를 DMD 모드라고 한다.

DMDior 방법을 정리하면 다음과 같다.

      1. 먼저 $m$ 개의 스냅샷 데이터 ${\mathbf{x}}_k$ 와 $(m-1)$ 개의 스냅샷 데이터 ${\mathbf{u}}_k$, ${\mathbf{y}}_k$ 를 수집하여 이를 다음과 같이 행렬 형태로 표현한다.

$$\begin{array}{rl}& \mathbf{X}=[{\mathbf{x}}_1{\mathbf{x}}_2\cdots {\mathbf{x}}_{m-1}]\in {\mathbb{R}}^{n\times (m-1)}\\ \\ & {\mathbf{X}}^{\prime }=[{\mathbf{x}}_2{\mathbf{x}}_3\cdots {\mathbf{x}}_m]\in {\mathbb{R}}^{n\times (m-1)}\\ \\ & \mathbf{U}=[{\mathbf{u}}_1{\mathbf{u}}_2\cdots {\mathbf{u}}_{m-1}]\in {\mathbb{R}}^{p\times (m-1)}\\ \\ & \mathbf{Y}=[{\mathbf{y}}_1{\mathbf{y}}_2\cdots {\mathbf{y}}_{m-1}]\in {\mathbb{R}}^{q\times (m-1)}\end{array}$$

 

      2. 행렬 $\mathbf{X}$ 의 SVD를 계산한다.

$$\begin{array}{rl}\mathbf{X}& =U\mathrm{\Sigma }{V}^T=\left[\begin{array}{cc}{U}_r& U_{rem}\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Sigma }}_r& 0\\ 0& {\mathrm{\Sigma }}_{rem}\end{array}\right]\left[\begin{array}{c}{V}_r^T\\ V_{rem}^T\end{array}\right]\\ \\ & \approx U_r{\mathrm{\Sigma }}_r{V}_r^T\end{array}$$

 

      3. 행렬 $\tilde{\mathbf{\Omega }}=\left[\begin{array}{c}{\mathrm{\Sigma }}_r{V}_r^T\\ \mathbf{U}\end{array}\right]$ 의 SVD를 계산한다.

$$\begin{array}{rl}\tilde{\mathbf{\Omega }}& =\hat{U}\hat{\mathrm{\Sigma }}{\hat{V}}^T=\left[\begin{array}{cc}{\hat{U}}_s& {\hat{U}}_{rem}\end{array}\right]\left[\begin{array}{cc}{\hat{\mathrm{\Sigma }}}_s& 0\\ 0& {\hat{\mathrm{\Sigma }}}_{rem}\end{array}\right]\left[\begin{array}{c}{\hat{V}}_s^T\\ {\hat{V}}_{rem}^T\end{array}\right]\\ \\ & \approx {\hat{U}}_s{\hat{\mathrm{\Sigma }}}_s{\hat{V}}_s^T\\ \\ & =\left[\begin{array}{c}{\hat{U}}_1\\ {\hat{U}}_2\end{array}\right]{\hat{\mathrm{\Sigma }}}_s{\hat{V}}_s^T\end{array}$$

 

      4. 축소 시스템을 구한다.

$$\begin{array}{rl}{\tilde{\mathbf{x}}}_{k+1}& =\tilde{A}{\tilde{\mathbf{x}}}_k+\tilde{B}{\mathbf{u}}_k\\ {\mathbf{y}}_k& =\tilde{C}{\tilde{\mathbf{x}}}_k+D{\mathbf{u}}_k\\ \\ {\mathbf{x}}_k& =U_r{\tilde{\mathbf{x}}}_k\\ \tilde{A}& =U_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T\\ \tilde{B}& =U_r^T{\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\\ \tilde{C}& =\mathbf{Y}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T\\ D& =\mathbf{Y}{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_2^T\end{array}$$

 

      5. 행렬 $\tilde{A}$ 의 고유값과 고유벡터를 구한다.

$$\begin{array}{r}\tilde{A}{\mathbf{w}}_i={\lambda }_i{\mathbf{w}}_i,i=1,...,r\end{array}$$

 

      6. DMDior 모드를 계산한다.

$${\varphi }_i={\mathbf{X}}^{\prime }{\hat{V}}_s{\hat{\mathrm{\Sigma }}}_s^{-1}{\hat{U}}_1^T{\mathbf{w}}_i$$