Reduced-order control using low-rank Dynamic Mode Decomposition

TL;DR

This paper introduces low-rank Dynamic Mode Decomposition (lrDMD), a data-driven model reduction technique that improves control of fluid flows by providing lower reconstruction error and better stabilization with reduced model complexity.

Contribution

The paper proposes lrDMD, a novel low-rank constrained DMD method that enhances model accuracy and control performance over standard DMD and OMD methods.

Findings

lrDMD achieves lower data reconstruction error than DMD and OMD.

lrDMD stabilizes complex fluid systems with lower-order controllers.

lrDMD provides better estimates of the adjoint system at a given rank.

Abstract

In this work we perform full-state LQR feedback control of fluid flows using non-intrusive data-driven reduced-order models. We propose a model reduction method called low-rank Dynamic Mode Decomposition (lrDMD) that solves for a rank-constrained linear representation of the dynamical system. lrDMD is shown to have lower data reconstruction error compared to standard Optimal Mode Decomposition (OMD) and Dynamic Mode Decomposition (DMD), but with an increased computational cost arising from solving a non-convex matrix optimization problem. We demonstrate model order reduction on the complex linearized Ginzburg-Landau equation in the globally unstable regime and on the unsteady flow over a flat plate at a high angle of attack. In both cases, low-dimensional full-state feedback controller is constructed using reduced-order models constructed using DMD, OMD and lrDMD. It is shown that lrDMD stabilizes the Ginzburg-Landau system with a lower order controller and is able to suppress vortex shedding from an inclined flat plate at a cost lower than either DMD or OMD. It is further shown that lrDMD yields an improved estimate of the adjoint system, for a given rank, relative to DMD and OMD.

Tables

Table 1. Table 1: Time taken in seconds by DMD, OMD and lrDMD for generating reduced-order models of rank r 𝑟 r with data matrices composed of n 𝑛 n snapshots. lrDMD (subProj) refers to subspace projection method implemented using Algorithm  1 and lrDMD (GradD) refers to gradient based method described in detail in  sashittal2018low

Method	$n=50$			$n=200$
	$r=10$	$r=20$	$r=30$	$r=10$	$r=20$	$r=30$
DMD	$0.1443$	$0.1317$	$0.1096$	$0.6965$	$0.6691$	$0.6464$
OMD	$3.2687$	$5.5517$	$10.5841$	$1.6759$	$17.4687$	$26.5706$
lrDMD (subProj)	$0.7731$	$0.4127$	$1.4995$	$2.3426$	$2.8213$	$1.4139$
lrDMD (GradD)	$19.9936$	$2.8781$	$0.8261$	$58.0509$	$17.1708$	$26.0866$

Table 2. Table 2: Error norm ϵ italic-ϵ \epsilon of reconstruction for reduced-order models of rank r 𝑟 r generated by DMD, OMD and lrDMD with data matrices composed of n 𝑛 n snapshots. lrDMD (subProj) refers to subspace projection method implemented using Algorithm  1 and lrDMD (GradD) refers to gradient based method described in detail in  sashittal2018low .

Method	$n=50$			$n=200$
	$r=10$	$r=20$	$r=30$	$r=10$	$r=20$	$r=30$
DMD	$14.7449$	$5.7669$	$3.3224$	$24.0859$	$2.2424$	$0.3296$
OMD	$4.2816$	$0.0921$	$0.0049$	$12.1151$	$0.6408$	$0.0176$
lrDMD (subProj)	$1.8183$	$0.0160$	$0.0015$	$6.3463$	$0.5206$	$0.0295$
lrDMD (GradD)	$1.2686$	$0.0031$	$0.0015$	$5.0654$	$0.1638$	$0.0058$

Table 3. Table 3: Ginzburg-Landau equation parameter values corresponding to the supercritical, globally unstable regime case in  chen2011h .

variable	description	value
$U$	advection velocity	$2$
$c_u$	most unstable wavenumber	$0.2$
$c_d$	dispersion parameter	$-1.0$
${\mu }_0$	overall amplification	$0.41$
${\mu }_2$	degree of non-parallelism	$-0.01$
${\mu }_t$	transitional $\mu$	$0.32$
${\mu }_c$	critical $\mu$	$0.4$

Table 4. Table 4: Numerical parameters for the actuator location and strength for feedback control of flow past inclined flat plate.

$x_{c,1}$	$y_{c,1}$	$x_{c,2}$	$y_{c,2}$	$a$	$c$
0	1.3423	0	0.89	20	2

Table 5. Table 5: Summary of the feedback control performance of DMD, OMD and lrDMD based controllers for flow over inclined flat plate. The time at which the controller was activated is indicated by t 0 subscript 𝑡 0 t_{0} and T = { 50 , 100 } 𝑇 50 100 T=\{50,100\} refers to the time separation of corresponding data snapshots in the data matrix used to construct the controllers. ✓  indicates that the flow converged to steady state and ✗  indicates that the flow became unstable.

Method	$t_0=170$		$t_0=190$		$t_0=210$
	$T=50$	$T=100$	$T=50$	$T=100$	$T=50$	$T=100$
DMD	✓	✓	✓	✗	✗	✗
OMD	✓	✓	✓	✓	✓	✗
lrDMD	✓	✓	✓	✓	✓	✗

Equations

$$x_{k+1}=f(x_k)$$

$$X:=(x_1|\cdots |x_n),Y:=(x_2|\cdots |x_{n+1}).$$

$$\hat{A}=\underset{A\in {\mathbb{R}}^{m\times m}}{\mathrm{argmin}}{\Vert Y-AX\Vert }_F^2,$$

$$X\approx U\Sigma V^T,$$

$$\tilde{A}=U^{T}YV{\Sigma }^{-1}.$$

$$\hat{A}=U\tilde{A}{U}^T.$$

$$\underset{L,M}{\min }$$

$$L^{T}L=I$$

$$M\in {\mathbb{R}}^{r\times r},L\in {\mathbb{R}}^{m\times r}.$$

$$M^{\ast }(L)=L^{T}Y{X}^{T}L(L^{T}X{X}^{T}L{)}^{-1}.$$

$$\underset{L,D,R}{\min }{\Vert Y-LD{R}^{T}X\Vert }_F^2,$$

$$D^{\ast }(L,R)=(L^{T}Y{X}^{T}R)(R^{T}X{X}^{T}R{)}^{-1}.$$

$$\underset{L,R}{\min }\left(-{\Vert L^{T}Y{Q}_R\Vert }_F^2\right)$$

$$\mathcal{M}:={\mathcal{G}}_{r,m}\times {\mathcal{G}}_{r,m}.$$

$$\underset{L,R}{\min }$$

$$(L,R)\in \mathcal{M}:={\mathcal{G}}_{r,m}\times {\mathcal{G}}_{r,m}$$

$$Q_R=X^{T}R(R^{T}X{X}^{T}R{)}^{-1}{R}^{T}X.$$

$$M=U\Sigma V^T.$$

$$L$$

$$R$$

$$\underset{\underline{L},\underline{R}}{\min }$$

$$(\underline{L},\underline{R})\in \mathcal{M}:={\mathcal{G}}_{r,p}\times {\mathcal{G}}_{r,p}$$

$$Q_{\underline{R}}={\underline{X}}^T\underline{R}({\underline{R}}^T\underline{X}{\underline{X}}^T\underline{R}{)}^{-1}{\underline{R}}^T\underline{X}$$

$$J(q,u)=$$

$$q_{k+1}=A{q}_k+B{u}_k,$$

$$A^{T}PA-P-(A^{T}PB)(S+B^{T}PB{)}^{-1}(B^{T}PA)+Q=0,$$

$$E(P)=A^{T}PA-P-(A^{T}PB)(S+B^{T}PB{)}^{-1}(B^{T}PA)+Q.$$

$$\bar{P}=WP{W}^T,$$

$$\frac{\partial q}{\partial t}+\nu \frac{\partial q}{\partial x}=\mu (x)q+\gamma \frac{{\partial }^{2}q}{\partial x^2},$$

$$\mu (x)={\mu }_0-c_u^2+{\mu }_2{x}^2/2,$$

$$q_{k+1}=A{q}_k+B{u}_k,$$

$$J=\sum _{k=1}^{\infty }(q_k^{\ast }Q{q}_k+u_k^{\ast }S{u}_k),$$

$$B=\exp \left(-\frac{(x-x_a{)}^2}{2{\sigma }^2}\right),$$

$$ϵ(q,V)=\frac{{\Vert q-{\mathbb{P}}_{V}q\Vert }_2}{{\Vert q\Vert }_2}$$

$$J=$$

Full text

1 Introduction

The dynamics of many flows of engineering interest are high dimensional. Solving optimal control and design problems for these high dimensional systems is challenging, even on very large computers kim2014adjoint . This creates a need for reduced-order models that can faithfully approximate the dynamics of the system. Model reduction and reduced-order control therefore has been an active area of research in fluid mechanics. The reader is referred to rowley2017model for a thorough review of model reduction techniques for flow analysis and control.

Data-driven model reduction techniques use high-fidelity numerical or experimental data to extract simplified models that capture essential features of the flow. Model-based control methods, such as Galerkin projection carlberg2011efficient ; willcox2002balanced ; rowley2005model , extract a reduced basis from the snapshots and project the full system dynamics onto the linear subspace formed by the reduced basis. Reduced order feedback control is then constructed using the resulting reduced order model. Model-based feedback control has been used in bagheri2009input ; semeraro2011feedback ; semeraro2013transition to control flat-plate boundary layers and in illingworth2016model to suppress vortex shedding over a circular cylinder. Although these methods can be very effective (see kim2007linear for a thorough review), they require prior knowledge of the full-state governing equations and require computation of the projection coefficients, which is inconvenient for complex nonlinear systems.

To overcome these issues, non-intrusive system identification methods are used to construct reduced order models using data without prior knowledge of the governing equations. The authors in huang2008control and herve2012physics used Autoregressive Modeling akaike1969fitting applied to flow over an airfoil and backward facing step respectively. In this method, a mathematical relationship between the input and output data is assumed with unknown parameters that are learned by minimizing the error between the data and model prediction. An alternative approach is to estimate the state-space model using ‘subspace identification methods’ such as N4SID van1994n4sid ; van2012subspace which have been used for control of transitional boundary layer in inigo2014dynamic . A thorough review of these methods is presented in qin2006overview . Finally, Eigensystem Realization Algorithm (ERA) juang1985eigensystem , which gives a balanced reduced order state-space model of the system, has been used extensively in flow control applications brunton2014state ; flinois2016feedback ; ahuja2010feedback ; illingworth2012feedback ; belson2013feedback ; illingworth2014active .

Dynamic Mode Decomposition (DMD) schmid2010dynamic ; tu2013dynamic is a data-driven model order reduction method that learns a linear approximation of the dynamical system such that the model prediction best fits the data snapshots in the $L_{2}$ norm. DMD and its variants goulart2012optimal ; chen2011h ; hemati2014dynamic ; hemati2017biasing ; dawson2016characterizing , have also been used for reduced-order control of systems with uncertain parameters kramer2017feedback , for suppressing flow separation deem2018experimental ; hemati2016improving , and for optimal actuator selection for airfoil separation control bhattacharjee2018optimal . These methods approximate the dynamics of the system by fitting an endomorphic linear function on some low-dimensional subspace. The linear approximation of the dynamical system is learned constrained to a particular low dimensional subspace. This is very useful in model order reduction since it provides a low order basis to describe the flow. However, in this study, we find that this constraint is restrictive for the reduced order control performance of these methods. We investigate the possibility of using a linear map between different subspaces to construct reduced-order models and controllers for unsteady fluid flows. To do this, we formulate a rank-constrained matrix optimization problem and propose two methods to solve it. We call this model reduction method low-rank Dynamic Mode Decomposition (lrDMD) sashittal2018low ; heas2016low . We apply full-state feedback control, constructed using the LQR framework, on model flow problems such as the Ginzburg-Landau equation and flow past an inclined flat plate. The performance of the controller is compared to controllers constructed using DMD and Optimal Mode Decomposition (OMD) goulart2012optimal ; wynn2013optimal methods.

The outline of the paper is as follows. Section 2 outlines the mathematical formulation of the problem and the three methods we investigate: DMD, OMD and lrDMD. Section 3 details the numerical implementation of subspace projection and gradient descent algorithms used to solve the lrDMD optimization problem, including the model reduction computational performance while Section 4 describes the controller design. The control application results are shown in Section 5 and Section 6 concludes the paper.

2 Mathematical Formulation

Consider a dynamical system with the state vector $x\in\mathbb{R}^{m}$ such that

[TABLE]

where the subscript denote time iteration number. We have access to a sequence of time snapshots of this state vector which we represent as a data matrices $X,Y\in\mathbb{R}^{m\times n}$ formed by $n$ pairs of data snapshots as follows,

[TABLE]

The matrix $Y$ is called the time-shifted data matrix associated with $X$. Our aim is to build a reduced-order model of the dynamical system (1) that captures the essential features of the function $f(\cdot)$. To this end, we construct a linear approximation of $f(\cdot)$ from the given data matrices. A natural choice of optimization problem to solve for this purpose would be

[TABLE]

where $\left\lVert\cdot\right\rVert_{F}$ is the Frobenius norm, $\hat{A}\in\mathbb{R}^{m\times m}$ is the inferred state transition matrix. However, for fluid systems $m$ is generally large, frequently in the range of $10^{6}-10^{9}$, which makes this optimization problem computationally intractable. To remedy this, low rank versions of this problem are adopted. We will review two popular methods: Dynamic Mode Decomposition (DMD) schmid2010dynamic ; tu2013dynamic and Optimal Mode Decomposition (OMD) goulart2012optimal ; wynn2013optimal , and introduce a generalized version referred to as low-rank Dynamic Mode Decomposition (lrDMD) sashittal2018low ; heas2016low .

2.1 Dynamic Mode Decomposition

DMD schmid2010dynamic ; tu2013dynamic provides the optimal linear representation of the dynamical system in the space spanned by the Proper Orthogonal Decomposition (POD) modes of the data matrix $X$. POD modes are the directions that maximally capture the energy or variance in a given data matrix. For a deterministic system, and taking the $L_{2}$ norm of the state vector as the energy measure, the POD modes are given by the left singular vectors of the data matrix $X$ ranked by the value of the corresponding singular values. For a $r$-rank representation of the dynamics we consider the space spanned by the $r$ leading POD modes of the data matrix. Let the rank $r$ reduced singular value decomposition of $X$ be

[TABLE]

where $U\in\mathbb{R}^{m\times r}$, $\Sigma\in\mathbb{R}^{r\times r}$ and $V\in\mathbb{R}^{n\times r}$. Then the DMD solution $\tilde{A}\in\mathbb{R}^{r\times r}$ is given by

[TABLE]

The rank $r$ linear approximation of the state transition matrix is given by

[TABLE]

2.2 Optimal Mode Decomposition

In OMD goulart2012optimal ; wynn2013optimal , the following optimization problem is solved to give an approximate solution of (2),

[TABLE]

where $r$ is the rank of the linear approximation of the state transition matrix $\hat{A}=LML^{T}$. This particular form of low rank approximation has useful implications in reduced-order modeling. The matrix $L$ provides the basis of a low dimensional subspace to approximate the trajectories of the dynamical system while $M$ provides the dynamical evolution of the state vector on this subspace. The solution of $M$ for a fixed $L$ is given by

[TABLE]

The OMD solution coincides with the DMD approximation when $L$ is chosen as the $r$ leading left singular vectors of the data matrix $X$. In this way OMD is a generalization of the DMD formulation.

2.3 Low-Rank Dynamic Mode Decomposition

We consider a rank-constrained solution of the optimization problem in (2). Recently the authors of heas2016low showed that there exists a closed form solution for this optimization problem. We solve the following equivalent optimization problem,

[TABLE]

where $L,R\in\mathbb{R}^{m\times r}$ and $L^{T}L=R^{T}R=I_{r}$ ($r\times r$ identity matrix) and $\hat{A}=LDR^{T}$ is the $r$-ranked matrix approximating the dynamics of the underlying system. Observe that (6) optimizes over all possible left and right subspaces of the state transition matrix in contrast to standard DMD and OMD. A key step to solving (6) is to find the optimal solution for $D$ as a function of $L$ and $R$. For a fixed $L$ and $R$ there exists a closed form optimal solution for $D$, namely

[TABLE]

Using (7) in the objective function and simplifying further the optimization problem can be reduced to

[TABLE]

where $Q_{R}=X^{T}R(R^{T}XX^{T}R)^{-1}R^{T}X$. We will denote the objective function (8) by $G(L,R)$. An important observation is that $G(L,R)$ is only a function of the spaces spanned by the columns of $L$ and $R$. To see this consider an orthogonal matrix $P\in\mathbb{R}^{r\times r}$ such that $P^{T}P=I_{r}$. For this objective function, $G(LP,R)=G(L,R)$ and $G(L,RP)=G(L,R)$. This implies that the objective function does not change values as long as the space spanned by the columns of $L$ and $R$ do not change. Hence the optimization can be equivalently performed over the set of $r$-dimensional subspaces in $\mathbb{R}^{m}$ instead of the set of orthogonal matrices in $\mathbb{R}^{m\times r}$.

The set of $r$-dimensional subspaces in $\mathbb{R}^{m}$ is known as the Grassmanian manifold $\mathcal{G}_{r,m}$ absil2009optimization ; edelman1998geometry . In matrix representation an element of $\mathcal{G}_{r,m}$ is specified by an orthogonal basis of the $r$-dimensional subspace. The manifold on which we perform the optimization (8) is a product manifold of two Grassmanian manifolds

[TABLE]

An element of $\mathcal{M}$, in matrix representation, will be characterized by $M=(L,R)$.

Remark 1

The $L$ and $R$ basis solutions of the optimization problem (8) need not be the best basis for a Petrov-Galerkin reduced-order model. These bases only provide the subspaces for the best low-rank linear description of the dynamical system in the least-squares sense. As will be made clear in Section 4, for the reduced-order model we use the input basis $R$ for both the trial subspace and the test subspace.

3 Numerical Methods for lrDMD optimization problem

In this section we present the numerical methods employed to solve the lrDMD optimization problem and compare the computational performance of the methods to OMD and DMD. Recall that the objective function and the optimization problem are,

[TABLE]

Two methods to solve this optimization problem are detailed in Appendix A. The first method in Appendix A.1 is a computationally efficient subspace projection method that provides an approximate solution to the problem. The second method in Appendix A.2 is a gradient-based method that is guaranteed to converge to a local minimum.

Remark 2

There is no guarantee that the gradient based method will converge to the global minimizer of the given optimization problem (9). The gradient based method can therefore be thought as a post-processor for the initial guess provided by either DMD or the subspace projection method described in Algorithm 1.

3.1 Computational performance

The advantage to using the form of the optimization problem in Eq. (6) instead of Eq. (2) is that we only need to solve for $\mathcal{O}(m)$ variables for the fully parameterized state transition matrix. For large $m$, even this can make the problem computationally intractable. An important observation that can help is that we only need to consider the solution of $R$ such that $\text{Im}(R)\subseteq\text{Im}(X)$. Consider $M$ to be the full ranked data matrix that contains all the data snapshots that form the data matrices $X$ and $Y$ and let there be $p$ such snapshots. Consider the singular value decomposition of the data matrix $M\in\mathbb{R}^{m\times p}$

[TABLE]

We assume that the basis of both $L$ and $R$ is contained within $U$ so that

[TABLE]

where $\underline{L},\underline{R}\in\mathcal{G}_{r,p}$. The optimization problem reduces to

[TABLE]

where $\underline{Y}=U^{T}Y$ and $\underline{X}=U^{T}X$. Since $p\ll m$, this saves significant computational time. In this study, we use this technique for lrDMD and OMD where gradient based methods are used to solve the optimization problem. Therefore, in this study, OMD and the gradient based implementation of lrDMD solve the optimization problem with the reduced data matrices $\underline{Y}$ and $\underline{X}$.

Table 1 compares the computational performance of DMD, OMD and the two algorithms to solve the lrDMD problem. It shows the time taken in seconds by the three methods in computing a reduced-order model for the given data matrices. The data snapshots have dimension $m=62001$ and describe the vorticity field for a flow past flat plate described in Section 5.2 with consecutive snapshots that are $100$ timesteps apart. The data matrices are generated with $n=\{50,200\}$ snapshots and reduced-order models of rank $r=\{10,20,30\}$ are computed. lrDMD (subProj) refers to the solution of lrDMD using the subspace projection method described in Algorithm 1 whereas lrDMD (GradD) refers to the gradient based method sashittal2018low . The gradient based methods used in OMD and lrDMD are implemented using ManOpt boumal2014manopt package on MATLAB using the trust-region absil2007trust algorithm. The initial condition for OMD is given by the DMD solution whereas the initial condition for the lrDMD is the solution from Algorithm 1. All computations are performed on a standard desktop PC with a 2.2 GHz quad-core Intel i7 processor and 16GB RAM running on Mac OS X 10.11.

Table 1 shows that time taken by DMD is always lower than OMD or either method of solving the lrDMD problem, with the subspace projection method incurring computational times $2-3\times$ larger than DMD but is faster than OMD for all cases. An interesting observation is that for certain cases like $(n,r)=(50,30)$, the initial condition provided to ‘lrDMD (GradD)’ by the subspace projection method is very close to the optimal and no iterations of the gradient based method are required for convergence. In almost all cases lrDMD with gradient descent takes around the same time as OMD except a few cases when OMD takes less time than lrDMD.

Table 2 shows the error norm defined by $\epsilon=\left\lVert Y-\hat{A}X\right\rVert_{F}$ for all the cases considered in Table 1. It can be seen that lrDMD with gradient descent and subspace projection shows error norm much less than OMD or DMD for all cases considered. Subspace projection in particular, shows lesser error than OMD in almost all cases while taking lesser time to compute the reduced-order model. For data snapshots that are spaced far apart, lrDMD with gradient descent is recommended for accurate low rank reconstruction of the system dynamics if computation time is not a factor.

4 Controller Design

Once the reduced order model is learned, the next step is to construct the feedback control law using the reduced order model. We use a Linear Quadratic Regulator kwakernaak1972linear (LQR) to construct a stabilizing feedback control independent of the initial condition. Consider a linear system with the state variable $q\in\mathbb{R}^{m}$ and control input $u\in\mathbb{R}^{p}$. The solution of the LQR problem provides a feedback gain $K\in\mathbb{R}^{p\times m}$ such that the control input $u_{k}=-Kq_{k}$ minimizes

[TABLE]

where the subscripts denote the temporal iteration number, $A\in\mathbb{R}^{m\times m}$ is the state transition matrix of the linear system and $B\in\mathbb{R}^{m\times p}$ describes the actuator. Under the assumption that $(A,B)$ is stabilizable, the feedback gain matrix $K$ can be computed by solving the discrete algebraic Riccati equation

[TABLE]

for $P$ and using $K=(S+B^{T}PB)^{-1}B^{T}PA$. Although most fluid flow systems are nonlinear, the LQR is a powerful tool used for stabilization about an unstable base flow when a feedback form of controller can be employed.

Direct methods of solving the Riccati equation (10) have time complexity of $O(m^{6})$. As noted in the introduction, fluid flows are governed by partial differential equations which result in high dimensional systems upon discretizations, making solving Eq. (10) computationally infeasible for practical fluid flow applications. To remedy this, we employ a Galerkin projection based method simoncini2016analysis ; alla2017order that yields a low-rank approximation of the solution to Eq. (10). Let $E$ be the residual of the discrete Riccati equation as follows,

[TABLE]

The projection method works by restricting the solution in some subspace and imposing an orthogonality of the residual to that subspace. Consider a $k$-dimensional subspace spanned by the basis $W\in\mathbb{R}^{m\times k}$. The solution of Eq. (10) is approximated as,

[TABLE]

and the orthogonality condition on the residual is $W^{T}E(\bar{P})W=0$. For all the reduced-order models in this study, we choose $W$ to be some basis of the row space of the low-rank approximation of $A$. This is given by the POD basis in case of DMD, the optimal solution for $L$ in case of OMD and the optimal solution of $R$ in case of lrDMD. This choice permits a natural comparison between the three reduction methods.

5 Results

In this section we compare the performance of the reduced-order LQR controllers built using DMD, OMD and lrDMD on unsteady dynamical systems. The first system we look into is the complex linear Ginzburg-Landau equation. We use the reduced-order models to find optimal feedback control as well as the optimal actuator location for the system in the unstable regime. We also employ LQR controllers on the incompressible flow past inclined flat plate at high angle of attack. The results show that the reduced-order linear feedback controllers can be effective in regions of phase-space with strong nonlinearities and can be used to suppress nonlinear vortex shedding.

5.1 Linearized Ginzburg Landau Equation

In this section we apply the control strategies described so far on the linearized complex Ginzburg-Landau (GL) equation, which is a well known model equation for fluid systems. Even though all the analysis in this paper has been on the real number field, the results can be easily extended complex numbers. The GL equation is as follows,

[TABLE]

where the real part of $q(x,t)$ represents velocity or stream function perturbation amplitude. The GL equation exhibits a variety of stability behaviors observed in fluid flows for different values of the constants $\nu,\,c_{u},\,\mu_{0}$ and $\mu_{2}$ in different regions of the spatial domain. For this study we have taken the value of these constants corresponding to the supercritical, globally unstable regime case in chen2011h , shown in Table 3.

In discrete space we use the spectral formulation of the derivative operators using Hermite polynomials evaluated at the Hermite nodes bagheri2009input . Time discretization is performed using forward Euler time-stepping scheme. In this study we consider $n=220$ grid points with the computational domain $x\in[-85,85]$ and a time-step of $dt=1$. For more details about the system the reader is referred to bagheri2009input .

We will describe this system as a discrete input-output system as,

[TABLE]

where $q_{k}\in\mathbb{R}^{m}$ is the state variable, $A\in\mathbb{R}^{m\times m}$ is the state transition matrix, $B\in\mathbb{R}^{m\times p}$ is the spatial support of the controller and $u_{k}\in\mathbb{R}^{p}$ is the control input. Our goal is to find $K\in\mathbb{R}^{p\times m}$ such that when $u_{k}=-Kq_{k}$, we minimize the following cost,

[TABLE]

for given positive definite matrices $Q$ and $S$. We denote the conjugate transpose of $q$ by $q^{*}$. The controller and the values of $Q$ and $S$ are taken from chen2011h . The controller we choose is a Gaussian centered in the convectively unstable region of flow,

[TABLE]

where we choose $x_{a}=8$ and $\sigma=5$ which places the actuator just inside the region of amplification which is $[-8.6,8.6]$. Using $15$ snapshots of the impulse response of the flow, we generate a low-rank approximation of $A$ which is used to construct reduced-order controllers using the LQR framework and to find the optimal actuator location.

5.1.1 Adjoint Reconstruction

It is known that state transition matrices arising from the discretization of the linearized governing equations for various fluid flow applications are nonnormal schmid2012stability . Nonnormal systems exhibit non-orthogonal eigenmodes which can differ significantly from the adjoint modes schmid2014analysis that are known to play a major role in flow control and optimization luchini2014adjoint . Therefore, it is crucial for effective flow control applications that the reduced order models not only accurately predict the flow field but also extract adjoint information from the data.

Figure 1a shows the eigenvalues of the discretized state transition matrix of the GL equation in the supercritical regime. It shows one unstable eigenvalue $0.8073-0.6109i$ outside the unit circle. Figure 1b shows the eigenmode $v$ and the adjoint mode $w$ corresponding to this unstable eigenvalue. For faithful reconstruction of these unstable modes, the eigenmode must lie in the column space of the reduced order model $\hat{A}$ and the adjoint mode must lie in the row space. We define the projection error

[TABLE]

where $\mathbb{P}_{V}$ is the orthogonal projection matrix in the column space of $V$. We look at the projection error of the unstable eigenmode in the column space of $\hat{A}$ given by $\epsilon(v,\hat{A})$ and the projection error of the unstable adjoint mode in the row space of $\hat{A}$ given by $\epsilon(w,\hat{A}^{*})$. For the DMD-based reduced order model, bases for the row space and column space are both given by the POD modes. For OMD, the solution $L$ of the optimization problem in Eq. (3) serves as the basis of both the column and the row space of the reduced order model. In the case of lrDMD, the basis of the row space is given by $R$ and the basis of the column space is given by $L$ from the solution of the optimization problem in Eq. (8).

Figure 1c shows the projection error of the unstable eigenmode for the reduced order models of different ranks obtained from DMD, OMD and lrDMD. At all rank approximations the error incurred by OMD and lrDMD are about the same and less than the error incurred by DMD. Similarly, Figure 1d shows the projection error of the unstable adjoint mode. In this case we see that lrDMD outperforms both OMD and DMD at all ranks. This shows that the additional degree of freedom of choosing separate input and output spaces in lrDMD allow for lower projection error in the adjoint modes of interest while also keeping the projection error in the eigenmodes as low as in OMD.

5.1.2 Optimal Control

In the supercritical regime, the GL system is globally unstable, with a single eigenvalue of the state transition matrix that lies outside the unit circle, as shown in Figure 1a. The goal of the feedback controllers is to stabilize the system while minimizing the cost in Eq. (12).

Figure 2 shows the eigenvalues of the controlled Ginzburg-Landau system using reduced-order LQR controllers constructed using DMD, OMD and lrDMD. Figures 2a and 2b illustrate the rank 5 approximations whereas Fig.  2c shows the rank 9 approximations. Figure 2d shows the eigenvalues of the controlled GL system with LQR optimal control constructed using the full-order system matrix. It can be clearly seen in Fig. 2b that for the rank 5 approximation, DMD is not able to stabilize the system and has one eigenvalue outside the unit circle. OMD also fails to stabilize the system with one eigenvalue outside the unit circle. This shows that $5$ modes are not sufficient for reduced order control even with OMD which finds the optimal basis for the low order subspace. However, the freedom of choosing separate input and output subspaces enables lrDMD based reduced order controller to stabilize the system at the same rank $5$ approximation. For the rank $9$ approximation in Fig. 2c, we see that the DMD, OMD and lrDMD controlled system eigenvalues are very close to each other and also close to the eigenvalues of the optimally controlled system. Thus with sufficient rank approximation all three methods can stabilize the system but lrDMD can stabilize the system at a rank lower than OMD or DMD.

5.1.3 Optimal actuator location

Optimal actuator placement is a challenging problem in stability theory especially when applied to large scale fluid systems. Reduced order models make the investigation of optimal actuator location computationally feasible. We use the reduced order feedback controllers to find the optimal actuator location $x_{a}$ (in Eq. 13) in the computational domain for a fixed variance in the spatial support of the actuator. The authors in chen2011h find the optimal actuator and sensor location that minimizes the $H_{2}$ norm of the controlled system natarajan2016actuator . We find the optimal actuator location that minimizes the supremum of the cost of the controlled system described in Eq. (12) for all possible initial conditions. We employ brute force sampling of actuator locations while using the reduced-order model approximations.

Let the control input for the GL system (Eq. (11)) be given by the feedback law $u_{k}=-Kq_{k}$. The controlled dynamics and the cost function can be re-written as

[TABLE]

for some given initial condition $q_{0}$. Let $\underline{Q}=(Q+K^{*}SK)$ and $\underline{A}=(A-BK)$. We are interested in the following infinite summation,

[TABLE]

where $F=\sum_{k=1}^{\infty}(\underline{A}^{*})^{k}\underline{Q}\underline{A}^{k}$. Notice that for a given dynamical system, the cost solely depends on the initial condition. When the controlled system governed by $\underline{A}$ is stable, i.e. all the eigenvalues of $\underline{A}$ lie within the unit circle, $F$ can be efficiently computed by solving the following discrete Lyapunov equation

[TABLE]

The maximum value of the cost function is given by the largest eigenvalue of $F$ denoted by $\lambda_{\text{max}}(F)$. The initial condition that yields this cost is the eigenvector of $X$ that corresponds to the largest eigenvalue. We use the reduced-order approximation of the state transition matrix $A$ to construct low-order feedback gain $K$ and compute the values of $\lambda_{\text{max}}(F)$. We then conduct an exhaustive search over all possible actuator location in an interval to find the optimal actuator location. The performance of the reduced-order method is measured by how close this location is to the true optimal actuator location.

Figure 3 shows the supremum of the cost over all possible initial conditions for different actuator locations $x_{a}\in[-7,1]$ using DMD, OMD and lrDMD approximation compared with the true full-order system. The optimal actuator location is given by the minima of this plot. Figure 3a shows the results for a rank 5 approximation. lrDMD outperforms both DMD and OMD in matching the cost over the range of actuator locations considered. The DMD approximation attains its minima at $x_{a}=0$ whereas the true optimal actuator location is $x_{a}=-2$. Optimal actuator locations given by OMD and lrDMD agree with the true optimal actuator location. Figure 3b shows the same results for a rank 9 approximation using DMD, OMD and lrDMD. At this rank all three methods perform equivalently and provide accurate costs for all actuator locations considered.

5.2 Flow past a flat plate

In this section we demonstrate the performance of the controllers on stabilizing the two-dimensional uniform flow approaching an inclined flat plate. The freestream flow is at a low Reynolds number of $100$ and the flat plate is inclined at an angle of $35^{\circ}$. At these conditions, it has been shown ahuja2010feedback that the steady state of the flow is unstable and the flow exhibits periodic vortex shedding. The goal of the reduced-order controllers is to bring the system back to steady state from different initial conditions in the transition process.

The flow is simulated using the fast immersed boundary method developed in colonius2008fast . It is an efficient method to solve incompressible Navier-Stokes equations based on an immersed boundary formulation. In order to achieve uniform flow conditions in the far field, a multi-domain approach is employed. The domain of interest is considered to be embedded in a series of domains, each twice-as-large as the preceding but with the same number of the uniform grid points. The numerical parameters of the flow are taken to follow the work of ahuja2010feedback . The grid size used is $250\times 250$ and the domain of interest is given by $[-2,3]\times[-2.5,2.5]$ where the lengths are non-dimensionalized by the chord length of the flat plate, $L$. The center of the flat plate is located at the origin. Five domains, each with the same number of grid points are used for an effective computational domain that is $2^{4}$ times larger than the domain of interest. The time-step is taken as $dt=0.01L/U_{\infty}$ where $U_{\infty}$ is the freestream velocity.

5.2.1 Steady state

The computation of the steady state of an unstable flow is generally more difficult than a stable flow configuration. We are interested in the initial condition that is a fixed point of the governing equations of the flow ahuja2009reduction . Using a computational Fortran wrapper around the immersed boundary method code, we perform Newton-GMRES iterations on the nonlinear solver. We find the vorticity field $q$ that satisfies

[TABLE]

where $\phi_{T}(\cdot)$ is the nonlinear solver that advances the solution $q_{k}$ at timestep $k$ to $q_{k+1}=\phi_{T}(q_{k})$ after $T$ timesteps. For the purpose of this study we choose $T=50$ and iterate until a convergence tolerance of $10^{-8}$ in the $L_{2}$ norm. Figure 4a shows the vorticity contours and velocity streamlines of the steady state while Figure 4b shows the coefficient of lift as instabilities grow with the unstable steady state as the initial condition.

5.2.2 Snapshots

The snapshots for the study are generated by the nonlinear impulse response of the flow to an actuator ahuja2010feedback located near the leading edge of the flat plate. The system can be described as

[TABLE]

where $q_{k}$ is the vorticity field at iteration $k$, $\phi_{T}$ is the nonlinear solver that advances the vorticity field by $T$ timesteps, $B$ is the actuator and $u_{k}$ is the control input at iteration $k$. The actuator is a simple model of localized body force ahuja2009reduction at the actuator location near the leading edge of the flat plate. The instantaneous vorticity field generated by impulse control input of the actuator is

[TABLE]

where $r_{i}^{2}=(x-x_{c,i})^{2}+(y-y_{c,i})^{2}$ for $i={1,2}$. The constants $a$ and $c$ determine the shape and strength of the control, respectively. Table 4 shows values of the constants used for this study. Figure 5 shows the vorticity field generated by the actuator relative to the position of the flat plate in the computational domain.

Getting a set of snapshots that captures the essential physics is crucial for the performance of data-driven methods. Therefore, an important step in building data-driven reduced-order controllers is to identify the nature of instabilities in the flow and select data that best captures the behaviour of interest. The goal is to select a set of snapshots that capture the near-linear behaviour of the flow before the nonlinearities dominate the system dynamics. Figure 6 shows the evolution with time of drag coefficient of the impulse response of the flat plate. The two regions marked as exponential growth and limit cycle can be approximated by linear dynamics. We are only interested in the exponential growth region and we choose to sample this region for our study.

5.2.3 Closed loop control

We model the system as linear about the steady state of the flow. To this end, we subtract the computed steady state from the flow snapshots and build a reduced-order model for the perturbations about the steady state. The reduced-order model can be described as

[TABLE]

where $\tilde{q}_{k}=q_{k}-\overline{q}$ and $\overline{q}$ is the steady state vorticity field. This can be considered as a model of the dynamics of perturbations in the flow around the steady state.

The control is performed using a Fortran wrapper around the nonlinear solver. The controller is activated at three different points in the transition process, $t_{0}=\{170,190,210\}$ (refer to Fig. 4b), in separate simulations, to show the effect of the nonlinear dynamics of the flow. The nonlinear effects grow larger as the flow deviates from the steady state with time. At $t_{0}=170$, the flow shows exponential growth in perturbation which is predominantly linear behavior. At $t_{0}=190$, the flow has entered an algebraic growth in perturbation magnitude and $t_{0}=210$ is just before the flow enters a limit cycle behavior. At $t=\{190,210\}$ the perturbations are too large for the linear assumption to be valid.

Equation (15) is used to generate snapshots to build the data matrix. We choose $T=20$ and collect $200$ snapshots to cover multiple period of oscillations of the exponentially growing perturbations. Linear approximations of rank $20$ using DMD, OMD and lrDMD are used to generate the LQR controllers. Since $T=20$, the feedback control is applied after every $20$ timesteps of the nonlinear flow solver. A rank of $20$ is high enough to faithfully reconstruct the dynamics but much lower than the dimension of the full order system which is of $O(10^{4})$. Figure 7 shows the evolution of $C_{D}$ as the controller is switched on at the three time points. Figure 7a shows DMD, OMD and lrDMD are all able to prevent the growth of perturbations in the flow when the controller is activated at $t_{0}=170$. This is expected since the perturbations are small and the flow is in the linear regime. In Fig. 7b at $t_{0}=190$, OMD and lrDMD show very similar control performance and are able to suppress the vortex shedding earlier than DMD. Finally in Fig. 7c, lrDMD clearly outperforms both OMD and DMD in bringing the flow back to the steady state. This shows that in the linear regime all three methods are equivalent in performance. As the perturbations grow and nonlinear effects in the flow become significant, the additional flexibility of lrDMD that allows a different input and output subspaces results in a better control performance compared to OMD and DMD.

Figure 8 shows the actuation of the controller with time for the three activation times of the controllers using DMD, OMD and lrDMD used in Fig. 7. The results follow the same trend as Fig. 7. In Fig. 8a when the controller is activated at $t_{0}=170$, all three controllers show similar performance. At $t_{0}=190$ in Fig. 8b, OMD and lrDMD incur much lower cost compared to DMD and in Fig. 8c the cost with lrDMD controller is clearly much lower than DMD and OMD controllers.

We repeat this numerical experiment using data matrices with consecutive snapshots that are $50$ and $100$ iterations apart. Table 5 shows the results of feedback based control at different initial conditions of the flow for DMD, OMD and lrDMD based controllers. As the separation between the flow snapshots increases, the actuation is being applied further apart in time which makes it harder to suppress the vortex shedding process. OMD and lrDMD controllers are successful in controlling the flow at all conditions except when the snapshots are $100$ timesteps apart and the initial condition is at $t_{0}=210$ at which the flow becomes unstable. DMD however fails more often especially when the snapshots are $100$ timesteps apart, it is only able to control when $t_{0}=170$.

6 Conclusion

To summarize, we introduce a method to approximate the dynamics of an unsteady fluid flow by a rank-constrained linear representation. We solve for the optimal linear map of a fixed user defined rank that reconstructs the data snapshots by minimizing the $L_{2}$ norm of the reconstruction error. Two methods of solving the optimization problem are presented and their limitations are discussed. The first method is a subspace projection method that is very fast and provides a good approximation to the optimal solution of the optimization problem. The second method employs a gradient descent approach which is guaranteed to converge at a local minimum. Reduced-order models generated using lrDMD are shown to produce lower reconstruction errors compared to Optimal Mode Decomposition (OMD) and Dynamic Mode Decomposition (DMD) while incurring comparable computational times. These reduced-order models are used to construct low-rank full-state feedback controllers to control model fluid flows.

We employ LQR based feedback control using the reduced-order models on the linearized Ginzburg-Landau equation in the globally unstable regime. lrDMD is able to stabilize the system at a much lower rank approximation as compared to DMD. The reduced-order models are also used to find optimal actuator location using brute force sampling approach. The true optimal actuator location is obtained using lrDMD at a rank 5 approximation, whereas DMD requires a rank 9 approximation to provide the true optimal actuator location. We also employed LQR based feedback control on unsteady flow over an inclined flat plate. OMD and lrDMD perform equally well in controlling the flow whereas DMD incurs much higher costs in the control for all cases.

This paper shows the potential of using a linear map that is a mapping between different subspaces to construct reduced-order controllers for unsteady high-dimensional systems. We observe that these linear maps fit the data better and can lead to better performing reduced-order controllers.

Acknowledgements

The authors would like to thank Mr. Daniel Floryan and Dr. Clarence Rowley for their help with the code used to simulate flow past an inclined flat plate. This work was sponsored, in part, by the Office of Naval Research (ONR) as part of the Multidisciplinary University Research Initiatives (MURI) Program, under grant number N00014-16-1-2617.

Appendix A Numerical methods

In this section we present the two methods to solve the lrDMD optimization problem.

A.1 Subspace projection method

In this method, we use iterative subspace projection to find a good approximation of the optimal solution. We first make the observation that $Q_{R}$ is an orthogonal projection matrix in the column space of $X^{T}R$. This means that there exists an orthogonal matrix $C_{R}$ such that $C_{R}C_{R}^{T}=Q_{R}$ and $\text{Im}(C_{R})=\text{Im}(Q_{R})$. Substituting $Q_{R}$ with $C_{R}C_{R}^{T}$, we get the following cost function,

[TABLE]

with the constraint that $\text{Im}(C_{R})=\text{Im}(Q_{R})$. If $C_{R}$ is fixed, the optimal solution for $L$ is given by the left singular vectors of $YC_{R}$. Finding $C_{R}$ for a fixed $L$ under the given constraint is not as trivial. If the left singular vectors of $X^{T}R$ span the same space as the left singular vectors of $Y^{T}L$, then optimal $C_{R}$ under the constraint will just be the left singular vectors of $X^{T}R$. As an approximation, we try to find $R$ that minimizes the distance between $X^{T}R$ and $Y^{T}L$ by solving the following ‘Orthogonal Procrustes Problem’ schonemann1966generalized

[TABLE]

The closed form optimal solution to this problem is $R=UV^{T}$ where $U$ and $V$ are the left and the right singular vector matrices respectively from the singular vector decomposition of $XY^{T}L$. $C_{R}$ is given by the orthogonal basis of the column space of $X^{T}R$ denoted by $\Pi(X^{T}R)$. Algorithm 1 describes all the steps for this method.

The solution provided by this algorithm relies heavily on the initial guess. We only need an initial guess for $R$ or, effectively, $C_{R}$. In our study, we choose the $r$ leading left singular vectors of $X$ as the initial guess for $R$. Note that this initial guess is the same as the projection subspace used in DMD. Due to the optimal choice of $L$ for a fixed $C_{R}$, the chosen initial guess ensures that the algorithm provides a solution with a reconstruction error at most as high as the error in DMD reconstruction of the same rank.

A.2 Gradient descent method

In this section, we describe the Riemannian gradient descent method employed to solve the optimization problem. For a thorough review on Riemannian optimization on the Grassmanian manifold the reader is referred to edelman1998geometry . For a function $F(L)$ defined on the Grassmanian manifold where $L\in\mathcal{G}_{r,m}$, the gradient $\text{grad}\,F(L)\in\mathcal{T}_{L}\mathcal{G}_{r,m}$ (tangent space to the manifold at the point $L$) is defined as

[TABLE]

Similarly, the action of the Hessian on any tangent vector $dL\in\mathcal{T}_{L}\mathcal{G}_{r,m}$ is defined as

[TABLE]

Further details of the method and its implementation can be seen in sashittal2018low . We employ the Trust-region algorithm absil2007trust to solve the optimization problem using the MATLAB package ManOpt boumal2014manopt .