Abstract
Accurate and rapid estimation of electromechanical mode plays an important role in sensing the security situation of power systems. In this paper, the Compressed Dynamic Mode Decompensation (Compressed-DMD) based estimation approach was proposed to extract the electromechanical mode from high-dimensional ambient data measured by the synchrophasor measurement unit. To improve the efficiency of DMD in processing high-dimensional ambient data under the premise of ensuring calculation accuracy, the Compressed-DMD was introduced to generate the approximation of the high-dimensional left and right singular vectors by employing the aggressive random test matrices and truncated eigendecomposition. Simulation examples of IEEE 16-generator 5-area system and real measurements verify the feasibility and effectiveness of the proposed method.
1. Introduction
Large-capacity long-distance power transmission is an important means to solve the unbalanced spatial distribution of electricity production and demand. Small signal stability is one of the key problems caused by large-capacity long-distance power transmission [1, 2]. Potential weak damping modes of small signal stability seriously threaten the safe operation of the power system [3, 4].
With the access of a high proportion of renewable energy, such as wind power and photovoltaic, etc. The random characteristics of the traditional load random fluctuation characteristics, the access of new energy brings increasingly uncertain factors [5–8]. The limitation of traditional fixed model-based electromechanical oscillation parameter identification method is prominent, and the method based on random signal is attracting increasing attention [9–13]. Electromechanical modal parameters, such as the damping ratio, are important indexes used to assess the small signal stability of large-scale power system. Data driven based approaches for electromechanical mode estimation from synchronous phasor measurements are the effective technology for small signal stability assessment. Among them, synchronized ambient data-based methods have received widespread attention because of its ability to identify the potential weakly damped modes.
Ambient data is a kind of seemingly irregular power system response under natural excitations, such as electrical loads, which vary randomly by nature. It has been proved that the electromechanical modes are contained in the ambient data under natural excitation. Prof. Venkatasubramanian and his group have done fruitful researches in the field of ambient data driven electromechanical mode estimation. The distributed frequency domain optimization (DFDO) [14], frequency domain decomposition (FDD) [15], recursive frequency domain decomposition (RFDD) [16], fast frequency domain decomposition (FFDD) [17], recursive adaptive stochastic subspace identification (RASSI) [18] and orthogonal wavelet bases method [19] were employed to extract the electromechanical modes from the ambient data by Prof. Venkatasubramanian’s researches. These methods have their own advantages in terms of calculation accuracy and speed, etc., and can be applied to different scenarios. For instance, due to the combination of recursive tracking technology, the RFDD and RASSI algorithms have high computational efficiency and are suitable for online applications.
The autoregressive (AR) [20] and the modified AR methods, such as autoregressive moving average (ARMA) [21], Least Mean Square (LMS) [22] and Robust Recursive Least Square (RRLS) [23] based optimal AR model, are also a typical ambient data driven modes estimation approaches. The Bayesian estimation theory was presented to estimate the frequency and damping ratio from ambient data in [24]. The AR, Bayesian method and orthogonal wavelet mentioned above only work well for signal data, consequently, the mode shape cannot be extracted.
System identification technologies represented by stochastic subspace identification (SSI) [20, 25–27] and eigensystem realization algorithm (ERA) [28] were introduced to estimate the electromechanical modes from ambient responses earlier, which is different from traditional algorithm such as Prony algorithm [29]. The classical ERA must be combined with the natural excitation technique (NExT) [30] or the random decrement technique (RDT) [31] because of the requirements of the impulse response as input. As the key step of the SSI algorithm, the singular value decomposition (SVD) (or LQ decomposition) and Kalman filtering are the main sources of the computation time.
Developed more recently, dynamic mode decomposition (DMD) for analyzing dynamical systems is a valuable and effective method that aims at extracting oscillation dynamics from data [32]. DMD was adopted to analyze the dynamic properties of complex power system in recent research. Nevertheless, the capability of the traditional DMD is limited in addressing handling ambient data. In particular, when processing high dimensional data, the computational time efficiency of the traditional DMD is difficult to guarantee.
Nowadays, there are many applications of data processing technology of high dimensional data in power systems such as neural networks and random matrices. Zhang and his group presented a generative adversarial network-based on trinetworks form (tnGAN), which can recover incomplete data by taking full advantage of the peer similarity of data features [33]. This method is effective in data recovery. They also established the spectral distribution change of the established multidimensional matrix to calculate the voltage change and detect hierarchical events based on spectral theory of multidimensional matrix (STMM) [34], which can accurately determine the change of active power and the location of the event source. These methods are useful but tend to classify and detect anomalies. In addition, there are some dimension reduction methods. Wang et al. used reduced second-order closed-loop transfer function to reduce the second-order closed-loop transfer function with the dynamic response of the original inverter, which is reduced to the lowest order [35], which has a high ability to evaluate the stability, maximum overshoot and settling time of the droop-controlled inverter. Then they proposed a reduced-order aggregate model based on a balanced truncation approach and an aggregated approach that involves independent reducing component responses and combining reducing component responses to reduce the input-output maps errors [36]. These methods are significantly efficient in their problem field. However, there is no pre-processing and analysis of random response data under ambient excitation for most methods. The problem solved by the above method is a different type compared with the dimension reduction of random response data with ambient excitation.
In generally, it is recommended that more than one type of electromechanical quantity of inputs from ambient data be used for high accuracy and reliability of the estimation modes. The data dimensions of the input of DMD will rapidly increase as the increase of types of the selected electromechanical quantities, the ensuing result is a marked reduction in computational efficiency. Therefore, the compressed dynamic mode decomposition (compressed-DMD) based mode estimation method, in which multiple types of electromechanical quantities containing the generator speed, power angle and active power are selected as the input data, was presented to improve the computational efficiency. The main contributions of this paper are as follows:(1)A new random test matrix was introduced to reduce the dimension of the input ambient data through multiplying by the original input data. Then, the compressed-DMD, in which the traditional singular value decomposition (SVD) was replaced by a compressed SVD, was proposed to ensure the accuracy of the estimation of the electromechanical mode after dimensionality reduction for designing the further steps.(2)A novel scheme based on the compressed-DMD is developed to extract the electromechanical oscillation modes from ambient data. In the proposed compressed-DMD based scheme, the ambient modes can be estimated highly efficient and precise from the high-dimensional multiple types input data.(3)The estimation results of multiple combinations of different types of electromechanical quantities containing the generator speed, power angle and active power as input were analyzed using the ambient data simulating on IEEE 16-generator 5-area test system, firstly. Next, the benefits of the proposed compressed-DMD based approach are illustrated using simulated data and a real measurement.
The remaining structure of this paper takes the form of three sections including the introduction section. Section 2 covers the establishment of review of traditional dynamic mode decompensation algorithm and presents the Compressed-DMD algorithm. The simulation case of IEEE 16-generator 5-area system and the real measurement case are presented in Section 3.
2. Materials and Methods
2.1. Mathematical Analysis of Dynamic Response of Power System under Ambient Excitation and Dynamic Mode Decompensation for Ambient Data
The dynamic model of the power system can be described by stochastic differential-algebraic equation:Where equation (1) is the system differential equation, which describes the dynamic process of the generator and corresponding control device in the system. Equation (2) is the algebraic equation of the system, which is generally composed of the system power flow equation and static equations such as generator and load. f and are continuous functions; x is the state variable of the system, such as generator power angle, etc. y is the algebraic variable of the system, such as bus voltage, bus phase Angle, etc. u is the random variable volatility variable.
Assuming that the random fluctuation obeys OU distribution, the system dynamic model of random fluctuation is shown in equation (3):Where, C is the diagonal element of the diagonal matrix by , and is the correlation time of random fluctuation process; δ is system noise intensity; ξ is an independent Gaussian random distribution vector.
Equations (1)–(3) together constitute Stochastic Differential Algebraic Equation (SDAE) of the power system under random load fluctuation. In the process of random disturbance, if the operating state of the power system does not change, that is, the operating equilibrium point of each generator remains unchanged, the dynamic process of the system state variables can be approximately analyzed by linearizing the system state. Linearization of equations (1)–(3) and elimination of algebraic variable y can be obtained, which is shown in equation (4).Where, and are the Jacobian matrices corresponding to x and y in equation (3) respectively. is the identity matrix. , and are the Jacobian matrices corresponding to x, y and u in equation (2) respectively. We choose:
Then equation (4) can be represented below:
Equation (8) is the linearized state-space model of the power system under random fluctuation. The eigenvalues corresponding to n oscillation modes of the system can be solved, then the analytical expression of the system state variables in the time domain is:Where, and ui are the right and left eigenvectors corresponding to the eigenvalue λi respectively; z(0) is the initial value of the state variable.
Equation (9) shows that the mathematical analytical expression of the random response of the power system under ambient excitation is composed of two parts: the first part is the oscillating component which contains the characteristic information of electromechanical oscillation of the system (oscillation frequency, damping ratio, etc.), and the second part is the random component which approximately obeys the Gaussian distribution. Therefore, the seemingly random response of the power system contains abundant mechanical and electrical dynamic characteristics in nature.
During the operation of the power system, natural excitations, such as small fluctuations in load and renewable energy, are always there. Responses of the power system under natural excitation are defined as ambient responses. The ambient data X = [x(t1), x(t2), …, x(tm)] discretized time domain ambient response can be collected by employing the PMU (Phasor Measurement Unit). Although the ambient data manifests itself as disorganized and is characterized as random Gaussian white noise, the electromechanical modes are contained. For instance, Figure 1(a) shows the ambient data of the active power collected from a real power system. Even though it does not show oscillations in Figure 1(a), the system has an obvious peak near 0.32 Hz with the power spectral density shown in Figure 1(b), which indicates an obvious electromechanical mode with the frequency of 0.32 Hz.

(a)

(b)
The dynamic mode decomposition (DMD) technique [37] is a powerful data-driven tool that extracts modes from measurements. After collecting the ambient data, the modes containing oscillation frequency, damping ratio and mode shapes can be estimated by using DMD as follows:
Forming the snapshot data matrix and :
Calculating the unitary matrices and through the singular value decomposition (SVD) of the original data matrix :where and satisfy ,; I is a unit matrix; Sr is a diagonal matrix with singular values.
Then the low-dimensional approximation matrix can be obtained by:
The eigenvalues and corresponding eigenvectors with the oscillation behavior of the dynamic system can be determined by the eigen-decomposition of the low-dimensional approximation matrix, by using formulation .
The DMD eigenvectors of the eigenvalues are given by the following formula, which is DMD mode.
The traditional DMD procedure is demonstrated in Algorithm 1. As can be seen, singular value decomposition is the core step in DMD. Nevertheless, the traditional SVD is very time-consuming for the high dimensional ambient data.
|
2.2. Ambient Modes Extraction Based on Compressed Dynamic Mode Decomposition
Although the traditional DMD algorithm has good computing power, it still has some challenges to deal with the calculation of big data matrix. Meanwhile, the truncated singular value decomposition algorithm in traditional DMD algorithm and other approximate singular value decomposition algorithms have the only disadvantage that they will produce error penalty factors.
For the high-dimensional input data, the computational efficiency on the premise of ensuring the computational accuracy is essential for the estimation electromechanical modes from synchronized ambient data of large-scale power system. In this section, to be computation efficient, the compressed-DMD will be implemented. Firstly, a random test matrix is employed to approximatively construct the compressed data to replace the high-dimensional input data. The purpose of matrix compression is to obtain a small compressed matrix that contains most of the information of the original data matrix. Similarity and difference can be determined through the geometric vector relationship between several data points [38]. Compressed matrix can ensure that the distance and angle between data points are preserved. The restricted isometric property (RIP) guarantees this property with high probability with a sufficient number of data points [39]. Next, we use the Compressed Dynamic Mode Decomposition (Compressed-DMD) to generate the singular vectors which are the key variables for forming the low-dimensional approximation matrix in the traditional DMD with multiple channels data.
2.2.1. Compressed Dynamic Mode Decomposition
To construct an effective measurement matrix, we rely on randomness. The most prominent choice is the Gaussian random test matrix, which consists of independent identically distributed (iid) standard normal terms. However, there are still some disadvantages Gaussian measurement matrices, such as high expensive cost for generating and sensing matrix multiplications, which are computationally intensive. In this paper, we use a compressed sparse row matrix with only a few nonzero iid random entries is constructed (test matrix is shown in Figure 2). which is shown below:Where:

is the parameter to control the sparsity. This paper uses a high aggressive sampling rate to get accurate results, which allows rapid execution of the compressed process. Therefore, the columns or rows are likely to be linearly independent.
At the beginning of the compressed-DMD, we provide a random test matrix to get a low-rank SVD approximation of as well as target rank , where and, represents a parameter for oversampling. For the oversampling parameter, we could observe a small amount of oversampling, which is usually sufficient in actual utilization. The compressed data matrix can be obtained as:
The approximate right singular vector is the final purpose, in general, we usually calculate the singular values and vectors from the inner and outer dot product
can be calculated by , which is singular, then we continue to obtain a small matrix :
To make sure the symmetry of the matrix , we need:
Then compute the eigen decomposition:
The approximate eigenvectors and eigenvalues are represented by the matrix and . Therefore, the approximate singular value is . Then, we truncate the decomposition if , which means that the first dominant eigenvectors and singular values can be extracted to be updated our prior decision of the target rank based on the approximated singular value spectrum. Thus, it is sufficient to use a smaller target rank because of its flexibility.
The approximate right singular vectors can be defined as follows
And the approximate left singular vectors are similarly defined as:
The idea for these two formulas comes from equation (12), which are and .
Finally, a second pass over the high-dimensional input matrix is involved. We often need to provide an additional updating step, because the columns of are only approximately orthogonal due to the approximation errors in the eigenvalues. Significantly, the SVD of the scaled right singular vectors (principal components ) can be calculated as:
Then, we continue to update the right singular vectors:
All steps are the whole approximate computation of the singular dynamic mode decomposition, which is shown in Algorithm 1. In addition, the main idea of this approach is motivated as equation (26), which is shown in Figure 2 at the same time. Then we can continue to replace equation (12) with equation (26) and perform the rest of the dynamic mode decompensation.
Then, define the same low-dimensional approximation matrix:
Compute the eigenvalues and eigenvectors of matrix , writing:
The eigenvalues and eigenvectors obtained from compressed-DMD are efficient to identify the electromechanical oscillation parameters. It can greatly improve the calculation speed and ensure the calculation accuracy. In addition, the whole process of compressed-DMD algorithm is shown in Algorithm 2.
Remark 1. In practice, the process of singular value decomposition is to decompose the matrix into two orthogonal matrices and a diagonal matrix. There is no need to do SVD for the orthogonal matrices, which can be written directly. Because at the beginning of SVD, and are identity matrices, the eigenvalues are 1. Thus, the of SVD is directly, the other two matrices are also identity matrices. In actual calculation, each column of matrix is approximately orthogonal with approximate errors of eigenvalue calculation. Therefore, even though the matrix decomposition cannot be written down directly, it is exceptionally fast for singular value decomposition of .
|
2.2.2. Ambient Modes Extraction
Using Compressed-DMD algorithm, some descriptions are presented to provide oscillation temporal characteristics (oscillation frequency , damping ratio and spatial characteristics (mode shape) from the measurement data. For example, we extract the eigenvalues , which are transformed from the measured data into a continuous system mode expression.
Then calculating frequency and damping ratio of the mode, which are:
The mode matrix is defined as
Which is corresponding dynamic mode .
3. Application of the Proposed Method on Mode Estimation for Power System
3.1. IEEE 16-Generator 5-Area Test System
Taking the IEEE 16-generator 5-area system as an example to illustrate the online applicability of the Compressed-DMD algorithm. The single-line diagram of the IEEE 16-generator 5-area system is shown in Figure 3 [40].

In the basic operation mode, the calculation results of the electromechanical small interference analysis show that the system has 15 electromechanical oscillation modes, the damping ratios of the small signal stability analysis (SSSA) for four inter area oscillation modes are: 10.35%, 7.25%, 2.80% and 6.99%, which includes one weakly damping mode (mode 3), showing in Table 1.
We use the generator speed, power angle and active power of each generator as the input data set of the ambient excitation. Based on the 400 s ambient excitation data of 16 generators of the system as an example, using the traditional DMD algorithm to calculate the frequency and damping ratio, one channel, two channels and three channels of input signal data were selected to carry out the identification calculation, the identification results of Compressed-DMD are used as the comparison with three channel data input, which is shown in Table 2, Table 3 shows the single calculating time of five situations.Where:
The input of Case 1 is Generator speed.
The input of Case 2 is the combination of Generator speed and Generator active power.
The input of Case 3 is the combination of Generator speed, Generator power angle, Generator active power.
With the increasing data channels are entered, the accuracy of DMD algorithm is getting higher and higher. Using Mode 3 as an example, if only angular velocity is used as input, the identification result is the normal damping mode, but the weak damping mode cannot be identified, which is 6.65%. When the input increases to the combination of the generator speed, power angle and active power, the accuracy of the identification results is greatly improved, which is close to the real value. But the corresponding calculation time has also increased significantly. High-performance workstations are needed to drive computing. Therefore, Compressed-DMD algorithm can be used to reduce the calculation time of the DMD algorithm with three sets of data, the calculation accuracy and efficiency of the two algorithms were compared. At the same time, a mature algorithm, SSI, was compared.
To make sure the accuracy and reliability of Compressed-DMD algorithm, we need to use the Monte Carlo method to make validation. Therefore, we choose the sliding window method to calculate the mean value and variance. Firstly, we need to choose the appropriate window length from two sides for Compressed-DMD algorithm and contrastive algorithm (Traditional-DMD algorithm and SSI-algorithm), which are the calculation time and accuracy.
On the one hand, for the IEEE 16-generator 5-area system, through the above analysis, we use the combination of generator speed, power angle and active power as input. At the same time, the traditional-DMD algorithm and the SSI algorithm are compared with the compressed-DMD algorithm. Using Lenovo Tianyi workstation to identify the oscillation parameters with different time window lengths. Figure 4 shows the single calculation time of window sliding once for the three identification methods.

From Figure 4, we can observe that the calculation time changes with different selection of window lengths. Compressed-DMD has better computing speed compared with SSI algorithm and traditional DMD. When the data segment length is 50 s, the computing time is 0.0138 s for Compressed-DMD algorithm. On the contrary, the SSI algorithm needs 0.2286 s. When the data segment length increases, the three algorithms need more time to calculate, but Compressed-DMD still need least time. Although the calculation time of Compressed-DMD algorithm is close to that of the traditional DMD algorithm, Compressed-DMD algorithm is still faster and closer to real-time calculation than the traditional DMD algorithm and SSI algorithm, as shown in Figure 4.
On the other hand, calculation platform performance and engineering requirements need to select a reasonable sliding calculation window as the input data window. To verify the relationship between data length and extraction results, we select data segments of different lengths to analyze and identify the accuracy of the results relationship. It can be seen that the extraction results of Compressed-DMD algorithm and traditional DMD algorithm have less fluctuation with the increase of window length compared with SSI algorithm, which means they have better robustness, using Mode 3 as an example, which is shown in Figure 5. It can be seen that the longer the window length is, the higher the identification accuracy of the three algorithms is.

Secondly, we need to verify the accuracy of the Compressed-DMD algorithm compared to the other two methods using Monte Carlo calculation. The larger the window length we choose, the results of the three methods will be more accurate and more comparable. However, considering the calculation speed, we compromise the window length of 400 seconds. The sliding speed is 10 seconds, sliding 1000 times, the total calculating time and identification results are shown in Table 4 and Table 5.
Table 4 presents that the total calculation time and speed are more obvious, Compressed-DMD with Case 3 (three channels) cost less total calculating time with the other two algorithms, which is close to total calculating time of DMD (Case 1). The data in Table 5 shows that the mean values of the four interval oscillation modes of the system and the damping ratio obtained by the proposed Compressed-DMD extraction are relatively close to the SSSA results. Firstly, the mean frequency of Compressed-DMD algorithm is closer to the SSSA. Using Mode 3 (weakly damping mode) as an example, the mean values of frequency are similar for these three algorithms, which is 0.5592 Hz, 0.5513 Hz and 0.5801 Hz respectively, but the result of Compressed-DMD algorithm is closer to the result of SSSA and the Compressed-DMD algorithm has smaller variance (). In addition, the maximum variance of the oscillation frequency is and the variance of the damping ratio is 0.1098, which indicates that the extraction results of Compressed-DMD of four modes fluctuate slightly near the mean, which ensures the reliability of the identification results. What’s more, the mean values of the four modes of SSI algorithm and traditional DMD algorithm have some deviations because of an increasing number of outliers. And the variance of the SSI algorithm results is significantly greater than that of the Compressed-DMD algorithm, and the reliability of the results is lower than that of the Compressed-DMD algorithm. At the same time, Figure 6 shows that the modal shape figures of the four modes are consistent with the results of SSSA, which further verifies the applicability of Compressed-DMD algorithm.

Figure 7 shows the Monte Carlo histogram for 4 mode damping ratios for the Compressed-DMD algorithm. It can be seen that the identification result of damping ratio presents a similar normal distribution near the mean value, which can prove the accuracy of the extraction result of damping ratio.

3.2. Real Measurements
The second case analyzes the performance of the Compressed-DMD in the actual measurement system for some actual system, which is recorded from PMUs during a real period normal operating resource from [41]. Damping ratio and frequency estimation results are achieved by analyzing using Compressed-DMD method and SSI method are presented in this section. Here another quantity is analyzed: Power Density Spectrum. The figure of Power Density Spectrum is shown in Figure 8.

Figure 8 presents the power spectrum of the selected bus frequency. The power spectra provide a critical mode at approximately 0.30 Hz, which means that the mode near 0.30 Hz can be observed strongly for all PMUs by the analysis of the behaviour. Another peak value can be observed between 0.8 Hz and 1 Hz.
The results of damping estimation and frequency estimation for the Compressed-DMD and SSI method are presented in Table 6. The analyzing time window length is 400 s, and the speed are chosen 10 s.
The damping ratio and frequency can be identified by the two methods respectively for real measurement system, the frequency of Compressed-DMD method for Mode 1 is 0.2883 and Damping ratio is 7.23%. The frequency is near the peak value of frequency in tie line power. Mode 2 shows another identified result corresponding to the second peak value, the frequency of Compressed-DMD method for Mode 2 is 0.8665, and Damping ratio is 7.48%. The variances of frequency and damping ratio for both two modes are very small, which means the Compressed-DMD method has better stability and can accurately identify real measured data. In addition, Table 7 shows that total calculation efficiency of the Compressed-DMD is better than the SSI algorithm.
However, the statistical average of the frequency identified results of SSI method is 0.2787 Hz and 0.9008 Hz. Although the statistical mean is also close to the true value, the variance is large, indicating that the value calculated by the sliding window fluctuates greatly. At the same time. The identified results of damping ratio from Table 6 are 15.63% and 10.95%. The statistical mean error rate of the identification result of the damping ratio is very large because it contains many outliers, which affect the result of the mean value. Therefore, the identification result of the damping ratio of the actual measured data by SSI is not reliable compared with Compressed-DMD method.
By comparing with the SSI method, the mean values of frequency and damping ratio are more accurate respectively with less variance in the results of Compressed-DMD method. Figure 9 shows the scatterplot and histogram of frequency and damping ratio. However, there are some high fluctuations in SSI method. Therefore, Compressed-DMD is more stable and has better robustness in real actual data.

The spatial characteristics of electromechanical oscillation for actual measured data can be described using the mode shape. After identification of the Compressed-DMD method, the eigenvector as the corresponding spatial characteristics can be determined. As shown in Figure 10, the obvious interarea swing pattern can be observed for the four areas. Therefore, the mode shape shows the PMUs behaviours, which are consistent with the results from the power density spectrum.

4. Conclusion
To eliminate the potential weak damping mode of the system, a new method of online damping monitoring based on ambient data is proposed. This paper proposes a combination of compressed SVD and traditional DMD, using a compressed sparse row matrix to multiply the input data for compression. The purpose of Compressed SVD is to ensure that the information of the input matrix after compression is more accurate. The function of the whole Compressed-DMD algorithm is to reduce the complexity of data and ensure the accuracy of calculation with multi-signal data input. It can greatly reduce the identification calculation time of traditional DMD algorithm, which can meet the requirements of DMD algorithm for multi-signal channel calculation.
Through this method, the Compressed-DMD algorithm is superior to the traditional DMD algorithm and SSI algorithm in terms of calculation speed and accuracy rate of extraction results. Compressed-DMD using multi-channels ambient data input ensures the high-speed calculation ability and high accuracy at the same time. In addition, the Compressed-DMD algorithm has good robustness and real-time performance. Finding the potential weakly damped oscillation mode with high-dimensional data in the system provides a stable and accurate monitoring platform for dispatch operators. The next step is to adopt a targeted control strategy based on the identification results based entirely on multi-signal input data.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported in part by the National Key Research and Development Program 2021YFB2400800.