DynamisLab/docs/My_Confirmation/Chapters/3.Methodology.tex

% \chapter{Methodology} \label{Chap: methodology}

\section{Numerical Environment}

In this study, the flow dynamics is described by the incompressible NS equations with a force term

\begin{equation}
\begin{gathered}
\nabla \cdot \boldsymbol{u}=0, \\
\rho_f\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \cdot \nabla \boldsymbol{u}\right)=-\nabla p+\mu \nabla^2 \boldsymbol{u}+\boldsymbol{f}_e,
\end{gathered}
\end{equation}

Where $\rho_f$ is the fluid density, $mu$ is the dynamic viscosity, $p$ the pressure, $\bf{u}$ the velocity, and $\bf{f}_e$ the body force acting on the fluid by the structure. $\rho_f$ is the fluid density, $mu$ is the dynamic viscosity, $p$ the pressure, $\bf{u}$ the velocity, and $\bf{f}_e$ the body force acting on the fluid by the structure.

The hydrodynamic simulations are performed using the lattice Boltzmann method (LBM), a mesoscopic approach particularly suited for unsteady flow simulations involving complex boundaries and moving interfaces. The method solves the discrete Boltzmann equation with collision and streaming processes governed by the multi-relaxation-time (MRT) formulation, which enhances numerical stability compared to single-relaxation-time models. The He–Luo velocity model is adopted to ensure strict incompressibility conditions at low Mach numbers, with the discrete evolution equation expressed as:

\begin{equation}
\begin{gathered}
f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\mathbf{M}^{-1} \mathbf{S} \left[ \mathbf{m}(\mathbf{x},t) - \mathbf{m}^{(eq)}(\mathbf{x},t) \right]
\end{gathered}
\end{equation}

where $f_i$ denotes the particle distribution function along the $i$-th discrete velocity direction $\mathbf{c}_i$, $\mathbf{M}$ is the transformation matrix mapping distribution functions to moment space, and $\mathbf{S}$ is the diagonal relaxation matrix. The equilibrium moments $\mathbf{m}^{(eq)}$ are constructed to recover the incompressible Navier-Stokes equations through Chapman-Enskog expansion. A D2Q9 lattice structure is employed for two-dimensional simulations, with the lattice speed $c_s = c/\sqrt{3}$ where $c = \Delta x / \Delta t$ defines the ratio of lattice spacing to time step.

\begin{figure}[t]
     \centering
     \includegraphics[width=0.50\textwidth]{Figures/d2q9.png}
     \caption{Lattice arrangement for D2Q9 scheme.}
     \label{fig:d2q9}
\end{figure}

Figure~\ref{fig:domain} illustrates the computational domain and boundary condition implementation. The primary simulation domain spans $1280 \times 512$ lattice units. Inflow boundary conditions at $x=0$ implement a Dirichlet velocity specification through a modified bounce-back scheme with momentum correction, enforcing $u_x = U_0 = 0.01c$ and $u_y = 0$. Top and bottom walls at $y = \pm 12.8D$ employ analogous bounce-back treatments for no-slip conditions. At the outflow boundary ($x = 64D$), a convective condition is imposed:

\begin{equation}
\begin{gathered}
\frac{\partial \phi}{\partial t} + U_0 \frac{\partial \phi}{\partial x} = 0
\end{gathered}
\end{equation}

where $\phi$ represents velocity components, effectively minimizing vorticity reflections. The cylindrical surfaces utilize curvilinear boundary treatments based on the ghost-node interpolation method, which has demonstrated second-order spatial accuracy in previous studies. This approach reconstructs distribution functions at boundary-intersected links using Lagrangian interpolation, maintaining consistency with the no-slip condition for arbitrarily rotating cylinders. The rotational velocity at cylinder surfaces is imposed through the interpolated boundary velocity in the bounce-back scheme, with the circumferential velocity $\mathbf{u}_{\text{wall}} = a_i (-\sin\theta, \cos\theta)^\top$ for the $i$-th cylinder as defined in last chapter.

% \begin{table}
%   \begin{center}
% \def~{\hphantom{0}}
%   \begin{tabular}{lccc}
%       $D/\delta x$  & $C_D$   &   $C_L$ & $St$ \\[3pt]
%        10   & $xxxx\pm~xx$ & $\pm~xxxx$ & xxxx\\
%        20   & $xxxx\pm~xx$ & $\pm~xxxx$ & xxxx\\
%        40   & $xxxx\pm~xx$ & $\pm~xxxx$ & xxxx\\
%   \end{tabular}
%   \caption{Mesh convergence results}
%   \label{tab:mesh}
%   \end{center}
% \end{table}

% Mesh convergence studies were conducted to determine appropriate spatial resolution, with results summarized in Table~\ref{tab:mesh}. Three resolutions were evaluated at $Re_D = 50$: $D/\delta x = 10$, $20$, and $40$, where $\delta x$ denotes lattice spacing. Key hydrodynamic parameters including mean drag coefficient $\overline{C_D}$, root-mean-square lift coefficient $C_{L,\text{rms}}$, and Strouhal number $St$ were compared against benchmark data from previous investigations. The medium resolution ($D/\delta x = 20$) provided optimal balance between computational efficiency and accuracy, with relative errors below x\% for all monitored parameters compared to high-resolution simulations. This resolution was consequently adopted for the extensive reinforcement learning training phase. For final deterministic validation runs, the high-resolution mesh ($D/\delta x = 40$) was employed to ensure solution fidelity.

\begin{figure}[t]
     \centering
     \includegraphics[width=0.85\textwidth]{Figures/architecture.png}
     \caption{Research environment framework.}
     \label{fig:architecture}
\end{figure}

To address the competing demands of computational accuracy and training efficiency in deep reinforcement learning applications, a specialized GPU-accelerated solver was developed. This implementation employs a hybrid architecture coupling CUDA C kernels for performance-critical computations with Python-based orchestration. During initialization, Python dynamically generates configuration-specific CUDA kernel code based on user-defined parameters, including collision operator selection (MRT in this study) and boundary condition types. The PyCUDA framework compiles these kernels to PTX code at runtime and manages all host-device memory transfers. During simulation execution, control variables are asynchronously transferred to the GPU via PyCUDA streams, allowing concurrent computation and data movement. This design achieves approximately 90\% GPU utilization efficiency by overlapping data transfers with kernel execution. The Python layer maintains full control over simulation state, enabling runtime modification of domain geometry, physical parameters, or even fundamental algorithms without process restart. This capability proves particularly valuable for reinforcement learning workflows, where environmental interactions require dynamic adaptation. On the specified hardware platform (AMD EPYC 7302 CPU with NVIDIA V100 GPU), this implementation completes one pinball flow-through period (approximately xxxx time steps) in under x minutes at training resolution, and approximately x minutes at high resolution for validation cases. While this framework offers significant advantages in flexibility and integration with machine learning ecosystems, certain limitations warrant acknowledgment. Dynamic kernel compilation introduces initial overhead during startup phases, though this becomes negligible in long-running simulations. Frequent host-device communication, while optimized through stream concurrency, may impose bottlenecks for extreme-scale problems exceeding available GPU memory. Furthermore, the dependency on PyCUDA may present compatibility challenges in high-performance computing environments requiring pure MPI implementations. Despite these constraints, the architecture delivers an effective compromise between computational performance and research flexibility, particularly suited for data-driven fluid dynamics investigations.

\section{Deep Reinforcement Learning}

The task of achieving hydrodynamic cloaking and illusion in the non-linear, multi-attractor regime of the fluidic pinball presents a formidable control problem. The system is characterized by a high-dimensional state space, strong non-linearities, and complex, time-delayed responses to actuation. Traditional model-based control methods are ill-suited for this challenge, as deriving an accurate and tractable analytical model of the wake dynamics is practically impossible. Consequently, a data-driven, model-free approach is not only advantageous but essential for discovering the sophisticated control policies required.

\begin{figure}[t]
     \centering
     \includegraphics[width=0.85\textwidth]{Figures/brunton2020machine2.png}
     \caption{Machine learning algorithms may be categorized into supervised,unsupervised,and semisupervised,depending on the extent and type of information available for the learning process. Abbreviations: PCA,principal component analysis; POD,proper orthogonal decomposition. \cite{brunton2020machine}}
     \label{fig:brunton2020machine2}
\end{figure}

The application of machine learning to fluid dynamics is not a recent phenomenon, with early uses in flow-field analysis dating back to the 1990s. In the context of feedback control, several data-driven paradigms have been explored. Pioneering work by Lee et al. \cite{lee1997application} successfully employed Neural Networks (NNs) for turbulence control, demonstrating their theoretical power to approximate arbitrary non-linear control laws. However, this flexibility often comes at the cost of exorbitant computational or experimental resources needed to optimize a vast number of parameters. Other methods, such as Genetic Algorithms (GAs), have also been deployed but typically require the structure of the control law to be prespecified, limiting their ability to discover truly novel strategies in problems with complex, unknown dynamics.

In recent years, Reinforcement Learning (RL), and particularly its deep learning variant, Deep Reinforcement Learning (DRL), has emerged as a transformative approach for such complex control problems. Its success has rapidly expanded from games to challenging physical systems, including robotics, autonomous flight, and the simulation of biological locomotion. Unlike supervised learning, which requires labeled data, or the aforementioned methods that may require some prior system knowledge, RL learns an optimal policy through direct, trial-and-error interaction with its environment. This "end-to-end" learning capability is perfectly suited for discovering control strategies in systems where the underlying physics are too complex to model explicitly. As noted by Brunton et al. \cite{brunton2020machine}, the success of RL in fluid mechanics is critically dependent on a thoughtful formulation of the problem, where the states, actions, and rewards are chosen to reflect the governing physical mechanisms.

In light of these advantages, this dissertation adopts a DRL framework to tackle the challenges of hydrodynamic cloaking and illusion. This approach allows the control system to autonomously discover effective manipulation strategies for the fluidic pinball's rotational actuators, navigating the complex flow physics to achieve the desired global wake structure without any prior knowledge of the governing equations.

The control strategy for hydrodynamic cloaking and wake illusion is formulated through a DRL paradigm, distinguished from conventional model-based approaches by its independence from prior fluid dynamic knowledge. This model-free methodology employs proximal policy optimization (PPO), an actor-critic algorithm particularly effective for continuous action spaces. The learning process evolves through iterative interactions between the DRL agent and the fluidic environment, commencing with randomized exploratory actions and progressively refining control policies via reward-driven optimization. As depicted in Fig.~\ref{fig:framework}, the framework comprises three fundamental components: (i) the environment state $s_t$ captured by velocity sensors (\ref{eq:sensors}) and instantaneous hydrodynamic forces acting on the pinball; (ii) the control action $a_t = [U_F, U_B, U_T]^\top$ dictating rotational velocities of the three cylinders; and (iii) the reward function $r_t$ quantifying control performance.

\begin{figure}[t]
     \centering
     \includegraphics[width=0.75\textwidth]{Figures/brunton2020machine1.png}
     \caption{Deep reinforcement learning scheme \cite{brunton2020machine}}
     \label{fig:brunton2020machine1}
\end{figure}

A critical innovation resides in the neural network architecture, where sinusoidal activation functions ($\sin$, $\cos$) replace conventional ReLU or tanh units. This design is motivated by the spectral characteristics of vortex-dominated flows, as demonstrated in Fig.~\ref{fig:act_func} through periodic signal reconstruction tests. Standard activations exhibit significant phase errors and amplitude decay when approximating vortical signatures, whereas trigonometric functions preserve spectral fidelity with mean squared errors reduced by x\% compared to next-best alternatives. The actor and critic networks thus adopt the form:

\begin{equation}
\begin{gathered}
h_1 = \sin(\mathbf{W}_1 s_t + \mathbf{b}_1), \quad h_2 = \cos(\mathbf{W}_2 h_1 + \mathbf{b}_2), \quad a_t = \mathbf{W}_a h_2 + \mathbf{b}_a
\end{gathered}
\end{equation}

where $\mathbf{W}$ and $\mathbf{b}$ denote trainable weights and biases, enabling explicit representation of phase-dependent flow features essential for vortex synchronization.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.75\textwidth]{Figures/act_func.pdf}
  \caption{Activation function affecting NN performance}
  \label{fig:act_func}
\end{figure}

The reward function constitutes the second pivotal innovation, addressing the fundamental challenge of perceptual manipulation in fluid flows. Rather than relying solely on delayed sensor-based metrics, we integrate direct force measurements with flow similarity assessments to accelerate policy convergence. For hydrodynamic cloaking, the reward combines force minimization and sensor alignment:

\begin{equation}
\begin{gathered}
r_{\text{cloak}} = \exp\left( -w_{C_D} \left| \sum_{i=1}^3 C_{D_i} \right| - w_{C_L} \left| \sum_{i=1}^3 C_{L_i} \right| - w_D \sum_{j=1}^3 \mathcal{D}(v_j, v_j^{\text{ref}}) \right)
\end{gathered}
\end{equation}

where $C_{D_i}$, $C_{L_i}$ are drag and lift coefficients of individual cylinders, $v_j$ denotes velocity at sensor $j$, and $\mathcal{D}$ quantifies the dynamic time warping (DTW) distance between measured and reference signals. The exponential form creates a steep reward gradient near optimal states, simultaneously satisfying physical constraints (vanishing net force implies minimal flow disturbance) and perceptual objectives (sensor alignment ensures invisibility to downstream observers).

For wake illusion, the reward structure maintains analogous components but replaces force minimization with target matching:

\begin{equation}
\begin{gathered}
r_{\text{illusion}} = \exp\left( -w_F \left| \| \mathbf{F}_{\text{pinball}} \| - \| \mathbf{F}_{\text{target}} \| \right| - w_D \sum_{j=1}^3 \mathcal{D}(v_j, v_j^{\text{target}}) \right)
\end{gathered}
\end{equation}

where $\mathbf{F}_{\text{target}}$ corresponds to hydrodynamic forces on the emulated cylinder. 

Training proceeds in discrete episodes, each spanning $400T_0$ ($T_0 = D/U_0$) to encompass multiple vortex shedding cycles. State observations and control updates occur at $0.8T_0$ intervals, balancing responsiveness with computational tractability. The PPO algorithm updates policy parameters $\theta$ by maximizing the clipped objective:

\begin{equation}
\begin{gathered}
\mathcal{L}(\theta) = \mathbb{E}_t \left[ \min\left( \rho_t(\theta) \hat{A}_t, \text{clip}(\rho_t(\theta), 1-\epsilon, 1+\epsilon) \hat{A}_t \right) \right]
\end{gathered}
\end{equation}

where $\rho_t$ is the probability ratio between new and old policies, $\hat{A}_t$ denotes advantage estimates from the critic network, and $\epsilon=0.2$ prevents destabilizing updates. Upon convergence (typically after $10^3$ episodes), deterministic control enforces the mean action $a_t = \mathbb{E}[ \pi_\theta(s_t) ]$ without exploratory noise, ensuring reproducible performance in validation studies. This framework effectively bridges short-term physical constraints with long-term perceptual objectives, enabling the discovery of non-intuitive control strategies for complex flow manipulation.

\begin{table}
  \begin{center}
\def~{\hphantom{0}}
  \begin{tabular}{lc}
        Parameter & Value/Method    \\[3pt]
        Network architecture (actor and critic) & $64 \times 64$ \\
        Activation function & Sin \\
        Actuations per $T_0$ & 1.25 \\
        Length of each episode & $600 T_0$ \\
        Optimizer & Adam \\
        Learning rate (actor) & $3 \times 10^{-4}$ \\
        Learning rate (critic) & $4 \times 10^{-4}$ \\
  \end{tabular}
  \caption{Hyper-parameters used in the DRL}
  \label{tab:drl}
  \end{center}
\end{table}

\section{Dynamic Time Warping}

The analysis of time-ordered data is a cornerstone of modern science and engineering, with applications spanning a vast array of disciplines. A fundamental challenge in this domain is quantifying the similarity between two temporal sequences, especially when these sequences are subject to temporal distortions such as stretching, compression, or phase shifts. While simple point-wise distance metrics like the Euclidean distance are effective for sequences of identical length and perfect alignment, they often fail in real-world scenarios where processes may unfold at varying rates. For instance, the same word may be spoken at different speeds, a person's signature may be written faster or slower, or a biological process may progress at different rates in different individuals. To address this challenge, a more robust and flexible similarity measure is required.

One of the most powerful and widely adopted solutions to this problem is the Dynamic Time Warping (DTW) algorithm. First introduced in the context of speech recognition by Sakoe and Chiba \cite{sakoe1978dynamic}, DTW is a dynamic programming-based method that finds the optimal non-linear alignment between two time series. Instead of comparing the i-th point of one sequence to the i-th point of another, DTW calculates an optimal "warping path" that maps the time axis of one sequence to the other, minimizing the cumulative distance between their corresponding, aligned points. This ability to non-linearly warp time allows DTW to find the true similarity between two sequences, irrespective of temporal misalignments.

The versatility and effectiveness of DTW have led to its adoption far beyond its origins in speech processing. It has become an essential tool in fields requiring robust pattern matching in time-series data. Notable applications include the recognition of online handwriting and signatures, where individual writing speeds can vary significantly \cite{tappert1990state}, and in gesture recognition from video feeds, where movements may not be perfectly synchronized \cite{corradini2001dynamic}. In bioinformatics, DTW has proven superior to simple clustering for aligning gene expression time series, enabling researchers to map corresponding biological states even when they unfold at different rates \cite{aach2001aligning}. Furthermore, the algorithm is a key component in modern information retrieval systems, such as query-by-humming song recognition \cite{zhu2003warping}, time-series database searching \cite{zhu2003warping}, and time-series clustering \cite{niennattrakul2007clustering}, solidifying its status as a fundamental algorithm for analyzing dynamic, time-ordered data.

The objective of DTW is to compare two sequences, \(X=\left(x_{1}, x_{2}, \ldots, x_{N}\right)\) and \(N \in \mathbb{N}\). In our context, $X$ represents the time series of velocity vectors from the target flow (e.g., the undisturbed vortex street or the wake of a target cylinder), and $Y$ is the corresponding time series from our DRL-controlled flow. The elements $x_n$ and $y_m$ are the velocity vectors $(u_x, u_y$ recorded by the sensors at discrete time steps.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.85\textwidth]{Figures/dtw1.png}
  \caption{(a) Illustration of a warping path and (b) its interpretation for some sequence $X$ of length $N = 9$ and some sequence $Y$ of length $M = 7$. Each cell $(n, m)$ belonging to the warping path is indicated by a red dot and corresponds to an alignment between the elements $x_n$ and $y_m$ indicated by a red bidirectional arrow. \cite{muller2015fundamentals}}
  \label{fig:dtw1}
\end{figure}

The first step is to define a local cost measure, $c(x_n,y_m)$, which quantifies the dissimilarity between any two points in the sequences. For our vector-valued velocity data, the Euclidean distance is a natural choice:

\begin{equation}
\begin{gathered}
c(x_n,y_m)=\|x_n-y_m\|
\end{gathered}
\end{equation}

Evaluating this cost for every pair of elements from $X$ and $Y$ yields and $N\times M$ cost matrix, $C$, where $C(n,m)=c(x_n, y_m)$. This matrix represents the landscape of local dissimilarities between the two signals. Regions of low cost (dark areas in visualizations) indicate moments where the two signals are locally similar. The goal of DTW is to find an optimal alignment, or warping path, through this cost matrix. A warping path $P=(p_1,\dots,p_L)$ is a sequence of index pairs $p_\ell=(n_\ell,m_\ell)$ that maps the indices of $X$ to the indices of $Y$. This path must satisfy three critical conditions:

\begin{equation}
\begin{aligned}
&\text{Boundary condition: }p_{1}=(1,1) \text{ and } p_{L}=(N, M).\\
&\text{Monotonicity condition: }n_{1} \leq n_{2} \leq \ldots \leq n_{L} \text{ and } m_{1} \leq m_{2} \leq \ldots \leq m_{L}.\\
&\text{Step size condition: }p_{\ell+1}-p_{\ell} \in\{(1,0),(0,1),(1,1)\} \text{ for } \ell \in[1: L-1].\\
\end{aligned}
\end{equation}

The total cost of a given warping path $P$ is the sum of the local costs of all the cells it passes through:

\begin{equation}
\begin{gathered}
c_{P}(X, Y)=\sum_{\ell=1}^{L} \mathbf{C}\left(n_{\ell}, m_{\ell}\right) .
\end{gathered}
\end{equation}

An optimal warping path, $P^\star$, is the path that minimizes this total cost. The DTW distance, denoted as $\operatorname{DTW}(X,Y)$, is formally defined as the total cost of this optimal path:

\begin{equation}
\begin{gathered}
\operatorname{DTW}(X, Y) =c_{P^{*}}(X, Y)=\min \left\{c_{P}(X, Y) \mid P \text { is an }(N, M) \text { warping path }\right\} .
\end{gathered}
\end{equation}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.85\textwidth]{Figures/dtw2.png}
  \caption{DTW algorithm based on dynamic programming. \cite{muller2015fundamentals}}
  \label{fig:dtw2}
\end{figure}

While the number of possible warping paths is exponential, the optimal path and its cost can be found efficiently in $O(NM)$ time using dynamic programming. This method works by building up an accumulated cost matrix, $\mathbf{D}$, where each element $\mathbf{D}(n,m)$ stores the minimum accumulated cost for aligning the prefixes $X(1:n)$ and $Y(1:m)$. The matrix $\mathbf{D}$ is computed recursively:

\begin{equation}
\begin{aligned}
\mathbf{D}(n, m)=\mathbf{C}(n, m)+\min \left\{\begin{array}{l}
\mathbf{D}(n-1, m-1) \\
\mathbf{D}(n-1, m) \\
\mathbf{D}(n, m-1)
\end{array}\right.
\end{aligned}
\end{equation}

The final DTW distance is the value in the top-right corner of this matrix, $\operatorname{DTW}(X,Y)=\mathbf{D}(N,M)$. Once the accumulated cost matrix is computed, the optimal warping path $P^\star$ can be found by backtracking from $(N,M)$ to $(1,1)$, at each step choosing the neighboring cell that led to the minimum accumulated cost.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=0.85\textwidth]{Figures/dtw3.pdf}
  \caption{The DTW algorithm exhibits broad sensitivity to various types of signal distortion.}
  \label{fig:dtw3}
\end{figure}