[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{align} \mathbf{x}_{k+1} &= A \mathbf{x}_k+B \mathbf{u}_k \tag{1} \\ \\ \mathbf{y}_k &= C \mathbf{x}_k+D \mathbf{u}_k \end{align} \]

 

여기서 \(\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≪n\) 이다.

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

 

\[ \begin{align} \mathbf{X} &=[\mathbf{x}_1 \ \ \ \mathbf{x}_2 \ \ \ \cdots \ \ \ \mathbf{x}_{m-1} ] \ \ \ \in \mathbb{R}^{n \times (m-1)} \tag{2} \\ \\ \mathbf{X}' &=[\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{align} \]

 

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

 

 

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

 

\[ \mathbf{x}_k=P \tilde{\mathbf{x}}_k \tag{3} \]

 

여기서 \(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{align} \mathbf{X} &= U \Sigma V^T = \begin{bmatrix} U_r & U_{rem} \end{bmatrix} \begin{bmatrix} \Sigma_r & 0 \\ 0 & \Sigma_{rem} \end{bmatrix} \begin{bmatrix} V_r^T \\ V_{rem}^T \end{bmatrix} \tag{4} \\ \\ & \approx U_r \Sigma_r V_r^T \end{align} \]

 

여기서 \(U \in \mathbb{R}^{n \times n}\), \(\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}\), \(\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{align} \tilde{\mathbf{x}}_{k+1} &= U_r^T \mathbf{x}_{k+1}=U_r^T AU_r \tilde{\mathbf{x}}_k+U_r^T B \mathbf{u}_k \tag{5} \\ \\ &= \tilde{A} \tilde{\mathbf{x}}_k + \tilde{B} \mathbf{u}_k \\ \\ \mathbf{y}_k &=CU_r \tilde{\mathbf{x}}_k+D \mathbf{u}_k \\ \\ &= \tilde{C} \tilde{\mathbf{x}}_k+D \mathbf{u}_k \end{align} \]

 

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

 

\[ \begin{align} \mathbf{X}' &= [ 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} \ ] \tag{6} \\ \\ &=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{align} \]

 

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

 

\[ \begin{align} \begin{bmatrix} \mathbf{X}'\\ \mathbf{Y} \end{bmatrix} &= \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{U} \end{bmatrix} \tag{7} \\ \\ & \approx \begin{bmatrix} U_r \tilde{A} U_r^T & U_r \tilde{B} \\ \tilde{C} U_r^T & D \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{U} \end{bmatrix} \\ \\ &= \begin{bmatrix} U_r & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} \tilde{A} & \tilde{B} \\ \tilde{C} & D \end{bmatrix} \begin{bmatrix} U_r^T & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} \mathbf{X} \\ \mathbf{U} \end{bmatrix} \end{align} \]

 

또는

 

\[ \begin{bmatrix} U_r^T \mathbf{X}'\\ \mathbf{Y} \end{bmatrix} \approx \begin{bmatrix} \tilde{A} & \tilde{B} \\ \tilde{C} & D \end{bmatrix} \begin{bmatrix} U_r^T \mathbf{X} \\ \mathbf{U} \end{bmatrix} \tag{8} \]

 

이다. 여기서 \(\mathbf{Z}= \begin{bmatrix} U_r^T\mathbf{X}' \\ \mathbf{Y} \end{bmatrix} \in \mathbb{R}^{(r+q) \times (m-1)}\), \(\mathbf{F}= \begin{bmatrix} \tilde{A} & \tilde{B} \\ \tilde{C} & D \end{bmatrix} \in \mathbb{R}^{(r+q) \times (r+p)}\), \(\mathbf{\Omega}= \begin{bmatrix} U_r^T \mathbf{X} \\ \mathbf{U} \end{bmatrix} \in \mathbb{R}^{(r+p) \times (m-1)}\) 로 놓는다.

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

 

\[ \begin{align} \mathbf{F} &= \arg \min_{\mathbf{F}} \lVert \mathbf{Z}- \mathbf{F} \mathbf{\Omega} \rVert _F = \mathbf{Z} \mathbf{\Omega}^+ \tag{9} \\ \\ &= \begin{bmatrix} U_r^T \mathbf{X}' \\ \mathbf{Y} \end{bmatrix} \begin{bmatrix} U_r^T \mathbf{X} \\ \mathbf{U} \end{bmatrix} ^+ \end{align} \]

 

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

 

\[ \begin{align} \mathbf{F} & \approx \begin{bmatrix} U_r^T \mathbf{X}' \\ \mathbf{Y} \end{bmatrix} \begin{bmatrix} U_r^T U_r \Sigma_r V_r^T \\ \mathbf{U} \end{bmatrix} ^+ \tag{10} \\ \\ &= \begin{bmatrix} U_r^T \mathbf{X}' \\ \mathbf{Y} \end{bmatrix} \begin{bmatrix} \Sigma_r V_r^T \\ \mathbf{U} \end{bmatrix}^+ \\ \\ &= \mathbf{Z} \tilde{\mathbf{Ω}}^+ \end{align} \]

 

여기서 \(\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{align} \tilde{\mathbf{\Omega}} &= \hat{U} \hat{\Sigma} \hat{V}^T = \begin{bmatrix} \hat{U}_s & \hat{U}_{rem} \end{bmatrix} \begin{bmatrix} \hat{\Sigma}_s & 0 \\ 0 & \hat{\Sigma}_{rem} \end{bmatrix} \begin{bmatrix} \hat{V}_s^T \\ \hat{V}_{rem}^T \end{bmatrix} \tag{11} \\ \\ & \approx \hat{U}_s \hat{\Sigma}_s \hat{V}_s^T \end{align} \]

 

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

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

 

\[ \tilde{\mathbf{\Omega}}^+ \approx \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_s^T \ \ \ \in \mathbb{R}^{(m-1) \times (r+p)} \tag{12} \]

 

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

 

\[ \begin{align} \mathbf{F} &= \begin{bmatrix} \tilde{A} & \tilde{B} \\ \tilde{C} & D \end{bmatrix} \approx \begin{bmatrix} U_r^T \mathbf{X}' \\ \mathbf{Y} \end{bmatrix} \begin{bmatrix} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T & \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \end{bmatrix} \tag{13} \\ \\ &= \begin{bmatrix} U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T & U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \\ \mathbf{Y} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T & \mathbf{Y} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \end{bmatrix} \end{align} \]

 

여기서 \(\hat{U}_s^T= \begin{bmatrix} \hat{U}_1^T & \hat{U}_2^T \end{bmatrix}\) 이고 \(\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{align} \tilde{A} &= U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \tag{14} \\ \\ \tilde{B} &= U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \\ \\ \tilde{C} &= \mathbf{Y} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \\ \\ D &= \mathbf{Y} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \end{align} \]

 

 

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

 

\[ \begin{align} \tilde{A} \mathbf{w}_i &= U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \mathbf{w}_i \tag{15} \\ \\ &= \lambda_i \mathbf{w}_i , \ \ \ \ \ i=1, \ ..., \ r \end{align} \]

 

이제 벡터 \(\phi_i\) 를 다음과 같이 정의하자.

 

\[ \phi_i = \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \mathbf{w}_i \tag{16} \]

 

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

 

\[ U_r^T \phi_i= \lambda_i \mathbf{w}_i \tag{17} \]

 

위 식의 양변에 \( U_r U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \) 을 곱하면,

 

\[ U_r U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T U_r^T \phi_i= \lambda_i U_r U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \mathbf{w}_i \tag{18} \]

 

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

 

\[ A \phi_i= \lambda_i U_r U_r^T \phi_i \approx \lambda_i \phi_i , \ \ \ \ \ i=1, \ ..., \ r \tag{19} \]

 

가 된다.

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

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

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

\[ \begin{align} & \mathbf{X} =[\mathbf{x}_1 \ \ \ \mathbf{x}_2 \ \ \ \cdots \ \ \ \mathbf{x}_{m-1} ] \ \ \ \in \mathbb{R}^{n \times (m-1)} \\ \\ & \mathbf{X}' =[\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{align} \]

 

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

\[ \begin{align} \mathbf{X} &= U \Sigma V^T = \begin{bmatrix} U_r & U_{rem} \end{bmatrix} \begin{bmatrix} \Sigma_r & 0 \\ 0 & \Sigma_{rem} \end{bmatrix} \begin{bmatrix} V_r^T \\ V_{rem}^T \end{bmatrix} \\ \\ & \approx U_r \Sigma_r V_r^T \end{align} \]

 

      3. 행렬 \(\tilde{\mathbf{\Omega}}= \begin{bmatrix} \Sigma_r V_r^T \\ \mathbf{U} \end{bmatrix}\) 의 SVD를 계산한다.

\[ \begin{align} \tilde{\mathbf{\Omega}} &= \hat{U} \hat{\Sigma} \hat{V}^T = \begin{bmatrix} \hat{U}_s & \hat{U}_{rem} \end{bmatrix} \begin{bmatrix} \hat{\Sigma}_s & 0 \\ 0 & \hat{\Sigma}_{rem} \end{bmatrix} \begin{bmatrix} \hat{V}_s^T \\ \hat{V}_{rem}^T \end{bmatrix} \\ \\ & \approx \hat{U}_s \hat{\Sigma}_s \hat{V}_s^T \\ \\ & = \begin{bmatrix} \hat{U}_1 \\ \hat{U}_2 \end{bmatrix} \hat{\Sigma}_s \hat{V}_s^T \end{align} \]

 

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

\[ \begin{align} \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}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \\ \tilde{B} &= U_r^T \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \\ \tilde{C} &= \mathbf{Y} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \\ D &= \mathbf{Y} \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_2^T \end{align} \]

 

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

\[ \begin{align} \tilde{A} \mathbf{w}_i = \lambda_i \mathbf{w}_i , \ \ \ \ \ i=1, \ ..., \ r \end{align} \]

 

      6. DMDior 모드를 계산한다.

\[ \phi_i = \mathbf{X}' \hat{V}_s \hat{\Sigma}_s^{-1} \hat{U}_1^T \mathbf{w}_i \]