Machine Learning for Engineering

Dynamical System: DMD Revisited

Instructor: Daning Huang

$$ \newcommand\norm[1]{{\left\lVert #1 \right\rVert}} \newcommand{\bR}{\mathbb{R}} \newcommand{\vf}{\mathbf{f}} \newcommand{\vg}{\mathbf{g}} \newcommand{\vv}{\mathbf{v}} \newcommand{\vw}{\mathbf{w}} \newcommand{\vx}{\mathbf{x}} \newcommand{\vy}{\mathbf{y}} \newcommand{\vA}{\mathbf{A}} \newcommand{\vB}{\mathbf{B}} \newcommand{\vU}{\mathbf{U}} \newcommand{\vV}{\mathbf{V}} \newcommand{\vW}{\mathbf{W}} \newcommand{\vX}{\mathbf{X}} \newcommand{\vtf}{\boldsymbol{\phi}} \newcommand{\vty}{\boldsymbol{\psi}} \newcommand{\vtF}{\boldsymbol{\Phi}} \newcommand{\vtS}{\boldsymbol{\Sigma}} \newcommand{\vtL}{\boldsymbol{\Lambda}} \newcommand{\cF}{\mathcal{F}} \newcommand{\cK}{\mathcal{K}} $$

TODAY: Dynamical System

• More on Dynamic mode decomposition

References:

• JFM2021 Data-driven resolvent analysis, 2021

Make the most out of DMD

What if we apply resolvent analysis to the linear system identified by DMD?

Extend the algorithm

The first two steps remain the same

• Step 1: Rank-r truncation $$ X=U_r\Sigma_r V_r^H $$

• Step 2: System matrix - we make the distinction that $A_d=\exp(A\Delta t)$ is for discrete-time system $$ A_d = X^1 V_r\Sigma_r^{-1}U_r^H \equiv BU_r^H $$ and the reduced system matrix $$ \tilde{A}_d = U_r^H A_d U_r = U_r^H B $$

For next step, let's now consider the more general Jordan form, i.e., when $\tilde{A}_d$ is not necessarily diagonalizable,

• Step 3.1: Reduced eigenvalue problem $$ \tilde{A}_d W = WJ_d $$ where $W$ contains the eigenvectors and generalized eigenvectors.

• Step 3.2: There is also the adjoint form - let $W^{-1}=S^H$, $$ S^H\tilde{A}_d = J_dS^H $$

Using exact DMD formulation,

• Step 4.1: Eigensystem recovery $$ \begin{align} U_r^H B W &= WJ_d \\ BU_r^H (BW) &= (BW)J_d \\ A_d (BW) &= (BW)J_d \end{align} $$ Again one can define $\Phi = (BW)J_d^{-1}$ as the (primal) DMD modes, so that $$ A_d\Phi = \Phi J_d $$

• Step 4.2: We can also do the recovery in adjoint form $$ \begin{align} S^H U_r^H B &= J_dS^H \\ S^H U_r^H B U_r^H &= J_dS^HU_r^H \\ (U_r S)^H A_d &= J_d (U_r S)^H \end{align} $$ Now one can define $\Psi = (U_r S)$ as the adjoint DMD modes, so that $$ \Psi^H A_d = J_d\Psi^H $$

One can verify that the primal and adjoint DMD modes are orthogonal to each other $$ \begin{align} \Psi^H\Phi = (U_rS)^H (BW) J_d^{-1} = S^H \underbrace{U_r^H B}_{\tilde{A}_d} W J_d^{-1} = J_d \underbrace{S^H W}_{I} J_d^{-1} = I \end{align} $$

• Step 5: Conversion to the continuous-time system
  • The modes $\Psi$ and $\Phi$ remain the same
  • The "eigenvalues" $J=\frac{1}{\Delta t}\log(J_d)$

Connect DMD to Resolvent

Now that one obtains a linear system from DMD, it is natural to apply resolvent analysis to extract out the gain characteristics over the spectrum. Due to the special structure in DMD operator, the resolvent can be computed efficiently.

We consider the following system with arbitrary forcing $$ \dot{x} = Ax+f $$ with DMD modes $\Phi$ and adjoint modes $\Psi$.

Understanding that the dynamics is happening in a reduced subspace spanned by the DMD modes, represent the states and forcing as $$ x=\Phi a,\quad f=\Phi b $$ Then $$ \Phi \dot{a} = A\Phi a + \Phi b = \Phi J a + \Phi b $$ and left-multiply by $\Psi^H$ $$ \dot{a} = J a + b $$

As usual, for a harmonic forcing $$ b = \hat{b}\exp(i\omega t) $$ and by Fourier Transform $$ \hat{a} = (i\omega I-J)^{-1} \hat{b} $$ In the full-order space $$ \hat{x}=\Phi\hat{a},\quad \hat{f}=\Phi\hat{b} $$

Note that $$ \norm{\Phi a} = \sqrt{a^H\Phi^H\Phi a} = \sqrt{a^H F^HF a} = \norm{F a} $$ where we used Cholesky decomposition $\Phi^H\Phi=F^HF$; $F$ is upper-triangular and invertible.

We want to maximize the gain $$ G = \max_{\hat{f}} \frac{\norm{\hat{x}}}{\norm{\hat{f}}} $$ The ratio of norms can be simplified $$ \frac{\norm{\hat{x}}}{\norm{\hat{f}}} = \frac{\norm{\Phi\hat{a}}}{\norm{\Phi\hat{b}}} = \frac{\norm{F\hat{a}}}{\norm{F\hat{b}}} = \frac{\norm{F(i\omega I-J)^{-1}\hat{b}}}{\norm{F\hat{b}}} = \frac{\norm{F(i\omega I-J)^{-1}F^{-1}(F\hat{b})}}{\norm{F\hat{b}}} $$ Then $$ G = \max_{F\hat{b}} \frac{\norm{F(i\omega I-J)^{-1}F^{-1}(F\hat{b})}}{\norm{F\hat{b}}} $$

Should not be surprising now, the optimal solution is given by the SVD $$ F(i\omega I-J)^{-1}F^{-1} = U(\omega)\Sigma(\omega)V(\omega)^H $$ by taking the maximum singular value $\sigma_1(\omega)$ and the associated singular vectors $u_1(\omega)$ and $v_1(\omega)$

But be careful as we optimized over $F\hat{b}$, we need to convert back to the reduced-order system $\hat{b}$, then the full-order system $\hat{f}$. The input and output modes are, respectively $$ V_f(\omega) = \Phi F^{-1} V(\omega),\quad U_f(\omega) = \Phi F^{-1} U(\omega) $$

Last comments

Nevertheless, the accurate characterization of the input and output modes using data-driven DMD really relies on the data quality

• Does it cover enough dynamics (e.g. start from different initial conditions)?
• If highly non-normal, it is possible that DMD may be inaccurate due to numerical issues
  • So it may fail

For highly non-normal cases

• If one already has the nonlinear solver (code) and its linearization
  • Directly do resolvent analysis for modal/nonmodal analysis
  • Use e.g. Balanced Proper Orthogonal Decomposition Rowley2005 to construct a linear system that captures the input/output modes
• If only data is available
  • Try something as simple as Eigensystem Realization Algorithm Juang1985 - but it does not produce physically meaningful states or input modes