Stable Dynamic Mode Decomposition Algorithm for Noisy Pressure-Sensitive-Paint Measurement Data

Abstract

Open AccessTechnical NotesStable Dynamic Mode Decomposition Algorithm for Noisy Pressure-Sensitive-Paint Measurement DataYuya Ohmichi, Yosuke Sugioka and Kazuyuki NakakitaYuya Ohmichi https://orcid.org/0000-0002-0812-1267Japan Aerospace Exploration Agency, Tokyo 182-8522, Japan*Associate Senior Researcher, Aeronautical Technology Directorate, 7-44-1 Jindaiji-Higashi, Chofu; .Search for more papers by this author, Yosuke SugiokaJapan Aerospace Exploration Agency, Tokyo 182-8522, Japan†Researcher, Aeronautical Technology Directorate, 7-44-1 Jindaiji-Higashi, Chofu. Member AIAA.Search for more papers by this author and Kazuyuki NakakitaJapan Aerospace Exploration Agency, Tokyo 182-8522, Japan‡Senior Researcher, Aeronautical Technology Directorate, 7-44-1 Jindaiji-Higashi, Chofu. Senior Member AIAA.Search for more papers by this authorPublished Online:13 Dec 2021https://doi.org/10.2514/1.J061086SectionsRead Now ToolsAdd to favoritesDownload citationTrack citations ShareShare onFacebookTwitterLinked InRedditEmail AboutI. IntroductionAN increasing number of large-scale time-series datasets are being generated with the development of numerical and experimental techniques. For example, pressure-sensitive paint (PSP) [1–3] and particle image velocimetry (PIV) [4–6] have been used for aerospace fluid analysis to obtain the spatial and temporal distribution of the flowfield. To gain deeper insights from multidimensional time-series data and use them for modeling fluid flow, data analysis techniques based on modal decomposition are being actively studied [7–9]. Dynamic mode decomposition (DMD) [10] is one of the most commonly used modal analysis methods along with proper orthogonal decomposition (POD) [11]. DMD has been applied in various studies [2,12–16] due to its advantages in extracting both spatial modes and their associated temporal behavior.In this study, we investigate the DMD method for noisy data, particularly for unsteady PSP measurement data. PSP is an optical measurement method for the pressure field based on oxygen quenching of luminescence. It has the advantage of enabling pressure measurements with high temporal and spatial resolution by using a fast-response porous binder. However, when the pressure fluctuation of the observation target is small, such as in a flowfield with low dynamic pressure, it is difficult to obtain a sufficiently large signal [17]. Therefore, when applying DMD to PSP data, a DMD algorithm is required that performs accurately and stably even in the presence of noise.In DMD analysis, we define the input data matrix, which is a sequence of m snapshots, as Ψ=[ψ1,ψ2,…,ψm], where ψi(i=1,2,…,m) is a column vector representing the ith snapshot, and subscript i corresponds to the time of each snapshot. In addition, we define two matrices, X=[ψ1,ψ2,…,ψm−1] and Y=[ψ2,ψ3,…,ψm]. In the standard DMD algorithm [10], the DMD mode is defined as the eigenmode of the matrix A=YX+, where X+ is the pseudoinverse of X; that is, A is computed as the least-squares (LS) solution of Y=AX.Various DMD methods using noise-aware [8,18–20], sparse representation [21,22], variable projection [23], and ODE-based [24] approaches have been proposed. In particular, Dawson et al. [18] and Hemati et al. [19] pointed out that the LS regression used in standard DMD implicitly assumes that noise is only included in Y and that there is no noise in X, which causes bias in the DMD eigenvalues if the input data matrix includes noise. They proposed total least-squares (TLS) DMD that computes A using TLS regression, which takes into account the noise in both X and Y. A number of studies [20,25] also reported that TLS DMD outperformed standard DMD in terms of the accuracy of DMD eigenvalues. However, it is well-known that TLS regression is prone to computational instability in ill-conditioned problems. Fierro et al. [26] mathematically derived the relationship between LS and TLS solutions using singular value decomposition (SVD) analysis and demonstrated that the contribution of components with small singular values can be larger in TLS solutions than in LS solutions for data with large noise. In experimental measurement data, such as unsteady PSP data, the noise component may be too large to be ignored (i.e., the signal-to-noise ratio is small) when the signal component of the observation target is small. In this case, the system to be solved by the DMD algorithm is likely to be ill-conditioned. However, the computational stability of the TLS DMD algorithm has not been sufficiently investigated.In this study, we investigate the stability of DMD algorithms to noisy data through DMD analysis of a numerical experiment and practical PSP measurement data. We also apply the truncated TLS (T-TLS) regression and optimal truncation level selection proposed by Fierro et al. [26] to DMD algorithm. We evaluate the effectiveness of the T-TLS DMD algorithm by comparing its results with those of other DMD algorithms.II. Algorithms for Dynamic Mode DecompositionA. Standard and Exact DMDThe typical algorithm for standard DMD is as follows: Take the SVD of X, X=UΣVT.Compute the reduced-order operator, A˜=UTYVΣ−1.Solve the eigenvalue problem of A˜, A˜ϕ˜=λϕ˜.Every nonzero λ is a DMD eigenvalue, and the corresponding DMD eigenvectors ϕ are given by ϕ=Uϕ˜.Tu et al. [27] pointed out that standard DMD eigenvectors lie in the column space of X but should lie in that of Y. This can be achieved by computing ϕ as ϕ=λ−1YVΣ−1ϕ˜. When this equation is used instead of step 4 above, the algorithm is called exact DMD. Note that the eigenvalues of the exact DMD are identical to those of the standard DMD. In this study, the exact DMD algorithm is used.B. T-TLS DMDThe algorithms for TLS DMD and T-TLS DMD are as follows: Perform dimensionality reduction so that r<m/2 is satisfied, where r is the reduced dimension. Letting Pr be the first r POD vectors of Ψ, the dimensionality reduction can be achieved by X˜=PrTX,Y˜=PrTY(1)Pr can be obtained by the SVD of Ψ. If Ψ is too large to apply the batch SVD algorithm, online algorithms, such as incremental POD [22], can be applied.Construct the augmented data matrix Z˜=[X˜TY˜T], which is an m−1 by 2r matrix, and take the SVD of Z˜, Z˜=UΣVT.Partition the matrix V such that V=[V11V12V21V22](2)Here, V11 and V21 are r by k matrices, and V12 and V22 are r by q matrices (with q=2r−k).Compute the reduced-order operator, A˜=V21V11+.Solve the eigenvalue problem of A˜, A˜ϕ˜=λϕ˜.Every nonzero λ is a DMD eigenvalue, and corresponding DMD eigenvectors ϕ are given by ϕ=Prϕ˜.When the regularization parameter k(=1,2,…,r) is k=r, the algorithm is identical to conventional TLS DMD. In the above T-TLS DMD algorithm, the T-TLS regression proposed by Fierro et al. [26] is used instead of TLS regression. By choosing the regularization parameter k in Eq. (2) appropriately, the component corresponding to the small singular values that causes instability is truncated. A similar truncation process has already been proposed by Hemati et al. [19], but they did not introduce the parameter selection algorithm. They used the SVD of the augmented data matrix Z=[XTYT] as a preprocessing (dimensionality reduction), whereas in the proposed algorithm, the SVD of the data matrix Ψ is used as a preprocessing, and the truncation including optimal parameter selection is performed in the low-dimensional space.The choice of the optimal k is nontrivial; in this study, we use k that minimizes the error E(k)=‖Y˜−A˜X˜‖2 as a simple criterion [26]. That is, we calculate E(k) for all k and adopt the k corresponding to the smallest E(k) as the optimal parameter. In the presence of noise, E(k) does not necessarily decrease monotonically as k increases because TLS regression does not necessarily minimize the LS error.III. Results and DiscussionA. Numerical TestFirst, the performance of the proposed and existing methods in the DMD analysis of the Kármán vortex shedding phenomenon behind a square cylinder was investigated. The unsteady flowfield obtained by a two-dimensional numerical fluid simulation was used as input data. The Mach number was 0.2 and the Reynolds number was 100. The computational grid was a Cartesian grid of 511×381 points (in the uniform flow and transverse directions), and the cylinder had its center located at the origin (x/L,y/L)=(0,0). Here, L is the length of the side of the square cylinder. The sixth-order-accurate compact finite difference method [28,29] and eighth-order-accurate filtering [30] were used for spatial discretization, and the third-order three-step TVD Runge–Kutta method [31] was used for time integration. As input datasets to DMD, the velocity field in the wake region x/L=[0,10] and y/L=[−5,5] were sampled into equally spaced 101×101 grid points. The number of snapshots was 400, and the time interval between the snapshots was Δt=0.25L/U∞. To investigate the effect of noise, random normal noise with variance σ2 was added to each snapshot as observation noise. The input data (i.e., velocity) and noise amplitudes were normalized by the uniform flow velocity U∞. The input data were projected onto the r POD vectors using Eq. (1).Figure 1 presents the eigenvalues λ of the first Kármán vortex (St=0.147) and second (St=0.296) modes obtained by the T-TLS, TLS, subspace, and exact DMD algorithms. As mentioned in Sec. II.A, the eigenvalues of the exact DMD are identical to those of the standard DMD. The subspace DMD [25] is one of the state-of-the-art DMD algorithms. The Strouhal number St represents the frequency nondimensionalized by L and U∞. The figure displays the 95% confidence ellipse and the average eigenvalue calculated based on 1000 random trials with observation noise of σ2=0.1. The 95% confidence ellipse indicates the magnitude of the variation of the computed eigenvalues, and the smaller the confidence ellipse, the more stable the algorithm. The eigenvalues of the noiseless case (obtained by standard DMD) are also plotted as true values. The optimal regularization parameters for T-TLS DMD were approximately k=45 and 5 for r=101 and 21, respectively. It was observed that for the case of the reduced dimension r=101, the eigenvalues of exact DMD had a larger deviation from the true eigenvalues than the other methods. The eigenvalues were shifted inward from the circumference of the unit circle, indicating that the growth rate decreased due to noise. In contrast, the average eigenvalues of T-TLS, TLS, and subspace DMD had small deviations from the true values; in particular, the average eigenvalue of TLS DMD was almost identical to the true value. However, the 95% confidence ellipse for TLS and subspace DMD was large, indicating that the variation of the eigenvalues caused by noise was large. In the T-TLS DMD, the variation of the eigenvalues was superior to that of TLS DMD, and it was confirmed that the T-TLS DMD algorithm was stable due to the effect of truncation. The 95% confidence ellipse for exact DMD was also relatively small. The eigenvalue of T-TLS DMD for the second mode (Fig. 1b) indicates that the eigenvalue was slightly shifted to the damping side. This was due to truncation in the T-TLS algorithm, which attenuated part of the signal representing the second mode. Figures 1c and 1d display the results for the reduced dimension r=21. Although the variation of the eigenvalues obtained by each DMD method was reduced by a decrease in the reduced dimension, the T-TLS DMD still had the smallest variation. Figure 2 presents the estimated eigenvalues at several noise levels. The larger the noise, the larger the variation in the eigenvalues, but the T-TLS DMD was the most stable at all noise levels. The eigenvalues of T-TLS and subspace DMD were found to be slightly shifted to the damping side for the case of the largest noise level σ2=0.5.Figure 3 presents a comparison of the DMD eigenvectors obtained by each algorithm. To compare the robustness of each DMD algorithm, the eigenvector with the largest error among the results of random trials with noise is displayed. For ease of comparison, the phases of all modes have been adjusted to match by multiplying by a complex number [32]. For the Kármán vortex mode, all methods exhibited similar distributions; however, the subspace DMD result was somewhat unclear. For the second mode, the results for TLS and subspace DMD had unclear distributions. That is, the coherent structures behind the cylinder were collapsed. In contrast, the eigenvector of T-TLS DMD was clear, indicating that by the regularization effect of truncation, the T-TLS DMD algorithm was stable even when the TLS DMD algorithm was unstable. The eigenvectors of the exact DMD were calculated stably, although there was bias in the eigenvalues. From this analysis, we found that T-TLS DMD calculated DMD eigenvalues and eigenvectors more stably and accurately than other methods even in the presence of observation noise.Fig. 1 Comparison of eigenvalues estimated by T-TLS, TLS, subspace, and exact DMD algorithms with σ2=0.1.Fig. 2 Effects of noise variance on eigenvalues of the Kármán vortex mode with r=101, estimated by T-TLS, TLS, subspace, and exact DMD algorithms.Fig. 3 Eigenvectors for the Kármán vortex mode (upper row) and second harmonic mode (lower row) estimated by T-TLS, TLS, subspace, and exact DMD algorithms with σ2=0.1. The eigenvectors estimated by standard DMD without noise are also shown as a reference.B. Pressure-Sensitive-Paint Measurement DataNext, we tested the effectiveness of the T-TLS DMD on practical PSP measurement data. We analyzed the transonic buffet phenomenon appearing on the NASA Common Research Model (CRM) wing surface [33,34]. The experiment was conducted in a 2 m×2 m transonic wind tunnel (JTWT1) at the JAXA Chofu Aerospace Center. The Mach number was set to 0.85, and the Reynolds number based on the mean aerodynamic chord was set to 2.27×106. The angle of attack was set to 3.71°. The unsteady pressure distribution over the wing surface was measured using polymer/ceramic PSP (PC-PSP) developed by Sugioka et al. [3]. Evaluating the noise level based on the previous study [33], the signal-to-noise ratio was approximately 20 dB for large-amplitude fluctuations at the foot of the shockwave, and O(1) dB for fluctuations downstream of the shockwave. For details of the experiment, please refer to [33].Figure 4 presents the instantaneous distribution of pressure fluctuation CP′ obtained from the experiment. CP′ represents the instantaneous pressure coefficient after subtracting the time-averaged value. Sugioka et al. [33] demonstrated that, in this flow, the shockwave generated on the wing oscillated and a pressure fluctuation pattern, the so-called buffet cell, appeared. The pressure pattern occurring near the spanwise position η=0.6 in Fig. 4 is a buffet cell. In addition, pressure patterns that are not buffet cells can also be seen in this figure. In particular, the mottled distribution over the entire wing and the relatively large fluctuation at the model edges are pseudo-pressure distributions (i.e., noise). The former was caused randomly by camera shot noise and was pronounced in regions with small pressure fluctuations. The latter was due to alignment errors in the PSP process, which were caused by vibrations of the model.The pressure coefficient distribution obtained in this experiment was input to the DMD analysis. The number of snapshots was 500, and the time interval between snapshots was Δt=0.38cMAC/U∞. The reduced dimension was set to r=101. The optimal regularization parameter for T-TLS DMD was k=49. The difference in computation time between DMD algorithms was about 1% of the total computation time. This is because in DMD analysis of high-dimensional data such as spatial distribution data, most of the computational cost is consumed by the SVD-based subspace projection, which is a process common to all DMD algorithms.Fig. 4 Pressure fluctuation (CP′) distribution of the instantaneous flowfield obtained by pressure-sensitive-paint measurement.Figure 5 presents the eigenvalues λ obtained using the T-TLS, TLS, subspace, and exact DMD algorithms. It can be seen that the distribution of eigenvalues differed greatly depending on the DMD algorithm used. Comparing the results of T-TLS and TLS DMD, the eigenvalues of TLS DMD had absolute values greater than 1, indicating that the growing mode appeared, whereas the eigenvalues of T-TLS DMD did not exhibit the growing mode. Most of the eigenvalues of T-TLS DMD were distributed near the unit circle; that is, they were expressed as modes close to steady oscillation. In contrast, the eigenvalues of subspace and exact DMD were distributed inside the unit circle, and many of the modes were damping modes. Figure 6 presents a dependence of T-TLS DMD eigenvalues on the regularization parameter k. When k=101, the algorithm is identical to TLS DMD algorithm. As k decreased, some eigenvalues appeared on the damping side, whereas most of the eigenvalues were still distributed near the unit circle.Figure 7 presents a comparison of the DMD eigenvectors corresponding to St≈0.4 indicated in Fig. 5. According to Sugioka et al. [33] and Ohmichi et al. [13], St≈0.4 is included in the range of the characteristic frequency of buffet cells. Because the eigenvalue distribution differed depending on the algorithm, as illustrated in Fig. 5, the modes close to St=0.4 for each DMD algorithm are displayed. Figure 7a demonstrates that a clear buffet cell pattern was captured by T-TLS DMD (k=49). That is, a pressure fluctuation pattern with periodicity in the spanwise direction appeared near η≈0.6. This pressure fluctuation pattern has also been observed in numerical fluid simulations [13]. The eigenvectors (St≈0.4) of T-TLS DMD with k=25 and 75 were similar to Fig. 7a, whereas a slightly larger noise appeared for k=75 (not shown here). Exact DMD also extracted a clear buffet cell pattern. Exact DMD was thus considered to be a relatively stable algorithm although it had the problem of eigenvalues shifting to the damping side. Although a similar buffet cell pattern appeared in the eigenvector of TLS and subspace DMD, it was more obscure than that of T-TLS DMD. This may be due to the overfitting of the TLS and subspace DMD algorithms to the noise. TLS and subspace DMD had other St≈0.4 modes not shown in Fig. 7, but their distributions also contained noise.Fig. 5 Eigenvalue distributions obtained by T-TLS, TLS, subspace, and exact DMD algorithms.Fig. 6 Eigenvalue distributions obtained by T-TLS DMD with changing the regularization parameter k.Fig. 7 Real parts of eigenvectors estimated by T-TLS, TLS, subspace, and exact DMD algorithms. Modes with the frequency of St≈0.4 are displayed.IV. ConclusionsIn this study, the stability of DMD algorithms to noisy data was investigated. To achieve a stable DMD algorithm, the truncated TLS (T-TLS) regression and optimal truncation level selection were applied to the TLS DMD algorithm. By adding truncation regularization to the TLS DMD algorithm, T-TLS DMD improves the stability of the computation while maintaining the accuracy of TLS DMD. The effectiveness of the T-TLS DMD was evaluated by the analysis of the wake behind a cylinder and PSP data for the buffet cell phenomenon. The results showed the importance of regularization in the DMD algorithm. With respect to the eigenvalues, T-TLS DMD was less affected by noise, and accurate eigenvalues could be obtained stably, whereas the eigenvalues of TLS and subspace DMD varied greatly due to noise. It was also observed that the eigenvalues of the standard and exact DMD had the problem of shifting to the damping side, as reported in previous studies. With respect to eigenvectors, T-TLS and exact DMD captured the characteristic flow patterns clearly even in the presence of noise, whereas TLS and subspace DMD were not able to capture them clearly due to noise.L. UkeileyAssociate EditorAcknowledgmentThis work was supported in part by the Japan Society for the Promotion of Science KAKENHI (Grant No. 20K14958).