본문 바로가기
유도항법제어/데이터기반제어

[DMD-1] 동적모드분해 (Dynamic Mode Decomposition)

by 세인트 워터멜론 2022. 10. 26.

전통적인 제어이론은 시스템의 수학적인 운동 모델을 요구한다.

 

 

운동 모델은 물리 법칙으로부터 해석적으로 유도할 수 있지만 입출력 데이터에 기반해서 수치적으로 얻을 수도 있다. 수치 데이터로부터 시스템의 운동 모델을 구하는 것을 시스템 식별(system identification) 또는 모델 식별이라고 한다. 시스템 식별 방법에는 ERA, OKID, QMC등 몇 가지가 있는데, 그 중 하나가 동적모드분해 (DMD, dynamic mode decomposition)이다.

DMD는 수치 시뮬레이션 또는 스냅샷(snapshot) 측정 데이터를 사용하여 선형 시스템의 수학적 모델을 식별하고 동적 특성을 추출하는 기법이다. 식별하고자 하는 미지의 이산시간 시스템이 식 (1)과 같이 표현된다고 하자. 일단 자율 시스템(autonomous system)부터 시작하고 차츰 입력이 있는 시스템과 출력까지 있는 시스템으로 확장하고자 한다.

 

\[ \mathbf{x}_{k+1}=A \mathbf{x}_k \tag{1} \]

 

여기서 \(\mathbf{x}_k \in \mathbb{R}^n\), \(A \in \mathbb{R}^{n \times n}\) 이다. 보통 시스템의 상태변수는 연속적이지만 연속 상태변수 \(\mathbf{x}(t)\) 의 측정값은 \(\mathbf{x}_k=\mathbf{x}(k \Delta t)\) 로 표시되는 규칙적인 시간 간격 \(\Delta t\) 에서 수집할 수 있다. 이와 같이 일정한 싯점에서 수집한 데이터 \(\mathbf{x}_k\) 를 스냅샷(snapshot)이라고 한다.

예를 들면 유체의 운동 해석 문제의 경우 공간을 그리드로 분할하여 일정한 시간 간격마다 수치해석으로 계산된 데이터나 또는 측정된 데이터를 얻을 수 있다. 이 때 유동 데이터(속도나 압력)를 전체 그리드 포인트에서 시간 \(t=t_1, \ t_2, \ ..., \ t_m\) 마다 수집하여 스냅샷 데이터 벡터 \(\mathbf{x}_k\) 로 만들 수 있다.

 

 

여기서 \(m\) 은 스냅샷 갯수다. 보통 시스템의 차수는 \(n \gt 10,000\) 인 정도로 매우 크고 스냅샷 갯수는 \( m \lt 1,000\) 정도로 상대적으로 작은 경우가 많다. DMD의 목적은 식 (1)에서 상태변수 \(\mathbf{x}_k\) 의 측정 데이터를 기반으로 시스템 행렬 \(A\)를 추정하는 것이다.

먼저 \(m\) 개의 스냅샷 데이터 \(\mathbf{x}_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)} \end{align} \]

 

여기서 \(\mathbf{X}'\) 은 \(\mathbf{X}\) 를 한 시간스텝 만큼 시프트(shift)시킨 행렬로서 식 (1)의 시스템 운동식을 이용하면 다음과 같은 관계식을 유도할 수 있다.

 

\[ \begin{align} \mathbf{X}' &=[ A\mathbf{x}_1 \ A \mathbf{x}_2 \ \cdots \ A \mathbf{x}_{m-1} ] \tag{3} \\ \\ & = A \mathbf{X} \end{align} \]

 

식 (3)에 의하면 \(A\) 는 다음과 같이 계산할 수 있다.

 

\[ A= \arg \min_A \lVert \mathbf{X}'-A\mathbf{X} \rVert_F = \mathbf{X}' \mathbf{X}^+ \tag{4} \]

 

여기서 아래첨자 \(F\)는 프로베니우스(Frobenius)놈을 뜻하며 \(\mathbf{X}^+ \in \mathbb{R}^{(m-1) \times n}\) 는 무어-펜로즈 유사 역행렬 (Moore-Penrose pseudo inverse matrix)이다. 유사 역행렬은 특이값 분해(SVD, singular value decomposition)를 이용하면 쉽게 계산할 수 있다. 행렬 \(\mathbf{X}\) 의 SVD는 다음과 같다.

 

\[ \begin{align} \mathbf{X} &= U \Sigma V^T =[U_r \ U_{rem} ] \begin{bmatrix} \Sigma_r & 0 \\ 0 & \Sigma_{rem}\end{bmatrix} \begin{bmatrix} V_r^T \\ V_{rem}^T \end{bmatrix} \tag{5} \\ \\ & \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 \mathbb{R}^{(m-1) \times r}\) 이다. \(r\) 은 행렬 \(\mathbf{X}\) 의 차원을 \(n\) 에서 \(r\) 로 줄이기 위한 절단값(truncation value)으로서 설계 변수가 된다.

 

 

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

 

\[ \mathbf{X}^+ \approx V_r \Sigma_r^{-1} U_r^T \tag{6} \]

 

식 (6)을 (4)에 대입하면 \(A\) 는 다음과 같이 된다.

 

\[ A \approx \mathbf{X}' V_r \Sigma_r^{-1} U_r^T \tag{7} \]

 

물리적 시스템을 모델링하면 일반적으로 고차의 동적 모델이 생성된다. 식 (7)의 경우도 \(n\) 차 모델인데 계산상의 이유 때문에 이러한 모델을 차수가 감소된 단순한 모델로 근사화 하는 것이 필요할 때가 있다. 물론 이 과정에서 원래 고차 모델의 중요한 동적 특성을 포착할 수 있도록 축소 모델 (ROM, reduced order model)을 구하는 것이 중요하다.

 

 

이제 \(n\) 차 모델을 \(r\) 차 모델로 축소하도록 한다. 여기서 \(r≪n\) 이다. 모델 축소화에서 일반적으로 쓰이는 방법은 투사(projection) 방법이다. 투사 방법은 선형 기저 변환(basis transformation)을 통해서 \(n\) 차 벡터 \(\mathbf{x}_k\) 를 \(r\) 차 벡터 \(\tilde{\mathbf{x}}_k\) 로 축소 변환하는 것이다.

 

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

 

여기서 \(P \in \mathbb{R}^{n \times r}\) 는 변환 행렬로서 \(P^T P=I\) 를 만족한다. 만약 \(P=U_r\) 로 놓으면 스냅샷 POD (proper orthogonal decomposition) 문제가 된다. 식 (8)을 식 (1)과 (7)에 대입하면 다음과 같이 된다.

 

\[ \begin{align} \tilde{\mathbf{x}}_{k+1} &= U_r^T \mathbf{x}_{k+1} =U_r^T AU_r \tilde{\mathbf{x}}_k \tag{9} \\ \\ & \approx U_r^T \mathbf{X}' V_r \Sigma_r^{-1} U_r^T U_r \tilde{\mathbf{x}}_k \\ \\ &= U_r^T \mathbf{X}' V_r \Sigma_r^{-1} \tilde{\mathbf{x}}_k \\ \\ &= \tilde{A} \tilde{\mathbf{x}}_k \end{align} \]

 

여기서

 

\[ \tilde{A} = U_r^T \mathbf{X}' V_r \Sigma_r^{-1} \ \ \in \mathbb{R}^{r \times r} \tag{10} \]

 

이다.

시스템의 동특성은 시스템 행렬 \(A\) 의 고유값과 고유벡터를 계산하면 알 수 있다. 이에 앞서 우선 축소 시스템 행렬 \(\tilde{A}\) 의 고유값을 \(\lambda_i\), 고유벡터 \(\mathbf{w}_i\) 를 먼저 계산한다. 축소 시스템 행렬 \(\tilde{A}\) 의 고유값과 고유벡터 식은 다음과 같다.

 

\[ \begin{align} \tilde{A} \mathbf{w}_i &= U_r^T \mathbf{X}' V_r \Sigma_r^{-1} \mathbf{w}_i \tag{11} \\ \\ &= \lambda_i \mathbf{w}_i , \ \ \ i=1, \ ... , \ r \end{align} \]

 

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

 

\[ \ \phi_i= \mathbf{X}' V_r \Sigma_r^{-1} \mathbf{w}_i \tag{12} \]

 

그러면 식 (11)은 다음과 같이 된다.

 

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

 

위 식의 양변에 \(\mathbf{X}' V_r \Sigma_r^{-1}\) 을 곱하면,

 

\[ \ \mathbf{X}' V_r \Sigma_r^{-1} U_r^T \phi_i= \lambda_i \mathbf{X}' V_r \Sigma_r^{-1} \mathbf{w}_i \tag{14} \]

 

이 된다. 식 (7)과 (12)에 의하면 위 식은 다음과 같이 된다.

 

\[ \ A \phi_i= \lambda_i \phi_i , \ \ \ i=1, \ ..., \ r \tag{15} \]

 

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

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

      1. 먼저 \(m\) 개의 스냅샷 데이터 \(\mathbf{x}_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)} \end{align} \]

 

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

\[ \begin{align} \mathbf{X} &= U \Sigma V^T =[U_r \ U_{rem} ] \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. 축소 시스템을 구한다.

\[ \begin{align} &\tilde{\mathbf{x}}_{k+1}= \tilde{A} \tilde{\mathbf{x}}_k \\ & \tilde{A}=U_r^T \mathbf{X}' V_r \Sigma_r^{-1}, \ \ \ \mathbf{x}_k=U_r \tilde{\mathbf{x}}_k \end{align} \]

 

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

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

 

      5. DMD 모드를 계산한다.

\[ \phi_i= \mathbf{X}' V_r \Sigma_r^{-1} \mathbf{w}_i \]

 

 

 

댓글