Analysis on numerical stability and convergence of RANS turbulence models from the perspective of coupling modes

Reynolds-averaged Navier-Stokes simulations are still the main method to study complex flows in engineering. However, traditional turbulence models cannot accurately predict flow fields with separations. In such situation, machine learning methods provide an effective way to build new data-driven turbulence closure models. Nevertheless, a bottleneck that the data-driven turbulence models encounter is how to ensure the stability and convergence of the RANS equations in posterior iterations. This paper studies the effects of different coupling modes on the convergence and stability between the RANS equations and turbulence models. Numerical results demonstrate that the frozen coupling mode, commonly used in machine learning turbulence models, may lead to divergence and instability in posterior iterations; while the mutual coupling mode can maintain good convergence and stability in the process of iterations. This research can provide a new perspective to the coupling mode for machine learning turbulence models with RANS equations in posterior iterations.

Introduction
Turbulence is considered as the last major unsolved problem in classical physics [1] . Numerical simulation methods play an increasingly important role in the research of turbulent problems. However, direct numerical simulation (DNS) and large eddy simulation (LES) are expensive in computational storage and cost for high Reynolds number flows. Therefore, Reynolds-Averaged Navier-Stokes (RANS) simulation is still the dominant tool for industrial problems in the near future.
Traditional RANS models, especially based on the eddy viscosity models (SA model [2] , k-ε/k-ω model [3]~ [5] , Menter's SST model [6] ), have been widely used to solve turbulence problems in engineering. These models can generally predict reliable aerodynamic coefficients for attached flows around complex configurations. However, the key difficulty that traditional RANS models encounter is that they cannot accurately simulate turbulent flows with separations.
Consequently, in recent years, researchers have developed many data-driven turbulence closures, including tensor based neural network (TBNN) models [8]~ [9] , machine learning augmented SA models [10]~ [13] , Reynolds stress discrepancy models [14]~ [16] and physics informed neural network (PINN) models [17]~ [18] . The models are trained by machine learning techniques based on DNS or LES data, which can significantly improve the predicted accuracy for separated flows compared with traditional RANS models. Although the recently developed models have made remarkable achievements in turbulence simulation, a serious problem is how to ensure the numerical stability and convergence when a machine learning turbulence model is coupled with the RANS equations in posterior iterations. The TBNN model, constructed by Ling et al. [8] , was able to reach convergence only based on the converged solution of traditional RANS models, and the authors only pointed out that future studies would investigate the feasibility of iterative convergence which the TBNN model was coupled at every iteration. Thompson et al. [19] found that in the channel flows directly substituting Reynolds stresses from DNS databases into RANS equations led to mean velocity fields with large numerical errors, and the errors would increase with the Reynolds numbers. Poroseva et al. [20] also confirmed that unphysical solutions would occur when the fixed Reynolds stresses were directly inserted to the RANS equations. Wu et al. [21] demonstrated that treating Reynolds stresses as an explicit source term in the RANS equations may result in the ill-conditioned problem, and they also considered that the implicit treatment of the Reynolds stresses (decompose Reynolds stresses into the eddy viscosity part and the nonlinear part) could effectively alleviate the convergence problem. Cruz et al. [22] proposed to replace the Reynolds stress tensor with its divergence, the Reynolds force vector, as the machine learning target, which can more accurately reconstruct the average flow fields from DNS data. More recently, Brener et al. [23] elaborately explored different methods to solve the ill-conditioned problem of the RANS equations according to the plane channel, the square duct and the periodic hill geometries, and results indicated that the RANS equations were still ill conditioned in the whole range of cases even when the Reynolds stress tensor was treated implicitly. They also recommended that the Reynolds force vector with implicit treatment using optimized eddy viscosity could mitigate the error propagation to the mean velocity field in RANS equations.
Guo et al. [24] considered that both the model error and the propagation error existed in RANS simulations, and the authors also mathematically derived the propagation error and claimed that it was an inevitable factor affecting the convergence of RANS simulations.
Mean velocity fields will produce large numerical errors when the Reynolds stresses extracted from DNS data with small errors are fixedly injected to the RANS equations, which is a key issue that the machine learning turbulence models are encountering. Many researchers proposed different strategies to circumvent the problem. Most existing work adopted the idea that inserting the fixed Reynolds stresses directly into the RANS equations, which did not consider the coupling effect between the Reynolds stresses and the RANS equations. The commonly used approach is to make the RANS solver firstly converge to the steady state with traditional turbulence models, and re-converge the simulation using the high-fidelity Reynolds stresses. However, this cannot ensure the iterative convergence if the simulation is started from uniform flow-fields. From a new perspective of the coupling mode between turbulence models and the RANS equations, this paper compares the convergence properties with different types of coupling modes, and we find out that the frozen coupling mode may result in numerical divergence and instability for the RANS simulation at high Reynolds numbers.
The remainder of this paper is organized as follows. Section 2 introduces the CFD governing equations and the methodology of RANS coupling modes. Section 3 presents the numerical examples to verify the proposed method, and finally the conclusions are drawn in Section 4.

Methodology
The steady-state RANS equations for incompressible flows can be written as: where u and p are the mean flow velocity and pressure;  is molecular viscosity; and ij R is the Reynolds stress tensor, which needs to be modeled. The term ij R can be decomposed into a linear part and a nonlinear part based on the eddy viscosity hypothesis, which is defined as: where t  denotes the turbulent eddy viscosity; is the velocity strain rate tensor; and ij R  is the nonlinear part of the Reynolds stresses. The momentum equation can be rewritten as: The main purpose to construct turbulence models is to resolve the RANS equations by modeling the unresolved term t  or ij R , which are injected into RANS equations in order to solve the mean velocity field. There are mainly two parts to We discover that the mutual coupling mode can ensure the numerical convergence and the stability in the process of iteration from uniform flows, while the frozen coupling mode may induce divergence and instability in the simulation of separated flows at high Reynolds numbers. This paper mainly investigates the influence of two coupling modes on numerical convergence and stability when the RANS solver is initiated from uniform flows.
According to the Boussinesq hypothesis of linear eddy viscosity, the nonlinear part of the Reynolds stresses is ignored, and Equation (3) can be simplified as first-order closure models: where M  denotes the modeling error introduced by the approximation of turbulence models. In the process of solving RANS equations by numerical methods, the propagation error P  will also occur, which can be expressed as follows: 6 The modeling error M  has no influence on the convergence and the stability since it is a statistic error if a turbulence model is determined. And the propagation error P  is produced during solving the RANS equations. P  approaches to zero if the solver converges; otherwise, the error will continue to increase with the iteration of the solver and eventually lead to numerical instability or divergence.
This paper mainly studies the influence on numerical stability and convergence of the propagation error P  with different coupling modes, and the model error M  is not taken into consideration currently. For the mutual coupling mode, the propagation error is denoted by MC P  , and Equation (5) can be rewritten as: For the frozen coupling mode, we can firstly solve the RANS equations and obtain the turbulent eddy viscosity in the convergent state. And then, re-compute the RANS solver by fixedly inserting the converged turbulent eddy viscosity. The propagation error is denoted by FC P  , that is * ( ) The RANS calculations are conducted with our in-house code which is based on the cell-centered finite volume method. In this paper, Roe scheme [29] is used to evaluate convective fluxes, and interface values are reconstructed by the second-order least-squares approach. The traditional SA model is used to solve the turbulent eddy viscosity. For the semi-discretized form of the governing equations, the implicit Symmetric Gauss-Seidel scheme [30] is adopted to march in time integral, and the local time stepping and residual smoothing techniques are also employed to accelerate convergence. Our in-house RANS solver has been validated completely by simulating many complex configurations in engineering at high Reynolds numbers, and the details can refer to our previous work [31]~ [33] .

Curved backward facing step
The first test case is the curved backward facing step, and the definition of the geometry is the same as Ref. [34]. The step height is H and the Mach number is set as 0.2, and the Reynolds number is 13,700 based on H and inflow velocity. The computational mesh near the step is displayed in Fig. 2, which is locally refined. The total number of elements is 24,849, and the mesh growth rate is set as 1.

Flow around a S809 airfoil at high Reynolds number
The second test case is the flow around a S809 airfoil at high Reynolds number.
The computational mesh is shown in Fig. 5 We further investigate whether the frozen coupling mode can ensure the numerical stability when the RANS solver is based on the converged flow. Firstly, the mutual coupling mode is adopted to iterate 30,000 steps and the turbulent eddy viscosity is obtained in the convergent state. And then, the RANS solver is re-computed by the frozen coupling mode to a large number of iterative steps.
Comparisons of convergence histories and lift coefficients for different coupling modes are shown in Fig. 11 and Fig. 12. With the increase of iterative steps, the mutual coupling mode can consistently maintain numerical stability and convergence.
However, for the frozen coupling mode, the propagation error of the RANS equations will be gradually enlarged with the process of iterations, and both the residual and lift coefficients will eventually present numerical oscillations and non-convergence. Fig.   13 depicts the comparisons of streamlines near the trailing edge at the end of iterative step, which clearly shows that the vortex in the separated region is unsteady and constantly shaking. The results demonstrate that the frozen coupling mode may induce numerical instability even though the RANS solver is based on converged solutions.
Step As for the frozen coupling mode, the turbulent eddy viscosity is fixed, and there is no information feedback during the whole process of numerical calculation, which cannot eliminate propagation errors effectively. As the iterative step continues to increase, propagation errors may gradually accumulate and are enlarged, which will eventually lead to numerical non-convergence and instability. On the contrary, the mutual coupling mode can dynamically adjust the balance between the unclosed term

Conclusions
Reynolds stresses from DNS data are injected to the RANS solver as closure terms, which can cause significant errors for the mean velocity flow-fields. This paper analyzes this problem from a new perspective of the coupling modes. Based on the RANS equations with the traditional SA model, we compare the influences of the frozen coupling mode and the mutual coupling mode on numerical convergence and stability. As a result, the widely used frozen coupling mode may induce divergence and instability for separated flows at high Reynolds numbers even based on converged mean flow-fields; while the mutual coupling mode can dynamically adjust the balance between unclosed terms and mean flows, which can effectively eliminate propagation errors during iterative steps and ensure numerical convergence and stability. The curved back facing step and the high Reynolds flow around the airfoil are used to verify the proposed idea. In conclusion, this paper demonstrates the problem from a novel perspective that the RANS equations cannot obtain accurate mean flow-fields by injecting high fidelity Reynolds stresses. The proposed viewpoint can provide a new breakthrough to solve convergence and stability problems currently encountered in the machine learning modeling.


Learnable Dynamic Temporal Pooling for Time Series Classification

With the increase of available time series data, predicting their class labels has been one of the most important challenges in a wide range of disciplines. Recent studies on time series classification show that convolutional neural networks (CNN) achieved the state-of-the-art performance as a single classifier. In this work, pointing out that the global pooling layer that is usually adopted by existing CNN classifiers discards the temporal information of high-level features, we present a dynamic temporal pooling (DTP) technique that reduces the temporal size of hidden representations by aggregating the features at the segment-level. For the partition of a whole series into multiple segments, we utilize dynamic time warping (DTW) to align each time point in a temporal order with the prototypical features of the segments, which can be optimized simultaneously with the network parameters of CNN classifiers. The DTP layer combined with a fully-connected layer helps to extract further discriminative features considering their temporal position within an input time series. Extensive experiments on both univariate and multivariate time series datasets show that our proposed pooling significantly improves the classification performance.

Introduction
In the last two decades, time series classification has been dominated by nearest neighbor classifiers which utilize handcrafted feature-based representations (Baydogan, Runger, and Tuv 2013;Schäfer 2015) or various distance measures between time series (Zhao and Itti 2018;Yuan et al. 2019b). Recently, there have been several attempts to exploit deep neural networks (DNN) for time series classifiers (Zheng et al. 2016;Zhao et al. 2017;Wang, Yan, and Oates 2017); they do not require heavy crafting on feature engineering or data preprocessing, and easily be applied to multivariate time series as well. In practice, fully convolutional networks (FCN) and residual networks (ResNet) designed for time series classification (Wang, Yan, and Oates 2017) showed the state-of-the-art accuracy among various DNN competitors (Fawaz et al. 2019).
However, the existing convolutional neural networks (CNN) are not able to fully utilize the temporal information of high-level features for their classification. On top of convolutional layers, the CNN classifiers adopt global average pooling (GAP) or global max pooling (GMP) that simply aggregates all hidden vectors along the time axis. By doing so, they can obtain a global representation for an input time series, as well as avoid the overfitting problem with the help of much fewer model parameters. Nevertheless, such global aggregation discards the temporal position of the hidden features, which makes the CNN learn only positioninvariant temporal features. Since each temporal position itself could be a useful feature for discrimination among different classes in time series classification, the global pooling layer eventually degrades the performance of the classifiers.
To address this limitation, we propose a novel pooling method that effectively reduces the temporal size (i.e., length) of network outputs while minimizing the loss of temporal information. Motivated by the observations that time series instances consist of multiple segments with distinct patterns, our dynamic temporal pooling (DTP) outputs a pooled vector for each segment rather than the one for a whole series. The DTP layer produces segment-level representations by aggregating hidden vectors in each segment, thus it enables to model the classification score based on segment-specific class weights. In other words, our CNN classifiers replace the global pooling with the segmentlevel pooling (being followed by a fully-connected layer), which allows extracting further class-discriminative features and improves the classification accuracy. We additionally present the class activation map (CAM) specifically for our DTP layer, indicating how much each temporal region contributes to predicting the class label of an input time series.
The challenge here is to find out consistent segments from input time series instances that are not temporally aligned with each other. To this end, the DTP layer performs semantic segmentation by using dynamic time warping (DTW). We first introduce trainable latent vectors as many as the number of segments to be identified, termed as prototypical hidden series, for encoding the prototypical features of each segment into them in a temporal order. Then, the DTP layer aligns the network outputs (i.e., the series of hidden vectors) with the prototypical hidden series while keeping their temporal order based on DTW; this generates the set of consecutive time points matching with each segment. In the end, we simultaneously optimize the network parameters of CNN classifiers and the prototypical hidden series, thereby both of them collaboratively improve with each other. That is, training the CNN assists to capture the prototypical features of the segments, and also learning the prototypical hidden series helps the CNN to extract discriminative features. Our empirical evaluation on extensive univariate and multivariate time series datasets demonstrates that the proposed pooling technique significantly improves the discrimination power of the CNN classifiers, regardless of its network architecture or pooling operation. The CNN classifiers with the DTP layer also beat nearest neighbor classifiers that use various distance measures by a large margin. Furthermore, we qualitatively show that the DTP layer is capable of providing an interpretable analysis on its classification result by localizing discriminative regions within a target time series.

Related Work Deep Learning for Time Series Classification
With the great success of deep learning, a variety of deep neural networks (DNN) have been applied to time series classification, and they achieved higher accuracy than conventional nearest neighbor classifiers (Fawaz et al. 2019). Among them, convolutional neural networks (CNN) have gained much attention, such as fully convolutional networks (FCN) and residual networks (ResNet) (Wang, Yan, and Oates 2017), because of their capability of capturing local patterns as well as efficient inference by parallel computations. As presented in Figure 1, the stack of multiple convolutional layers outputs the hidden vector at each time point, 1 which eventually encodes high-level features about the local temporal context 2 surrounding the time point. Based on global average pooling (GAP) or global max pooling (GMP), all the hidden vectors are summarized along the time axis into a single vector (Figure 1a), and it is finally used for computing classification scores.
However, such global pooling causes the loss of information about temporal dynamics of the high-level features, and this leads to the limited performance. In case of time series classification, local temporal patterns can have different meanings depending on their temporal positions where they occur (i.e., position-variant), unlike image classification where the position of visual semantic features does not much affect its class label (i.e., position-invariant). Nevertheless, using the fully-connected layer (with local pooling after each convolution) instead of the global pooling layer, e.g., Time-LeNet (Le Guennec, Malinowski, and Tavenard 2016), not only requires a large number of parameters increasing with the length of time series but also makes the network overfitted to the training data (Fawaz et al. 2019).

Differentiable Dynamic Time Warping
Dynamic time warping (DTW) is a popular technique for measuring the distance between two time series of different lengths, based on point-to-point matching with the temporal consistency. Given two times series X and Y of length M and N , the (m, n)-th entry of its cost matrix ∆(X, Y ) ∈ R M ×N represents the distance (or alignment cost) between X m and Y n . The DTW distance between X and Y is defined by the minimum inner product of the cost matrix and any binary alignment matrix A, where A ⊂ {0, 1} M ×N is the set of possible binary alignment matrices whose (m, n)-th entry indicates whether X m and Y n are aligned or not. Each alignment matrix corresponds to a warping path that connects the upper-left (1, 1)th entry to the lower-right (L, T )-th entry using ↓, →, moves. That is, DTW searches for the optimal warping path that minimizes the total alignment cost, and the path eventually represents the best temporal alignment between the two series. The DTW distance and its alignment matrix can be efficiently obtained by dynamic programming based on Bellman recursion, which takes a quadratic O(M N ) cost.
Recently, the continuous relaxation of DTW (Cuturi and Blondel 2017), named as soft-DTW, has been proposed in order to calculate the gradient of DTW with respect to its input series. Instead of the discontinuous (i.e., nondifferentiable) hard-min operation taking only the minimum value in Equation (1), soft-DTW utilizes the soft-min operation with a smoothing parameter γ by adopting the concept of global alignment kernels (Cuturi et al. 2007). The soft-DTW distance between X and Y is calculated by (2) Note that the original DTW is the special case of soft-DTW with γ = 0. The differentiability of soft-DTW facilitates the optimization of distance (or similarity) among time series in deep learning frameworks.

Dynamic Temporal Pooling Problem Formulation
Given a training set of N time series instances from C classes, we aim to train a CNN classifier that accurately predicts the class label of an input time series instance. For each

Temporal Pooling based on Segmentation
The purpose of our temporal pooling is to reduce the temporal size T of the hidden representation (i.e., the output of f ), while minimizing the loss of temporal information in the time series. The key idea is to partition a series of hidden vectors into L segments (L T ) then generate pooled representations by summarizing the vectors in each segment. Formally, the temporal pooling layer outputs the series of pooled vectorsH = [h 1 , . . . ,h L ] ∈ R K×L of length L, whose l-th vector is obtained as follows.
where φ is the pooling operation, and T l = {t l−1 +1, . . . , t l } is the set of consecutive time points belonging to the l-th segment (t 0 = 0, t L = T ). Three functions can be used for the pooling operation: calculating the average value (denoted by avg), the summation value (denoted by sum), and the maximum value (denoted by max) for each latent dimension. A straightforward strategy for segmentation is to partition a whole time series into shorter ones of the same length in a static manner (Figure 1b), but it has several limitations that have to be addressed. First, time points from different time series instances are not temporally aligned in general, which makes it difficult to find absolute temporal locations for segmentation (Figure 2, upper). Furthermore, considering a segmentation task aims to discover distinct temporal patterns which are internally homogeneous, the length of each optimal time series segment cannot be the same as the others' in most cases (Figure 2, lower).
To tackle these challenges, we perform semantic segmentation by matching each time point with its semanticallyclosest segment in a temporal order. To this end, we first introduce a prototypical hidden series P = [p 1 , . . . , p L ] ∈ R K×L of length L which best summarizes the high-level features of L segments. Our proposed method, termed as dynamic temporal pooling (DTP), temporally aligns the series of prototypical hidden vectors (i.e., P) with that of target hidden vectors (i.e., H) by using DTW. Based on the result of DTW alignment, the series of hidden vectors is partitioned into L segments (Figure 3), and each of them is pooled by Equation (3). In this case, the optimal alignment matrix A * between P and H can be obtained by ∆ ∈ R L×T is the alignment cost matrix whose (l, t)-th entry δ(p l , h t ) encodes the distance between p l and h t . We define the cost by using their cosine similarity as follows. 4 δ(p l , ht) = 1 − p l · ht p l 2 ht 2 .
Unlike eligible warping paths of conventional DTW, we need to impose an additional constraint that each time point should be aligned only with a single segment. Thus, we limit A to the set of possible binary alignment matrices representing a path that connects the upper-left (1, 1)-th entry to the lower-right (L, T )-th entry using only →, moves.
As we discussed, soft-DTW is differentiable with respect to its input, thus P can be easily optimized using the gradients ∇ P DTW γ produced by Equation (6). Similar to the original DTW, the soft-DTW distance (and the soft alignment matrix) can be obtained by solving a dynamic program based on Bellman recursion. Algorithm 1 describes the process of forward and backward recursions to compute the alignment cost DTW γ (P, H) and its gradient ∇ P DTW γ (P, H). Please refer to (Cuturi and Blondel 2017) for more details of deriving the algorithms. Note that a single time point should not be aligned with multiple consecutive segments. For this reason, our forward recursion does not allow the ↓ relation in its recurrence (i.e., R l,t depends on only R l−1,t−1 and R l,t−1 ), and accordingly, the backward recursion does not consider the ↑ relation (i.e., E l,t is obtained from E l,t+1 and E l+1,t+1 ).
Learning the parameters of the CNN classifier Instead of the class weight vector w (c) ∈ R K that learns the importance of each latent dimension for class c, we introduce a class weight matrix W (c) = [w L ] ∈ R K×L in order that the class weights are independently modeled for each segment. Using the pooled vectors and the class weight matrices, we calculate the classification score s (c) = L l=1h l · w (c) l , and the posterior probability that an input time series instance belongs to class c is defined as follows: The network parameters W and class weight matrices {W (c) } are optimized by the following classification loss.
log P (y = y n |X n ) (8) To sum up, our CNN classifier adopts the segment-level fully-connected layer for time series classification, by combining Equation (7) with the DTP layer.
Optimization All the parameters including the prototypical hidden series, the network parameters, and the class weight matrices are effectively optimized at the same time, by minimizing the corresponding losses: Note that we consider H = f (X; W) rather than X for DTW alignment (Equations (4) and (6)), so as to perform semantic segmentation in the latent space where the high-level features are embedded. For this reason, P and H collaboratively improve with each other during the training, and each of them respectively captures stereotypical macro-patterns and discriminative micro-ones. Training the network parameters of CNN classifiers (Equation (8)) helps to learn the prototypical high-level features of each segment, while optimizing the prototypical hidden series (Equation (6)) facilitates the CNN to learn further discriminative features by providing more consistent segments.

Class Activation Map for Temporal Pooling Layer
For further explanation on how much each temporal region (or time point) contributes to the classification of a target time series, we tailor a class activation map (CAM) (Zhou et al. 2016;Wang, Yan, and Oates 2017) for CNN classifiers with the DTP layer. CAM has been used to localize the regions relevant to the target class within a time series, but it cannot take into account the segment-specific class weights and also requires the class weights trained for the GAP layer.
Inspired by Grad-CAM (Selvaraju et al. 2017), we generalize the activation map for class c at time point t by The point-wise gradient ∂s (c) /∂h t,k , which implies the position-specific class weights, is used as the weight for each hidden feature h t,k . For example, M (c) t becomes equivalent to w (c) l · h t , if h t is aligned with the l-th segment by the temporal sum pooling. Using the CAM, we can identify discriminative regions or temporal patterns while considering different class weights for each segment.  from the UCR/UEA repository Dau et al. 2018). As the datasets are publicly available as well as collected from a wide range of domains, they have been widely used for the time series classification task.
Baselines As our main baselines, we use CNN classifiers equipped with different temporal pooling layers. For the experiments, we build three types of temporal pooling with avg, sum, and max operations: global temporal pooling (GTP), static temporal pooling (STP) with the fixed pooling size, and the proposed dynamic temporal pooling (DTP). We remark that several DNN architectures that have been studied for time series forecasting, for instance, dilated CNN (Oord et al. 2016), autoregressive model based on recurrent neural networks (Lai et al. 2018;Rangapuram et al. 2018), and transformer (Li et al. 2019), cannot be directly applied to time series classification (i.e., instance-level label prediction). All of them need to be customized or tuned to obtain the representation for an input time series instance; in this sense, our DTP layer can be easily generalized for the architectures as well, but we leave this for future work.
In addition, we compare our classifiers with nearest neighbor classifiers using various distance measures, especially for univariate time series. 5 Since the purpose of our work is to enhance the discrimination power of a single classifier, we exclude ensemble classifiers Lines and Bagnall 2015;Lines, Taylor, and Bagnall 2016) that require highly intensive computations for exploiting dozens of classifiers together. Implementation details We employ two CNN architectures, FCN and ResNet, specifically designed for time series classification (Wang, Yan, and Oates 2017); they showed the best accuracy among various types of DNNs (Fawaz et al. 2019). We implement all the CNN classifiers using PyTorch, and make use of the Numba compiler to compute the forward and backward recursions of the DTW (Algorithm 1) in parallel. 6 In order to eliminate the benefit from hyperparam- eter tuning, we train our classifiers without introducing the weight that balances the two losses, and also fix the number of segments L to 4 for STP and DTP. 7 Table 1 describes the details of hyperparameters used in our experiments.
Evaluation strategy For quantitative evaluation, we conduct the pairwise posthoc analysis (Benavoli, Corani, and Mangili 2016) that statistically ranks different classifiers according to their accuracy over multiple datasets, as done in (Fawaz et al. 2019;Yuan et al. 2019b). We visualize the results by the critical difference (CD) diagram (Demšar 2006) which indicates the average rank of each classifier with thick horizontal lines showing a group of classifiers that are not significantly different (p = 0.05). For all the datasets, we repeatedly train each classifier three times with different random seeds, and report the median accuracy.

Comparison of Different Pooling Layers
We first directly compare the classification accuracy of our proposed CNN classifier (using DTP) with that of the baseline CNN classifier (using GTP). In Figure 4, a single dot represents each dataset, thus how far each dot is located from the y = x line indicates the performance gap between the two pooling methods. For all the cases, we observe that most of the datasets are dotted in the lower right side of each figure, showing that DTP outperforms GTP regardless of its CNN architectures and pooling operations. In particular, the performance improvement of DTP over GTP becomes larger when it is used with the max operation. As the max pooling is effective to detect specific features in general,  the DTP-MAX becomes good at discovering such features segment-specifically, which makes the final representation further class-discriminative. For more statistical evaluation, we also compare different temporal pooling methods (i.e., GTP, STP, and DTP) based on pairwise statistical tests. Figure 5 presents their average rank over a bunch of datasets (85 for the univariate case, and 30 for the multivariate case) with the pairwise statistical differences. Specifically, DTP consistently performs the best among all types of temporal pooling while GTP performs the worst. In case of STP, it shows slightly better but not much statistically different performances than GTP, even though it considers the same number of segments with DTP; this implies that pooling the vectors from same-length segments at fixed temporal positions cannot effectively boost the accuracy of CNN classifiers. On the contrary, our DTP method, which uses variable-length segments identified by DTW, is capable of modeling high-level features depending on each segment, thus its discrimination power improves a lot. We can conclude that DTP successfully utilizes temporal information for its classification compared to GTP and STP.

Comparison with Nearest Neighbor Classifiers
As the nearest neighbor classifier has been one of the powerful benchmarks in the field of univariate time series classification, we compare their performances with ours. For the nearest neighbor classifier, we consider 8 different dis-   (Yuan et al. 2019a), and LSDTW (Yuan et al. 2019b). Figure 6 clearly shows that the CNN classifiers (DTP-MAX) beat all the other classifiers by a large margin. In addition to its higher accuracy, the deep learning approach has some benefits in terms of efficiency and scalability. To predict the class label of an input time series instance, our classifiers require only a single CNN inference with an additional cost of O(LT ) for the DTW alignment between P and H. By contrast, the nearest neighbor classifiers need to perform O(N ) comparisons to find the closest instance from the input time series, and each comparison takes O(T 2 ) for computing the DTW-based distance based on point-to-point matching. This highly limits the scalability of the nearest neighbor classifiers with respect to N and T .

Parameter Analysis on L
We investigate the performance changes of the CNN classifiers (DTP-AVG and DTP-MAX) on the multivariate datasets, increasing the number of segments (i.e., L) from 2 1 to 2 4 . In Figure 7, a single line denotes each dataset, and its color is determined by the optimal value of L that achieves the highest accuracy. The tendency of performance changes is not consistent over most of the datasets, and also the optimal number of segments varies depending on the dataset. As illustrated in Figure 2, time series instances in a dataset share their own temporal patterns and lengths, which are distinguished from the ones in other datasets. For this reason, the optimal number of segments largely depends on such properties of each dataset. This result strongly indicates that finding the optimal L value for a target dataset can further enhance the performance of the CNN classifiers in practice, compared to the case of using the fixed L in our experiments.

Qualitative Analysis
To qualitatively compare the localization performance of GTP and DTP, we visualize their CAM scores by highlighting the input time series proportionally to the scores. In Figure 8a and 8b, the time series instances of class 1 and 2 are colored in red and blue, respectively, and we use the CAM scores after normalizing them in the range of [0, 1] for each time series instance. Although the two CNN classifiers achieve almost the same accuracy for both the datasets (i.e., GunPoint and TwoLeadECG), they highlight different regions as the discriminative temporal patterns that most contribute to predicting its class label. It is worth noting that the localized regions discovered by DTP are more clearly distinguishable between the classes compared to GTP, which results in better interpretability of the CNN classifiers. Furthermore, Figure 8c provides segmentation results obtained from the DTP layer; each time series instance is divided into four segments by vertical lines. The l-th segments (l = 1, . . . , 4) of all instances share similar (or consistent) temporal patterns, even where their original input series are not temporally aligned. The results support that our prototypical hidden series P successfully encodes the prototypical high-level features of the segments in a temporal order, thereby the DTP layer can perform semantic segmentation based on the DTW alignment between P and H.

Conclusion
This paper proposes a dynamic temporal pooling, termed as DTP, which outputs the pooled vectors from temporallyordered segments; this enables CNN classifiers to make use of the segment-level fully-connected layer for time series classification. We present a learning framework to simultaneously optimize the network parameters of a CNN classifier and the prototypical hidden series that encodes the latent semantic of the segments. By finding the optimal alignment between the prototypical hidden series and the hidden representation of a target time series, the DTP layer is able to partition the whole series into a fixed number of segments. Our extensive experiments show that the proposed DTP layer outperforms other baseline pooling layers, and it successfully identifies consistent segments from time series instances that are out of temporal alignment.


extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
Data Freshness and Energy-Efficient UAV Navigation Optimization: A Deep Reinforcement Learning Approach

In this paper, we design a navigation policy for multiple unmanned aerial vehicles (UAVs) where mobile base stations (BSs) are deployed to improve the data freshness and connectivity to the Internet of Things (IoT) devices. First, we formulate an energy-efficient trajectory optimization problem in which the objective is to maximize the energy efficiency by optimizing the UAV-BS trajectory policy. We also incorporate different contextual information such as energy and age of information (AoI) constraints to ensure the data freshness at the ground BS. Second, we propose an agile deep reinforcement learning with experience replay model to solve the formulated problem concerning the contextual constraints for the UAV-BS navigation. Moreover, the proposed approach is well-suited for solving the problem, since the state space of the problem is extremely large and finding the best trajectory policy with useful contextual features is too complex for the UAV-BSs. By applying the proposed trained model, an effective real-time trajectory policy for the UAV-BSs captures the observable network states over time. Finally, the simulation results illustrate the proposed approach is 3.6% and 3.13% more energy efficient than those of the greedy and baseline deep Q Network (DQN) approaches.

I. INTRODUCTION
The rapid deployment of the Fifth-generation (5G) wireless network sets an unparalleled criteria for high-quality wireless connectivity and services [1]. As a result, conventional cellular networks face enormous challenges to meet the stringent requirements for different 5G application types, such as enhanced mobile broadband (eMBB), ultra-reliable and lowlatency communications (URLLC), and massive machine-type communications (mMTC) applications [2] [3]. One of potential solutions is to deploy unmanned aerial vehicles (UAV) in the 5G network environment where UAVs serve the network applications and users as aerial base stations (UAV-BS) [4]. Unlike the conventional wireless network infrastructure, UAV-BSs are more agile and capable of not only providing better coverage but also significantly strengthening the network capability to meet the stringent demands of high capacity with wide coverage, and low delay constraints. However, the deployment and autonomous navigation of multiple UAV-BSs in 5G networks are still challenging due to the limited energy capacity of UAV-BSs. Moreover, in recent years, the concept of edge computing with 5G [5] has also emerged to complement the need for a remote cloud environment for enabling computation oriented communications (COC) applications such as virtual and augmented reality (VR and AR) [6], real-time monitoring and surveillance [7]. As the demand for the upcoming COC applications [8] becomes more prevalent, the need for receiving fresh information data update from different futuristic applications [9] requires a new metric, age of information (AoI) [10], to measure data freshness apart from the traditional performance metrics for the 5G application types.
In case of the UAV-BS navigation, the existing research [11][12][13][14][15][16] mainly focuses on the path planning and placement of UAV-BSs in the network along with communication and energy constraints. However, to incorporate the computation oriented applications at the edge network requires an agile UAV-BS navigation that not only enhances the energy efficiency [17] but also ensures the up-to-date data delivery on time for seamless operation of COC applications at the edge computing platform.
Under the above circumstances, we focus on optimizing the UAV navigation considering energy efficiency and the AoI context in the 5G enabled edge computing environment. The main contributions of the paper are summarized as follows, • First, we formulate an energy-efficient UAV-BS navigation optimization problem in the edge computing network under the contextual constraints such as energy, navigation, and AoI metric, and then we show that the formulated problem is NP-hard. • Second, we employ a deep reinforcement learning technique, deep Q-network with experience replay memory, which can achieve energy-efficient UAV navigation under the contextual constraints. With the proposed trained model, an effective real-time trajectory policy for the UAV-BSs can be obtained that captures the observable network states over time. As a result, we design the state, observation, action space and reward explicitly for the proposed deep Q network (DQN) with experience replay which can effectively solve the trajectory optimization problem under the AoI and energy-efficiency constraints. Also, unlike the traditional deep Q network, the proposed model utilizes the benefit of experience replay memory to obtain the optimal trajectory policy while coordinating multiple UAV-BS locations. • Finally, we perform an extensive experimental analysis to evaluate the performance of the proposed approach. We also conduct an extensive simulation analysis to find the appropriate system parameters to empower the learning model with the proper discount factor and AoI performance metric threshold. The results show that the navigation policy that is obtained by applying the proposed DQN with experience replay achieves significant energy efficiency and data freshness compared to the baseline approaches. The remainder of the paper is organized as follows. In Section II, we present an extensive literature review based on the current research. In Sections III and IV, we present the system model and problem formulation, respectively. Section V explains in detail how we solve the proposed optimization problem with deep Q-learning with experience replay. In Section VI, we present the simulation analysis to validate the performance and efficiency of our proposed approach for UAV-BS navigation. Finally, in Section VII we conclude the discussion.

II. LITERATURE REVIEW A. UAV Navigation
In [18], the authors proposed an online deep reinforcement learning approach to enable UAV navigation in a largescale complex environment. The problem formulation for the UAV navigation is based on the partially observable Markov decision process (POMDP) where the authors proposed the actor-critic framework to solve the problem. In [19], the authors focused on capturing the UAV motion while planning the trajectory for UAV navigation through massive multiple-inputmultiple-output (MIMO) technique. The UAV agent makes the navigation decision on the basis of the received signal strengths which are used to train the proposed DQN. A trajectory planning method for UAVs in urban environments is proposed in [20], where the authors considered the UAV's three-dimensional (3-D) environment map to enable navigation to fuse global navigation satellite (GNSS) signals with ambient cellular signals of opportunity. A novel Deep Reinforcement Learning (DRL) algorithm is proposed in [21] for nonholonomic robots with continuous control in an unknown dynamic environment with moving obstacles. A distributed sense-and-send protocol to coordinate the UAVs for sensing and transmission is proposed in [22]. A reinforcement learning technique is therefore applied to solve some of the UAV's key problems related to trajectory control and wireless resource management.

B. Energy Efficiency
In [23], the authors strive to reduce the total energy consumption of UAVs, including both power propulsion and communication-related energy, while meeting the requirement of multiple ground node (GN) communication throughput. The problem formulation considers the issue of energy minimization by jointly optimizing the allocation of UAV trajectory and contact time between GNs, as well as the overall mission completion time. Using the successive convex approximation (SCA) technique, an effective iterative algorithm is proposed to update the UAV trajectory and contact time allocation at the same time at each iteration, which may converge to a solution that satisfies the Karush-Kuhn-Tucker (KKT) conditions. The goal of the work in [24] is to reduce the propulsion energy consumption of the UAV while meeting the requirement of throughput by optimizing the trajectory of the running track. To transform the problem into a discrete counterpart, a variable discretization approach is used and later, the problem is transformed into a problem of convex optimization where the proposed method can obtain a locally optimal solution. In order to address the critical issue of insufficient on-board UAV energy and CE transmission energy, the authors in [25] focused on maximizing energy efficiency (EE) by jointly optimizing the scheduling of the backscatter devices (BD), the power reflection coefficients of the BDs, the transmission power of the CEs and the trajectory of the UAV. Moreover, the authors considered the BD's throughput and other realistic constraints for the problem formulation. In [26], the authors considered minimizing UAV and user equipment (UE) weighted sum energy consumption subject to task constraints, information-causality constraints, bandwidth allocation constraints, and trajectory constraints of the UAV. The UAV-BSs energy-efficient repositioning trajectories are designed in [27] using the Kuhn-Munkres based algorithm, where an Echo State Network based algorithm enables user equipment (UEs) to predict future trajectories.

C. Age of Information (AoI)
In [28], the authors proposed a dynamic programming approach to study a problem of UAV path planning and data acquisition with the concept of AoI metric. In the proposed approach, the authors jointly considered the selection of data acquisition mode, energy consumption at each node, and age evolution of the information collected by the UAVs. The UAV flight trajectory and status update packet scheduling are jointly configured in [29] to achieve the required weighted sum for the age-of-information (AoI) values of various processes at the UAV, referred to as weighted sum-AoI. A deep reinforcement learning (RL) algorithm is proposed to achieve the optimal policy that minimizes the weighted sum-AoI, called the ageoptimal strategy. The authors suggested a combination of ground sensor nodes (SN) and trajectory planning strategy in [30] to strike a balance between the upload time of the SNs and the flight time of the UAVs using dynamic programming in different scenarios. In [31], the authors sought to minimize UAV's total energy consumption by jointly optimizing the association of the Internet of Things Devices (IoTD), the allocation of computer resources, the UAV hovering time, the wireless power duration and the IoTD service sequence. The authors proposed a UAV trajectory planning model for data collection in [32], where the objective is to minimize expired data packets across the entire sensor network system. The authors also simplified the original problem into a minmax-AoI-optimal path scheme due to complex constraints and proposed a reinforcement learning-based strategy for the solution.
Unlike the existing works, in this paper, we design the trajectory policy of UAV-BSs considering the trade-off

III. SYSTEM MODEL
In Fig. 1, we consider a set of given trajectory points, P = {1, 2, · · · , P }. The trajectory points in P are covered by a set of battery-powered UAV-BS, U = {1, 2, · · · , U }, which act as relay for a set of IoT devices, I = {1, 2, · · · , I}. For simplicity, we consider the set of IoT devices are randomly located at different trajectory points. In this paper, we also consider a single ground station (i.e., ground base station) b which is equipped with the multi-access edge computing (MEC) server (BS-MEC) and acts as the information fusion center that receives information updates from the IoT devices through the UAV-BS relays to support the computation oriented communications applications. Therefore, we assume that the BS-MEC is a point in the network where various communication resources are available to achieve a certain computational accuracy where the timely information updates from the different network sources are essential. Moreover, the BS-MEC b is also considered as a trajectory point within the set P where the set of neighboring trajectory points of b (exclude itself) is denoted by P b = {p ∈ P : (b, p) ∈ } where ⊆ P × P is the set of trajectory links between the trajectory points. The UAV-BS u ∈ U traverses within different trajectory points in P over a finite observation time T . As UAV-BS u ∈ U travels with a trajectory p ∈ P, it gathers the information data packets from the active IoT device i ∈ I located near the trajectory point p ∈ P using the uplink communication channel. Moreover, UAV-BSs u ∈ U uses the backhaul communication link when the BS-MEC b is within the transmission range of UAV-BS u ∈ U. As a result, BS-MEC can receive a fresh information update from different trajectory points and calculate the AoI metric. We assume that the BS-MEC is equipped with an array of mmWave directional antennas and provides a dedicated mmWave spectrum for back-haul communication for the UAV-BS. Moreover, the IoT devices and the UAV-BS are also equipped with directional antennas so that the IoT devices can transmit information updates to the UAV-BS using the non-mmWave spectrum. In this paper, we limit the scope by focusing on the deployment, communication, and navigation of the UAV-BSs u ∈ U for the data relaying from different sources i ∈ I at trajectory points p ∈ P to the ground BS b. Also, for the user association between the IoT devices and the corresponding UAV-BSs at different trajectory points, we apply the default max-signal-tointerference-plus-noise-ratio (SINR) [33] based approach.

A. IoT-to-UAV-BS Communication Model
At the trajectory point p ∈ P, the air-to-ground path loss probability of the UAV-BS u ∈ U with IoT devices i ∈ I is calculated as [34], (1) Here α andα are the environment dependent constants for the LoS and NLoS channels, respectively, where Θ u is the elevation angle of UAV-BS u ∈ U. The intuition of calculating the air-to-ground LoS and NLoS path loss probabilities using (1) is that, in urban/sub-urban environment, the uplink communication link between the the UAV-BSs and the IoT devices may be hindered (i.e., multi-path fading) by the surrounding obstacles (e.g., buildings) unlike IoT devices deployed in the rural environment. In addition, the path loss in decibel (dB) is calculated as [35], ( Here f c is the uplink channel frequency, and and¯ are the attenuation factors for the LoS and NLoS channels, respectively. Using (2), the received signal power from the IoT device i ∈ I at trajectory point p ∈ P to UAV-BS u ∈ U is calculated as, whereP u i,p is the transmit power of IoT device i ∈ I for offloading the data to UAV-BS u ∈ U. At time t, the received SINR for UAV-BS u ∈ U with IoT device i ∈ I at trajectory point p ∈ P is calculated as, Here 10 ) −1 is the received interference of UAV-BS u ∈ U from the other UAV-BSs u ∈ U, u = u which is serving IoT i ∈ I, i = i that is located in different neighboring and overlapping trajectory points p ∈ P, p = p and σ 2 is the noise power. Using (4), the channel capacity at time t is defined as, Here β u is the fixed non-mmWave uplink channel bandwidth that is equally distributed to the IoT devices I at the trajectory point p ∈ P, γ th is the SINR threshold for ensuring successful uplink transmission between IoT devices and UAV-BSs.

B. UAV-BS-to-BS Communication Model
The received power of the ground BS b from UAV-BS u ∈ U is calculated as [36], Here P tx b,u is the transmit power of UAV-BS u ∈ U to BS b, δ b,u is the distance between the UAV-BS u ∈ U and ground BS b, c is the speed of light, f mmW ave c is the carrier frequency of the mmWave back-haul link, G tx u and G rx b are the antenna gains of the transmitter UAV u ∈ U and receiver ground BS b, receptively. At time t, the back-haul capacity of the channel between UAV-BS u ∈ U and ground BS b at time slot t is calculated as, between the UAV-BS and the ground BS b is less than a thresh-hold distanceᾱ, the UAV-BS transmits the information update to the ground base station using β mmW ave b,u . Using (7), the transmission energy of UAV-BS u ∈ U while using back-haul link at time t is calculated as, E mmW ave C. UAV-BS Relay Network Energy Efficiency Metric Design UAV-BS u ∈ U covers the observation area horizontally at a constance altitude h u where different UAV-BSs may maintain different altitudes. The assumption is practical for UAV-BSs according to the Federal Aviation Administration (FAA) regulations for small unmanned aircraft (UAS) operations. Moreover, the UAV-BS trajectory at time t is defined as, . Therefore, the time varying distance covered by UAV-BS u ∈ U horizontally at constant altitude h u is defined as [37], The total mobility energy cost of UAV-BS u ∈ U for covering distance δ u at time t is calculated as, Here is the upper bound of the propulsion power consumption where k 1 and k 2 depends on the UAV-BS design and g = 9.8 m/s 2 is the gravitational acceleration [38].
The total energy efficiency for UAV-BS u ∈ U that covers trajectory points P to serve IoT devices in I over time T is defined as, .

D. Age of Information Model for Ground Station
The AoI metric is used to measure the freshness of information collected by the UAV-BSs from the trajectory points in P where the UAV-BSs act as relay node for the IoT devices. Therefore, at the BS b, the AoI of the trajectory p ∈ P at time t is calculated as, Here ∆ u (p, t) is denoted as the most recent received data packet from the trajectory point p ∈ P at the base station b by the UAV-BS u. At t = 0, we assume ∆ u (p, t) = 0 where we adopt the the just-in-time transmission policy [39]. The average AoI is calculated at the base station b for trajectory points P over time slot T as,

IV. PROBLEM FORMULATION
To formulate the UAV-BS navigation optimization problem under contextual constraints (i.e., trajectory, AoI, energy efficiency constraints), first, we consider each UAV-BS u ∈ U can cover only a sub-set of trajectory points in a given time window T . Moreover, we consider sub-sets for each u ∈ U comprised of trajectory points which are denoted as, P u ⊂ P, P u ∩P u = ∅ where u = u . As a result, the objective of energy efficient UAV-BS navigation optimization problem is to find the cooperative trajectory path configuration of the UAV-BSs that maximizes the total energy efficiency of the UAV-BSs relay network subject to the energy and AoI metric. Therefore, the optimization problem formulation is represented as follows, subject to In the above formulated problem, the constraints (15) (17) and (18) are coupled with the decision variable P u where both the energy efficiency and AoI metric are the functions of P u . Constraint (17) ensures the total energy efficiency of the UAV-BSs where the communication and mobility energy should be greater than a minimum energy efficiency threshold η th . Finally, constraint (18) ensures the average freshness of information updates by configuration P u should be less than an AoI threshold∆ th b . Due to constraint (18) in problem (14), the UAV-BSs jointly navigate different trajectory points under not only the energy efficiency constraint but also the AoI constraint where the performance of the computation oriented communication applications depend on the up-to-date information update from different trajectory points. Therefore, problem (14) is different from the traditional energy efficiency maximization problem for UAV navigation.
The decision problem in (14) can be reduced to a base problem of vertex cover problem (i.e., Maximum Clique Problem) [40] with the corresponding constraints (15)- (18), which is NP-Complete. Similar to the maximum clique problem, problem (14) is combinatorial in nature. Moreover, there is no known polynomial algorithm that can tell, given a solution of (14), whether it is optimal. As a result, we can infer that the decision problem in (14) belongs to the same category of the problem of the vertex cover problem, which is proven to be NP-hard. In the next section, we solve problem (14) with the corresponding constraints (15)-(18) (14), we apply a deep reinforcementlearning model which is the combination of a deep neural network and a reinforcement learning algorithm. Specifically, the proposed DQN approach is comprised of three components: (i) a deep neural network to reduce the dimension of the state space that is used to extract the contextual features (e.g., AoI, energy consumption) necessary for UAV-BS navigation, (ii) An experience replay memory to store the state transitions that the UAV-BS agents observe, and (iii) an reinforcement learning (RL) framework to find the best trajectory policy that achieves the objective of problem (14) with the corresponding constraints (15)- (18). Unlike the state-fo-the-art method such as control methods, the DQN does not need a network dynamic model as it is model free. Moreover, in the proposed approach, the use of experience replay ensures stability by breaking the temporal dependency among the observations used in the training of the deep neural network.
In Section V(A), first, we model the state and action space of problem (14). After that, in Section V(B), we model the reward and control policy based on problem (14) with the corresponding constraints. Finally, in Section V(C), we provide the proposed training and testing model for UAV-BS trajectory policy.

A. State and Action Space
The state space for trajectory policy of the UAV-BS is a four-dimensional state space. At each time step t = {1, 2, · · · , T }, the state or joint observation space of the learning agents (i.e., UAV-BS) is denoted by, S = {s t = (p u current , p end , η, ∆)|η ∈ [0, η th ], ∆ ∈ [1,∆ th b ]}, which corresponds to the current positions of UAV-BS u ∈ U with individual heights h u , target position, energy efficiency of UAV-BSs, and average age for the navigation optimization.
Moreover, the trajectory position for navigation of the UAV-BS u ∈ U is comprised of x u ∈ [0, X u ] and y u ∈ [0, Y u ], where X u and Y u are the maximum coordinate of a particular geographic location. Furthermore, the initial position of each of the UAV-BS is randomly assigned for each trial along with the number of IoT devices. The lower and upper bounds of continuous state variables η and ∆ in the state space are calibrated from the real-world trajectory data.
The action space of the UAV-BSs is the trajectory planning each of the UAV-BS's navigation from one feasible state (i.e., position) to the next state while satisfying the trajectory and communication constraints (i.e., constraints (15)-(18) of problem (14)). The learning agent selects an action a t from the set of available actions upon state s t where a t ∈ A st ⊂ A, and A = {a 1 , · · · , a U } = {P u } u∈U is the configurations of the UAV-BS navigation.

B. Reward and Control Policy
When a learning agent implements action a t , the environment moves to a new state s t+1 and the immediate reward R t+1 with the transition (s t , a t , s t+1 ) is associated and the learning agent receives the reward through feedbacking. In other words, at each state transition, the agent receives the immediate reward which is used to form the trajectory control policy for navigation. For future usage, the control policy is used by the learning agent that maps the current state to optimal control action. The immediate reward is formulated by the instantaneous energy efficiency metric of the UAV-BS's and defined as follows, Here α 1 is a coefficient multiplied to the energy efficiency function and also used to penalize the agent when the constraints are violated.
The objective of the learning agent over T time slots is, therefore, to maximize the future reward which is defined as, Here γ = [0, 1] reflects the trade-off between the importance of immediate and future rewards which, in turn, converges to the optimal control policy. Moreover, the reward function (20) is obtained at time t after learning the current state of the UAV-BSs over the last T time steps duration. Therefore, we define a control policy as π for the agent where the Q-function or the action-value function is defined as, Q π (s, a) =R(s, a) + γ s∈S P s,s V π (s ), Here P s,s is the transition probability of the states in the environment where s = s t+1 , π is the control policy, and action a is enforced though the environment simulator. Here, the state and reward update is based on the information received by the ground BS b. Using (21), we set the goal Algorithm 1: DQN with experience replay for UAV-BS Trajectory Policy Optimization for Navigation 1 Step 1: Initialization 2 Initialize Q(s, a; θ), M, target DQN parameters θ − and construct DQN 3 Step 2: Training DQN with experience replay 4 for e = 1, · · · , E do 5 Initialize S 6 for t = 1, · · · , T do 7 Calculate the energy efficinecy metric of the UAV-BSs using (11) 8 Calculate instant reward R t using (19) 9 Select action a t with given probability . 10 Observe instant reward R t and next state s t

11
Store experience (s t , s t , a t , R t , R t ) in the experience replay memory M

12
Randomly sample minibatch of experiences from M

13
Adopt stochastic gradient descent (SGD) to train the DQN using loss function in (27) 14 Update θ and Q(s, a; θ) 15 Store the Q-network 16 Step 3: Testing UAV-BS trajectory policy for joint UAV-BS navigation 17 Load the stored Q-network of Step 1 18 Retrieve R t of the UAV-BSs at time slot t 19 Retrieve and select joint UAV-BS action a t = max at Q π opt (s t , a; θ) 20 Update trajectory of UAV-BSs based on joint action index and target values of DQN of our model which is to obtain the best control policy π opt . Therefore, the maximum Q-function is defined as, where the discounted cumulative state function is, To derive the optimal control policy π opt , the Q-function is updated as, Here t = t+1 and a = a t+1 where the Q-function is updated using the recursive mechanism and ψ is the learning rate.

C. Training with Experience Replay and Testing UAV-BS Trajectory Policy
The proposed DQN approach learns how to optimally control the trajectory configurations of the UAV-BSs for navigation during the simulation. Therefore, it is vital for the simulation process to train a Q-network where the target value for each trajectory observation environment state is given as, Here we introduce θ k which is the network weight obtained by the training during the k th iteration. Hence, using (25), the loss function of the training network is designed as, Here ρ(s, a) is the probability distribution over the sequences s and actions a, y is the target value of the training network which is derived from (25), and the optimal network weights θ opt are obtained by training. Furthermore, to enhance and stabilize the training of the DQN, we apply the mini-batch method to randomly collect examples from all the training episode steps e t = (s t , a t , R t , s t ) in a fixed size replay memory M t = {e 1 , · · · , E}. As a result, one sample is used multiple times in the training that improves the data efficiency significantly. Therefore, using (25), the loss function in (26) is represented with a uniform distribution over M as, Here U (M) is the uniform distribution over the experience replay memory M and θ − is the stored weight parameters of the target DQN network. In step 1 of Alg. 1, we first initialize the network parameters randomly and introduce the target DQN with the same network structure as the original DQN network (lines 1-2 in Alg. 1). Then, at each training episode steps e t , the energy efficiency metric is calculated using (11) considering the navigation and communication parameters in ((1) -(11)) (line 7 in Alg. 1). As a part of the -greedy policy framework, in the exploration stage, the energy efficiency metric is used to calculate the reward function where the action is derived from the current DQN with the exploration probability . The the reward function is observed considering the AoI context along with other constraints in problem (14) and transit to the next state s t where t = t + 1. Moreover, the exploration stage enables the UAV-BSs to explores all the joint actions for achieving the better reward values that lead toward choosing the appropriate action with the highest energy efficiency. Subsequently, we adopt the mini-batch approach that shuffles the experience from the replay memory buffer at random to remove the correlation in the observation sequence, and thus, smoothing the changes in the observation data distribution (lines 11-12 in Alg. 1). To train the DQN, we adopt the stochastic gradient descent (SGD) algorithm using the training loss function and update the network parameter θ, and network bias (lines 13-14, in Alg. 1). The training process stops when the UAV-BSs arrive at the terminal trajectory and the DQN network is finally stored for testing.
In step 3 of Alg. 1, the stored UAV trajectory policy network in the training phase is used where R t of the UAV-BS at time slot t is retrieved (lines 17-18 in Alg. 1). As a part of the exploitation in the testing phase, the action a t is selected for joint UAV-BS navigation where both the trajectory of the UAV-BSs and the target values of DQN are updated (lines 19-20 in Alg. 1).

VI. PERFORMANCE EVALUATION
In this section, we first address the performance analysis experiment environment through various key metrics. We then describe the outcomes obtained from the experiment and finally provide an in-depth discussion and key observations from the results of the simulation. In order to train the deep neural network (DNN), we consider a neural network architecture with two fully connected (FC) hidden layer with 100 hidden nodes. We also set the experience replay memory size, M = 200. The simulation results are obtained by averaging and normalizing the values over 100 episodes.

A. Experiment Setting
For the performance evaluation of the proposed approach, we consider the simulation settings which is summarized in Table I. In addition, we compare the proposed approach with two baseline approaches which are, • Baseline DQN: The structure of the baseline DQN is different from the proposed DQN with experience replay approach in terms of not having the experience replay memory. • Greedy: In case of the greedy approach, at each timeslot t in an episode, each UAV-BSs u ∈ U co-operatively finds the trajectory paths for navigation that may provide the maximum immediate reward. In addition, the approach applies penalty for violating the system constraints in order to make fair comparison between the proposed and the baseline DQN approaches.

B. Finding Appropriate Discount Factor and AoI Threshold
We present the experimental results for finding the appropriate discount factor (i.e., γ) and AoI threshold (i.e., ∆ th b ) for the proposed DQN with replay memory approach. Here, we first find the discount factor while considering the energy efficiency metric, and then we fix the discount factor to find the appropriate AoI threshold concerning the trade-off between the average AoI and average rewards metrics.
In Table II and III, we evaluate the appropriate discount factor and AoI threshold considering the average energy efficiency metric and we set the number of trajectory points, |P| = 14. In Table II, we observe how the discount factor affects energy efficiency. More specifically, when the discount factor increases in the case of the proposed DQN with replay memory, the energy efficiency metric goes up at 1.73% at the discount factor is γ = 0.5. On the other hand, the energy efficiency metric for the baseline DQN slightly drops by 0.23% at discount factor γ = 0.5 and the trend continues until discount factor γ = 0.6. Similarly, the energy efficiency metric for the proposed approach also slightly fluctuates at discount factor γ = 0.6. However, up to this point, the proposed DQN with replay memory still outperforms the baseline DQN approach. At γ = 0.7, the energy efficiency metric is at the peak for the proposed DQN with replay memory and the performance of the baseline DQN also improves significantly. However, with a higher value of γ, the performance of both the approach decreases drastically. The explanation behind such phenomenon is that with a lower discount factor value (i.e., γ < 0.8), the future reward does not matter much and the agents take action to increase the immediate reward with fewer steps forward. In our case, this is not suitable for reaching toward reasonable energy efficiency. On the other hand, with increasing value of discount factor (i.e., γ > 0.7), the agents give more attention to the benefit that future actions may bring which captures the temporal behavior of the system. When the discount factor is γ > 0.7, the UAV agents care too much for the future reward that eventually leads toward neglecting the immediate reward. So, the performance is degraded and sometimes leads to penalty as the battery of the UAV-BSs drain quickly. Therefore, in our simulation, we set the discount factor γ = 0.7.
In Table III, after fixing the discount factor at γ = 0.7, we investigate the effects of different AoI thresholds to the normalized average reward and average AoI that infer to choose the appropriate AoI threshold for finding the trajectory policy for the UAV-BSs. We investigate the trade-off between the two performance metrics by setting four different levels of AoI threshold which are∆ th b = [0.3, 0.5, 0.7, 0.9]. At ∆ th b = 0.3, we observe that with strict AoI threshold, the normalized average reward is at the minimum whereas the AoI is at the peak point (disregarding the normalized average AoI metric value at∆ th b = 0.9). This is normal because the strict AoI threshold leads to an increasing penalty, and therefore, the performance of the system in terms of average cumulative reward decreases. As we relax the AoI threshold (i.e.,∆ th b > 0.5), the average AoI decreases significantly and the average reward increases. However, at∆ th b = 0.9, the reward value is the maximum but at this point the effect of the AoI threshold is trivial. More specifically, the system at this point the UAV-BSs operate disregarding the freshness of information updates at the ground BS which is not desirable. Therefore, we set the acceptable AoI threshold at∆ th b = 0.7 to balance between the average reward and the AoI metric.

C. Experiment Results
In Fig. 2, we analyze the performance of the proposed DQN with a replay memory with two baseline approaches in terms of average cumulative reward with the increasing number of trajectory way-points. With a small number of trajectory points (i.e., up to |P| = 8), the performance of the approaches is not distinct since the density of the waypoints in the geographic environment is less. However, as the number of trajectory points increases, the performance gaps between the approaches increase significantly. More specifically, the proposed DQN with the replay memory approach outperforms the greedy and baseline DQN correspondingly up to 2.84% and 1.91%. The performance gaps between the approaches further increase with dense trajectory way-points (i.e., |P| = 14) where the proposed DQN with replay buffer outperforms the greedy and the baseline DQN by 5.08% and 4.02%, respectively.
One of the key factors in the proposed model is to consider minimizing the AoI metric while finding a policy for UAV-BSs navigation. Therefore, in Fig. 3, we evaluate the performance of the proposed DQN with replay memory with the greedy and baseline DQN approaches. As we can see from Fig. 3 that the proposed approach outperforms the greedy and baseline DQN by reducing the average AoI correspondingly up to 1.21% and 1.17%. The performance gaps between different approaches increase gradually up to |P| = 8. At |P| = 10, the performance of the approaches slightly fluctuates since in the experiment we independently run the simulations with the different numbers of trajectory points. Therefore, this has some impacts on the experiment results. However, the trend of the proposed approach to reduce the AoI metric than the baseline approaches continues to carry on as the number of trajectory way-points increases. At |P| = 14, the performance gaps among the proposed DQN with replay memory, greedy, and baseline DQN approaches are increased by 0.72% and 1.29%, respectively.
In Fig. 4, we show the performance of the proposed DQN with replay memory with the baseline approaches concerning the average energy efficiency. We observe, the proposed DQN with replay memory gradually performs better than that of the baseline approaches. When the number of trajectory way-points is relatively high (i.e., |P| = 10), the proposed DQN with replay memory is proven to be slightly energy efficient by 0.32% and 1.97% than the greedy and baseline DQN, respectively. However, the proposed approach is significantly energy efficient with the higher number of trajectory way-points (i.e., |P| = 14) where the proposed approach outperforms the greedy and baseline DQN correspondingly up to 3.6% and 3.13%.
The front-haul capacity is limited with an increasing number of IoT devices at the trajectory points and therefore, the bandwidth should be utilized efficiently. In Fig. 5, we compare the efficacy of the proposed approach with the greedy and baseline DQN with the varying number of trajectory points and different IoT device density. The proposed DQN with replay buffer efficiently utilizes the front-haul and back-haul bandwidth while traversing across different trajectory waypoints. The average bandwidth efficiency is quite similar when the number of trajectory points is less dense in the environment and the distance between the points is large. Therefore, the received interference level at the IoT devices which are served   by different UAV-BS is significantly less in all the approaches. However, with a lightly dense trajectory way-point network with the increasing number of IoT devices (i.e., |P| = 10), all the approaches face interference and we observe a slight decrease in bandwidth efficiency. Nevertheless, the proposed DQN with replay buffer still outperforms the greedy and baseline DQN by 2.41% and 2.87%, respectively in terms of ensuring bandwidth efficiency. Fig. 6 illustrates the utilization of the UAV-BSs or network resources under the proposed approach and the other two baseline approaches. As the number of trajectory waypoints increases, the number of IoT devices using the network resources of the UAV-BSs also increases due to the increased number of associations per UAV-BSs. For a fixed number of UAV-BSs (i.e., |U| = 3), the IoT devices at different trajectory way-points tend to utilize the maximum network resource provided by the UAV-BSs. However, since the proposed DQN with replay buffer covers the trajectory way-points more efficiently than that of the baseline approaches, the network resource provided by the UAV-BSs are utilized 9.26% and 4.71% more efficiently compared to the greedy and baseline DQN. Fig. 7 depicts the average energy efficiency of the UAV-BSs that operate under different AoI threshold values. We observe that with relatively relaxed AoI threshold∆ th b = 0.7 the energy efficiency of the UAV-BSs using the proposed approach is 1.63% and 0.95% more energy efficient than that of the greedy and baseline DQN, respectively. The performance gap between the proposed approach and the baseline DQN is relatively close because we use the same discount factor γ = 0.7. However, the performance gap between the proposed approach and the greedy approach is significant over all the threshold values.

D. Discussion
From the experiment results described above, we find some important observations to prove the efficacy of the proposed approach than that of the baseline approaches. The indepth discussion on the experiment results can be summarized as below, • The proposed approach can significantly enhance the UAV-BS trajectory decision where unlike the baseline DQN, the proposed approach can effectively store the transition (i.e., experience) of different environment states to reuse the transition data by random sampling. This stabilizes and improves the DQN training which eventually leads to better trajectory policy that considers the energy consumption of the UAV-BSs, data freshness and bandwidth utilization. • We can observe from Table II and III that the objective of the UAV-BSs is largely dependent on the setting of the appropriate discount factor and the AoI threshold. Especially, the appropriate discount factor value can effectively enhance the training of the DQN network by not only providing better convergence but also giving a chance of improving the training of the DQN network which has a replay memory. • The greedy approach sometimes performs better than that of the proposed DQN with replay memory and baseline DQN. However, the performance gain is limited to a small network and in the case of a large network, the proposed DQN with replay memory significantly outperforms the greedy approach.

VII. CONCLUSION
In this paper, we focused on developing the UAV-BS navigation policy to improve data freshness and accessibility to the IoT network. As a result, we have introduced an agile deep learning reinforcement with an experience replay model that is well-suited to solving the energy-efficient UAV-BS navigation problem under trajectory and AoI constraints. We also performed a comprehensive simulation study to determine appropriate system parameters with the applicable discount factor and AoI efficiency metric threshold to empower the learning model. The simulation results show a strong correlation between energy efficiency and AoI thresholds whereby setting the proper threshold values can effectively enhance the energy efficiency and data freshness for the COC applications. The simulation findings also confirmed the effectiveness of the proposed DQN with experience replay memory under different network conditions.