{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1eb0b409",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from warnings import filterwarnings\n",
    "filterwarnings('ignore')\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "422d8c84",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## AERSP 597 - Machine Learning in Aerosapce Engineering\n",
    "## Lecture 22E, Dynamical System: DMD Revisited\n",
    "### Instructor: Daning Huang"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0068a128",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "$$\n",
    "\\newcommand\\norm[1]{{\\left\\lVert #1 \\right\\rVert}}\n",
    "\\newcommand{\\bR}{\\mathbb{R}}\n",
    "\\newcommand{\\vf}{\\mathbf{f}}\n",
    "\\newcommand{\\vg}{\\mathbf{g}}\n",
    "\\newcommand{\\vv}{\\mathbf{v}}\n",
    "\\newcommand{\\vw}{\\mathbf{w}}\n",
    "\\newcommand{\\vx}{\\mathbf{x}}\n",
    "\\newcommand{\\vy}{\\mathbf{y}}\n",
    "\\newcommand{\\vA}{\\mathbf{A}}\n",
    "\\newcommand{\\vB}{\\mathbf{B}}\n",
    "\\newcommand{\\vU}{\\mathbf{U}}\n",
    "\\newcommand{\\vV}{\\mathbf{V}}\n",
    "\\newcommand{\\vW}{\\mathbf{W}}\n",
    "\\newcommand{\\vX}{\\mathbf{X}}\n",
    "\\newcommand{\\vtf}{\\boldsymbol{\\phi}}\n",
    "\\newcommand{\\vty}{\\boldsymbol{\\psi}}\n",
    "\\newcommand{\\vtF}{\\boldsymbol{\\Phi}}\n",
    "\\newcommand{\\vtS}{\\boldsymbol{\\Sigma}}\n",
    "\\newcommand{\\vtL}{\\boldsymbol{\\Lambda}}\n",
    "\\newcommand{\\cF}{\\mathcal{F}}\n",
    "\\newcommand{\\cK}{\\mathcal{K}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2f8d40f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## TODAY: Dynamical System - II\n",
    "\n",
    "* More on Dynamic mode decomposition\n",
    "\n",
    "### References:\n",
    "\n",
    "+ [JFM2021](https://arxiv.org/pdf/2010.02181.pdf) Data-driven resolvent analysis, 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e63934f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Make the most out of DMD\n",
    "> What if we apply resolvent analysis to the linear system identified by DMD?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c29ae6db",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Extend the algorithm\n",
    "The first two steps remain the same\n",
    "+ **Step 1**: Rank-r truncation\n",
    "$$\n",
    "X=U_r\\Sigma_r V_r^H\n",
    "$$\n",
    "\n",
    "+ **Step 2**: System matrix - we make the distinction that $A_d=\\exp(A\\Delta t)$ is for discrete-time system\n",
    "$$\n",
    "A_d = X^1 V_r\\Sigma_r^{-1}U_r^H \\equiv BU_r^H\n",
    "$$\n",
    "and the reduced system matrix\n",
    "$$\n",
    "\\tilde{A}_d = U_r^H A_d U_r = U_r^H B\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33a0383a",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For next step, let's now consider the more general Jordan form, i.e., when $\\tilde{A}_d$ is not necessarily diagonalizable,\n",
    "+ **Step 3.1**: Reduced eigenvalue problem\n",
    "$$\n",
    "\\tilde{A}_d W = WJ_d\n",
    "$$\n",
    "where $W$ contains the eigenvectors and *generalized* eigenvectors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a2f346c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "+ **Step 3.2**: There is also the adjoint form - let $W^{-1}=S^H$,\n",
    "$$\n",
    "S^H\\tilde{A}_d = J_dS^H\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2af6e786",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using exact DMD formulation,\n",
    "+ **Step 4.1**: Eigensystem recovery\n",
    "$$\n",
    "\\begin{align}\n",
    "U_r^H B W &= WJ_d \\\\\n",
    "BU_r^H (BW) &= (BW)J_d \\\\\n",
    "A_d (BW) &= (BW)J_d\n",
    "\\end{align}\n",
    "$$\n",
    "Again one can define $\\Phi = (BW)J_d^{-1}$ as the (primal) DMD modes, so that\n",
    "$$\n",
    "A_d\\Phi = \\Phi J_d\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cdfea94",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "+ **Step 4.2**: We can also do the recovery in adjoint form\n",
    "$$\n",
    "\\begin{align}\n",
    "S^H U_r^H B &= J_dS^H \\\\\n",
    "S^H U_r^H B U_r^H &= J_dS^HU_r^H \\\\\n",
    "(U_r S)^H A_d &= J_d (U_r S)^H\n",
    "\\end{align}\n",
    "$$\n",
    "Now one can define $\\Psi = (U_r S)$ as the adjoint DMD modes, so that\n",
    "$$\n",
    "\\Psi^H A_d = J_d\\Psi^H\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1949913e",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "One can verify that the primal and adjoint DMD modes are orthogonal to each other\n",
    "$$\n",
    "\\begin{align}\n",
    "\\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\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d4bf7b",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "+ **Step 5**: Conversion to the continuous-time system\n",
    "  - The modes $\\Psi$ and $\\Phi$ remain the same\n",
    "  - The \"eigenvalues\" $J=\\frac{1}{\\Delta t}\\log(J_d)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde0f422",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Connect DMD to Resolvent\n",
    "\n",
    "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.\n",
    "\n",
    "We consider the following system with arbitrary forcing\n",
    "$$\n",
    "\\dot{x} = Ax+f\n",
    "$$\n",
    "with DMD modes $\\Phi$ and adjoint modes $\\Psi$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74f0905e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Understanding that the dynamics is happening in a reduced subspace spanned by the DMD modes, represent the states and forcing as\n",
    "$$\n",
    "x=\\Phi a,\\quad f=\\Phi b\n",
    "$$\n",
    "Then\n",
    "$$\n",
    "\\Phi \\dot{a} = A\\Phi a + \\Phi b = \\Phi J a + \\Phi b\n",
    "$$\n",
    "and left-multiply by $\\Psi^H$\n",
    "$$\n",
    "\\dot{a} = J a + b\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2bac86f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As usual, for a harmonic forcing\n",
    "$$\n",
    "b = \\hat{b}\\exp(i\\omega t)\n",
    "$$\n",
    "and by Fourier Transform\n",
    "$$\n",
    "\\hat{a} = (i\\omega I-J)^{-1} \\hat{b}\n",
    "$$\n",
    "In the full-order space\n",
    "$$\n",
    "\\hat{x}=\\Phi\\hat{a},\\quad \\hat{f}=\\Phi\\hat{b}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "911b33d7",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that\n",
    "$$\n",
    "\\norm{\\Phi a} = \\sqrt{a^H\\Phi^H\\Phi a} = \\sqrt{a^H F^HF a} = \\norm{F a}\n",
    "$$\n",
    "where we used Cholesky decomposition $\\Phi^H\\Phi=F^HF$; $F$ is upper-triangular and invertible."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c51e4ab9",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We want to maximize the gain\n",
    "$$\n",
    "G = \\max_{\\hat{f}} \\frac{\\norm{\\hat{x}}}{\\norm{\\hat{f}}}\n",
    "$$\n",
    "The ratio of norms can be simplified\n",
    "$$\n",
    "\\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}}}\n",
    "$$\n",
    "Then\n",
    "$$\n",
    "G = \\max_{F\\hat{b}} \\frac{\\norm{F(i\\omega I-J)^{-1}F^{-1}(F\\hat{b})}}{\\norm{F\\hat{b}}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40a51dd9",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Should not be surprising now, the optimal solution is given by the SVD\n",
    "$$\n",
    "F(i\\omega I-J)^{-1}F^{-1} = U(\\omega)\\Sigma(\\omega)V(\\omega)^H\n",
    "$$\n",
    "by taking the maximum singular value $\\sigma_1(\\omega)$ and the associated singular vectors $u_1(\\omega)$ and $v_1(\\omega)$\n",
    "\n",
    "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\n",
    "$$\n",
    "V_f(\\omega) = \\Phi F^{-1} V(\\omega),\\quad U_f(\\omega) = \\Phi F^{-1} U(\\omega)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "198b9e4e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example\n",
    "Switching to PPT slides..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5c6562a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Last comments\n",
    "Nevertheless, the accurate characterization of the input and output modes using data-driven DMD really relies on the data quality\n",
    "- Does it cover enough dynamics (e.g. start from different initial conditions)?\n",
    "- If highly non-normal, it is possible that DMD may be inaccurate due to numerical issues\n",
    "  + So it may **fail**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b42cc5",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For highly non-normal cases\n",
    "- If one already has the nonlinear solver (code) and its linearization\n",
    "  + Directly do resolvent analysis for modal/nonmodal analysis\n",
    "  + Use e.g. **Balanced Proper Orthogonal Decomposition** [Rowley2005](https://www.worldscientific.com/doi/abs/10.1142/S0218127405012429) to construct a linear system that captures the input/output modes\n",
    "- If only data is available\n",
    "  + Try something as simple as **Eigensystem Realization Algorithm** [Juang1985](https://arc.aiaa.org/doi/abs/10.2514/3.20031) - but it does **not** produce physically meaningful states or input modes\n",
    "  + If there are many data snapshots, try also the **Spectral Proper Orthogonal Decomposition** [Towne2018](https://arxiv.org/pdf/1708.04393.pdf) - but it only identifies the output modes"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}