A Levitation Condition Awareness Architecture for Low-Speed Maglev Train Based on Data-Driven Random Matrix Analysis

Modern low-speed maglev trains typically use multi-node decentralized levitation control modules, which results in a complex levitation control system with coupling interaction. Conducting systematic levitation condition awareness of the levitation control system is still a promising challenge. In this paper, under the hypothesis of levitation residuals following normal distribution, a levitation condition awareness architecture for the levitation control system is proposed based on data-driven random matrix analysis. The proposed architecture consists of an engineering procedure followed by a cascaded mathematical procedure. In the decentralized engineering procedure, the data-driven modeling for individual levitation control modules is achieved by nonlinear autoregressive modeling with an exogenous input neural network, and the unknown parameters are identified by a modified combinatorial genetic algorithm. On this basis, high-dimensional analysis of streaming residual random matrices for the levitation control system is conducted aided by large-dimensional random matrix theory, and the control limits of the constructed indicators are well-designed using the theorical distributions. Based on the comparative analysis of the experimental datasets, the proposed awareness architecture is verified to show the effectiveness of the systematic condition evaluation of the levitation system, and incipient train-guideway interaction vibration abnormalities can be detected in a timely manner.


I. INTRODUCTION
As a promising transportation technology, low-speed maglev train has the advantages of low vibration and noise, no friction loss, environmental friendliness, and low maintenance costs [1]- [5]. In the past decades, several studies and demonstration projects have been carried out on medium-low speed maglev train, such as the Korean Urban Transit Maglev (UTM) [3], and the Japanese High-speed Surface Transport (HSST) [4]. Changsha Maglev Express (CME), in particular, which was officially put into operation The associate editor coordinating the review of this manuscript and approving it for publication was Eunil Park . in 2016 in China, is the longest low-speed maglev commercial demonstration line in the world [5].
The levitation system is critical to long-term reliable and stable operation of maglev trains. For low-speed maglev trains, electromagnetic suspension (EMS) is utilized for levitation. EMS achieves train levitation through the electromagnetic attraction between the track and the electromagnet on the bottom of the train.
Taking the CME as an example, a typical structure of such a low-speed maglev train, which is composed of three carriages, is presented in Fig. 1. Each carriage is equipped with five maglev frames and twenty independent levitation control modules (LCMs) with ten on each side. Fig. 2 is VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/  the schematic of an individual LCM. The integrated sensor located at the bottom of the maglev frame collects the levitation gap and acceleration data. On this basis, the levitation controller controls the output current of the chopper according to the levitation gap data, so that the electromagnetic attraction generated by the suspension electromagnetic coils can ensure that the levitation gap is within a preset range [8], [9]. LCMs fulfill basic requirements for maintaining the expected levitation gap against disturbances from the train body vibration and track unevenness during train operation through the use of advanced control strategies [10]- [13], which results in a complex levitation control system with coupling interaction. On top of these hardware subsystems, an effective online condition monitoring (CM) system enables real-time levitation condition awareness of the maglev train. Information such as the floating clearance, vertical acceleration, and electromagnet current are uploaded to the CM system, which directly reflects the levitation of the train. Contemporary CM systems, which utilize the uploaded information from these control modules, can alarm the operation abnormality using the predefined levitation gap threshold. Nevertheless, the contemporary CM systems fail to provide timely warning information about train levitation condition against incipient faults. As the levitation control system (LCS) is a typical multi-node complex close-loop system, if the incipient faults cannot be detected or tracked in a timely manner and addressed effectively during the maglev train normal operation, they will accumulate continuously and propagate inside the integrated system through feedback control, thereby triggering chain reactions that cause abnormal system shutdown or damage [14], [15]. Therefore, operating condition awareness is a significant practical merit to enhance the reliability and durability of the LCS, which is one critical element of maglev train.
Fortunately, driven by the development of big data technology [16], [17], data correlations indicated by statistical parameters in multi-node complex control system can be computed, and insight into inherent mechanisms can be attained. As an effective universal data analysis tool, random matrix theory (RMT) can be applied to the operating condition awareness of the decoupled multi-node LCS. RMT has been introduced in the big data architecture design of smart girds for correlation analysis and anomaly detection purposes [18]- [20].
The basic principle of the original large-dimensional random matrix theory is that, if the elements of a random matrix that characterize the operating state of a system satisfy the independent identical distribution (i.i.d) of Gaussian nature, the eigenvalue distribution of the corresponding covariance matrix will exhibit nonrandom characteristics [18]. Then, the variation of the eigenvalue distribution characteristics of the random matrix resulting from system abnormality can be utilized to achieve system condition evaluation and abnormality/potential fault detection.
However, the original random matrices characterizing actual operating condition of LCS generally cannot satisfy the hypothesis of strict i.i.d Gaussian distribution, as the element channels formulating the streaming random matrices have inherent spatiotemporal correlations. To eliminate the spatiotemporal correlations of streaming random matrices, a modified random matrix analysis theory based on ARMA modeling has been proposed in recent years. This theory allows obtainment of the analytical solution to the eigenvalue distribution of large-dimensional random matrices that are constructed with autocorrelation random variables [22]- [25].
Nevertheless, individual LCMs can be considered as a typical multi-input-single-output (MISO) time-varying dynamic system. The dynamics of the associated controlled levitation gap not only have distinct autocorrelation characteristics, but are also closely related to the input current controlled by the levitation chopper. Additionally, the neighborhooddistributed LCMs show complicated coupling vibration characteristics. This indicates that there is complexity in the spatiotemporal correlations among multi-point levitation gaps that can characterize the operating condition of LCS. Quantifying the spatiotemporal correlations among the elements of the characteristic random matrices remains difficult with the existing large-dimensional random matrix analysis methods. As a result, different operating conditions (normal/minor fault) cannot be distinguished effectively, and an associated system operating condition awareness cannot be achieved.
The dynamic relationship between the suspended air gap and the controlled input current can be partially characterized by the state-space equations-based mechanism modeling of the suspension control module. On this basis, an iterative estimation of the levitation gap can be obtained by employing model-driven methods, and the streaming random matrices can be constructed using the multi-node levitation gap estimation residuals to eliminate the inherent influence of the spatiotemporal correlations. However, it is difficult to implement the analysis of levitation gap correlation based on the mechanism modeling of LCMs. On the one hand, the coupling vibration characteristics of neighborhood LCMs are not accounted for in the mechanism models, and additional residuals can be introduced by the unmodeled dynamics. On the other hand, the levitation gap control is based on the fast-dynamic modulation of the levitation chopper, and the frequency of the levitation control is usually 10000 Hz, while the upload frequency of the status monitoring information is merely 10 Hz [29]. The expected signal sampling boost for the levitation gap estimation will result in significant computational overhead to the configured low-cost levitation controller. Moreover, an additional communication load can also be imposed, which is detrimental to the stable operation of the existing monitoring based control architecture.
In contrast, the data-driven system modeling and identification method can not only mine the time-dependent dynamic correlation between input and output of the system effectively but also reduce the transient impact of the fast-dynamic control process. This method is more suitable for the correlation identification analysis of the LCS modeling.
To effectively characterize the complex spatiotemporal correlation between the inputs and outputs of a dynamic system, the autoregressive correlation characteristics between the input and output variables are further considered in some studies based on the identification of nonlinear neural networks [30], [31]. This type of black-box models, or so-called data-driven system identification models, can better approximate real dynamic systems, without the necessity of knowing all the physical parameters or coefficients. Therefore, it is promising that, the abnormality of individual LCM can be reflected by the model residuals, which remains an open problem that needs to be resolved.
All things considered, a novel systematic condition awareness architecture for the LCS of low-speed maglev train is proposed in this paper. The major contributions of this paper are summarized as follows: 1) A nonlinear autoregressive with exogenous input neural network modeling method is proposed to determine the levitation gap residual generation for individual LCMs. On this basis, the Gaussian distribution characteristics of streaming data matrices under normal system operation are enhanced, while the signal correlation feature can be retained under abnormal condition.
2) A modified combinatorial genetic algorithm (GA) is proposed for the parameter identification of decentralized LCM modeling. The structure parameters of the neural network-based LCM model, which are the number of neurons in hidden layer and the weights-bias parameters of neural network, can be optimized simultaneously using the proposed algorithm. This allows deep mining of inherent dynamic correlation characteristics of the decentralized LCMs.
3) Aiming to address the difficulty of analytically solving the random characteristics of the large-dimensional time series with complex spatiotemporal correlation, a twostep decoupling analysis framework is proposed, which is based on the decentralized identification of LCMs dynamics and the centralized large-dimensional residual correlation analysis of the LCS. On this basis, systematic condition awareness of multi-node coupled LCS for low-speed maglev train is achieved for the first time, which can further improve the state-aware performance for incipient train-guideway interaction vibration abnormality.
The reminder of this paper is organized as follows. Section II describes the experiment based on the CME and the associated data acquisition system is described, as well as the experiment carried out to collect the data for the architecture implementation. In Section III, the proposed levitation condition awareness architecture for low-speed maglev train is described, the mathematical fundamentals of the RMT are introduced, and the modeling methodology using GA based NARX neural network structure is presented. A brief statistical analysis of the experimental data and the results obtained with the proposed architecture combining centralized RMT evaluation with decentralized levitation gap residual generation are presented in Section IV. Finally, in Section V, conclusion is drawn, highlighting the main contributions of this paper.

II. EXPERIMENTAL SETUP AND DATA ACQUISITION
The CME stretches over 18.55 km and serves Changsha International Airport, Langli, and Changsha South Railway Station. The maglev train operational experiments and data analysis were carried out on this line. The train has three carriages: the head car (MC1), the middle car (M), and the tail car (MC2), as shown in Fig. 1.
In the experiment, the monitoring signals include the levitation gaps, the vertical accelerations and the suspension electromagnet currents of all LCMs, which characterize the levitation condition of each levitation point. All data are uploaded to the CM system and stored. The maximum speed of the train is 100 km/h, and the expected levitation gap is set at 8.5 mm.
The levitation gap and the levitation acceleration are sampled by the integrated sensor located at the bottom of the train, while the suspension electromagnet current is sampled by the chopper current transformer.
Due to the distribution of the train structure, the levitation coupling between the three cars is weak, while neighborhood LCMs located at the same car have strong vibration coupling characteristics. Therefore, it is possible to assess the levitation condition of each car separately. The middle carriage M, which is equipped with 20 LCMs, was chosen for this study, and two pieces of data with durations of 1500s and 1200s were selected from all the test data. The data from the No. 5 LCM during the 1500s test period are shown in Fig. 3. In this period, all the stages of maglev train operation are included, train way-in and way-out, acceleration, deceleration, and static levitation.
In Fig. 3, the levitation gap fluctuates from 6 mm to 12 mm, and the vertical acceleration fluctuates approximately 10 m/s 2 under normal condition. It is noteworthy that, the train speed is zero during the interval of [900s, 1000s], indicating that the train enters the station while remaining levitated. At this time, the train is less affected by external disturbances, and the levitation gap can be maintained preferably at approximately 8.5 mm as expected, with small vibration acceleration and electromagnetic current fluctuations. However, this is not the case for the interval of [150s, 500s], during which the train is also in the same levitated condition with zero speed. At the beginning, the fluctuations of levitation gap, vertical acceleration and current are also small, while the levitation gap becomes larger after a few seconds, and obvious abnormal fluctuation can be observed. The fluctuation of the vibration acceleration varies between −50 m/s 2 and 40 m/s 2 and is much more severe than that during normal operation. This suggests that the train is in the train-guideway interaction vibration abnormal condition during this interval, but the alarm is not triggered in the CM system since the levitation gap does not exceed the alarm threshold and the levitation acceleration has not been configured as an abnormal diagnosis indicator in the CM system.
It should be noted that, the levitation abnormality represented by the train-guideway interaction vibration can be considered to be incipient faults for the low-speed maglev train and will lead to more serious anomalies if it is not detected in a timely manner. However, according to the above analysis, the configured levitation gap threshold-based diagnosis is far from competent, and there is an urgent need for a novel levitation condition awareness architecture with enhanced robustness.

III. LEVITATION CONDITION AWARENESS ARCHITECTUE
The proposed levitation condition awareness architecture for low-speed maglev train is illustrated in Fig. 4.
This system is comprised of two parallel processing procedures. One procedure is continuous random matrix modeling that can be considered as an engineering procedure based on decentralized data-driven approaches. The other procedure is a formulated random matrix analysis, which can be regarded as a mathematical procedure for achieving systematic levitation condition awareness for the maglev train. During the engineering procedure, critical input and output signals of each LCM are sampled, and the levitation gap residuals are obtained from the decentralized data-driven NARX neural network (NARX-NN) modeling. During the RMT-based mathematical procedure, meticulous design of the effective eigenvalue indicators is needed to reflect the system operating condition. It is noteworthy that, the functions of the anomaly diagnosis, location, and fault-tolerant control are included to complete the proposed architecture, as shown in Fig. 4. The systematic condition awareness-based functions shown with a dashed box, however, are beyond the scope of this paper.
This system is comprised of two parallel processing procedures. One procedure is continuous random matrix modeling that can be considered as an engineering procedure based on decentralized data-driven approaches. The other procedure is a formulated random matrix analysis, which can be regarded as a mathematical procedure for achieving systematic levitation condition awareness for the maglev train. During the engineering procedure, critical input and output signals of each LCM are sampled, and the levitation gap residuals are obtained from the decentralized data-driven NARX neural network (NARX-NN) modeling. During the RMT-based mathematical procedure, meticulous design of the effective eigenvalue indicators is needed to reflect the system operating condition. It is noteworthy that, the functions of the anomaly diagnosis, location, and fault-tolerant control are included to complete the proposed architecture, as shown in Fig. 4. The systematic condition awareness-based functions shown with a dashed box, however, are beyond the scope of this paper.
In the following subsections, the data-driven random matrix modeling based on the NARX-NN structure for the levitation gap residual signals from decentralized LCM is described, and the RMT-based condition awareness method is explained in detail.

A. MODELING METHODOLOGY
Since the NARX structure has feedback connections, it can be used to build autoregressive models, i.e. models whose current output depends on past values of the output itself. The model presented in this paper is based on the idea that the changes in the LCM conditions can be reflected by the levitation gap itself. Hence, it is intuitively reasonable to apply an autoregressive model employing the past values of the levitation gap as part of the model's input. In addition, the input also includes past and present signals of the control current and acceleration of vibration.
As shown in Fig. 5, the NARX-NN structure provides a prediction y(n) of the output based on the past values of the output and on the past and present values of an external input, where the time delay (TDL) represents the time delay series. In this NARX-NN model, the external input vector X is composed of two elements, namely, the levitation control current and the vibration acceleration.
In matrix notation, the two-layer NARX-NN model can be integrated represented aŝ where y denotes the iterated predicted stack voltage; ϕ T ∈ R 1×(m+p−1) represents the regression vector formulated by measured current, past inputs, and past outputs; θ is the parameter vector corresponding to bias {b H , b o } and weights {W Hx , W Hy , W Ho } of the NARX-NN model; and f (·) is the union of the nonlinear functions of the hidden-layer f H (·) and output-layer f o (·) determined by a feedforward neural network (FNN) with r neurons.
It is noteworthy that the model signal correlation involved in the prediction residuals can be minimized if the system dynamics under normal operating condition can be well characterized by the parameter vector θ. Accordingly, the additional signal correlation resulting from an abnormal condition can be reserved in the residuals and further identified. Therefore, model identification is essential.
The identification procedure can be initiated when the experimental dataset is ready. It is divided into the following two stages: Minimum time delay identification, and model training and validation.

1) MINIMUM TIME DELAY IDENTIFICATION
As a first step, the number of minimum TDL needs to be identified based on the prepared experimental dataset (or real-time sampled dataset in the case of on-line monitoring). This provides reference for the determination of the actual TDL applied. The sampled levitation gap signal can be considered stationary as it is controlled to track the reference levitation gap. At this point, the autocorrelation function (ACF) and the partial ACF (PACF) can be applied to evaluate the dependence of the present output on the past samples. Similarly, the cross-correlation function (XCF) can be applied to evaluate the influence of the past input series samples on the present output [32].
Based on the ACF, PACF, and XCF results, it can be confirmed that the lags in the correlation plot with higher lengths of line segments correspond to the samples with stronger influence on the present output. Hence, the minimum TDLs of the model input and output can be identified straightforwardly.
It is noteworthy that the identified minimum TDLs are references for the selection of the actual TDLs. Considering that the model is built into LCMs with upload rate of 10 Hz, the number of delayed samples should be a tradeoff between the estimated residuals and the processing time, which requires optimization.

2) MODEL TRAINING AND VALIDATION
Generally, the model training and validation stage for each NARX-NN model corresponding to an individual LCM are carried out off line using a typical test dataset. The test dataset is classified into a training dataset and a validation dataset.
In the training sub-stage, the training dataset is selected from the prepared test dataset. Training is performed with the series-parallel architecture using NARX, so that input to the feedforward network is more accurate and that the resulting network has a purely feedforward architecture for which static backpropagation can be used. Furthermore, the validation sub-stage is carried out with the trained NARX network by applying the validation dataset.
Once the basic structure of the NARX-NN-based model is fixed as shown in Fig. 5, it is necessary to identify the structure parameters {m, p, r} and to optimize the weight-bias parameters {W Hx , W Hy , W Ho , b H , b O } of the FNN, to minimize both the testing root-mean-squared error (RMSE), and the structure complexity. This helps guarantee modeling accuracy and generalization capability with a compact NARX-NN structure.
Although traditional backpropagation error algorithms, such as Levenberg-Marquardt (LM) and Scaled Conjugate Gradient (SCG), are suitable for function approximation problems [33], they can only optimize the weight-bias parameters of the FNN with predefined structure parameters based on a priori knowledge or exhaustive searching [32]. In contrast, evolutionary algorithms, especially GA, are preferable for the mixed integer optimization with integer structure parameters and fractional weight-bias parameters [34], [35]. Hence, an improved GA algorithm combined with backpropagation is proposed to simultaneously optimize the structure and weight-bias parameters. The optimization process of the combinatorial GA algorithm is modified based on the improved GA algorithm proposed in [35]. The differences between the two GA algorithms are described as follows.
In the encoding/decoding stage, the encoding of the i th chromosome is designed as follows: where N is the population size; C i is a 3-by-6 matrix, in are encoded by binary number c M 1 · · · c M 6 , c P 1 · · · c P 6 , c R 1 · · · c R 6 ; N m , N p , and N r are empirical values of the maximal lags of X(k), y(k) and the number of neurons in hidden layer, respectively. The decoding of m, p and r is then derived as follows: where · rounds the element to the nearest integer. After decoding each chromosome to obtain structure parameters of the NARX-NN model, weight-bias parameters can be calculated based on the backpropagation error algorithm in terms of the training dataset.
In the fitness evaluation, the fitness function minimizing both the modeling errors of the training dataset output (Y tr ), the testing dataset output (Y test ), and the complexity of the NARX-NN structure can be expressed as where y tr (i) and y test (i) are the samples in Y tr and Y test , respectively. When the NARX-NN model is identified, y tr (i) are the predictions based on the training dataset, and the generalization capability is validated by the produced y test (i).
In addition, the structure complexity (m+p+r) is considered, and λ is a weighting coefficient. It is noteworthy that the lower the value of λ, the weaker the effect of structure complexity. During the model training and validation stage, it is preferable to adopt k-fold cross validation for dividing the input/output data into Y tr and Y test . This helps reduce the influence of overfitting and enhance the validity of the identified models. This identification strategy has often been used in the evaluation of intelligent optimization algorithms.
Afterwards, the crossover and mutation operators as described in [35] can be executed to generate the offspring.
To guarantee a compact structure of the NARX-NN model produced by the combinatorial GA and to improve the associated global search capability, the following operator maintenance needs to be performed. For the current population {C i } and next population C i , if the fitness function difference d J (C i ) = J (C i ) − J C i between two individuals from the two populations with the same structure is excessively small, e.g. d J (C i ) ≤ ε 1 , then, the initial values of the backpropagation error algorithm for the offspring with the same structure is selected randomly. Otherwise, the associated initialization is inherited from the individual of the former two populations with smaller object function value.
In a real levitation condition awareness application, it appears more efficient to deploy the model training and validation sub-stages separately. This means that the training sub-stage can be performed just using typical test dataset in an off-line centralized manner, while the validation sub-stage can be carried out online in each decentralized LCM with the trained model providing the levitation gap predictions.
Afterwards, the levitation gap residuals from these decentralized LCMs can be uploaded periodically to the centralized levitation condition awareness system for subsequent data-driven random matrix analysis.

B. DATA-DRIVEN RANDOM MATRIX ANALYSIS
The basic data-driven random matrix analysis can be summarized in the following three stages: streaming data based random matrix formulation, recursive eigenpairs updating of streaming covariance matrix (SCM), and asymptotic characteristics analysis of the SCMs.

1) STREAMING DATA MODELING
For maglev train typically configured with three carriages, 20 levitation gap residual signals for each carriage uploaded from decentralized LCMs are readily accessible. In general, the sampling process is synchronized. Fig. 6 shows the conceptual representation of the streaming data structure Specifically, at the kth time-stamp, the formulated matrix is X k ∈ R m×n t , formed by the data x k 1 , x k 2 , . . . , x k n t , which are modeled as vectors x i ∈ R. As the size of X k (rows m for number of dimensions, and columns n t for number of samples) increases, so do the complexity and the relationship underneath the data, and concentration of the spectral measure will occur according to the random matrix theory [21].
Since the sampling matrix X k may be non-Gaussian standardized, a data preprocessing procedure needs to be adopted to obtain the normalized non-Hermitian matrixX k ∈ R m×n t . For i th row, 1 ≤ i ≤ m, wherex i and σ x i are the mean and standard variance ofx i , in whichx i = x i + ε i ; and ε i ∈ R n t is an additive white noise vector with a sufficiently large noise-to-signal ratio (NSR). Then, the normalized SCMM k ∈ R m×n t is constructed in the form below

2) RECURSIVE EIGENPAIRS UPDATING
For SCM analysis based systematic condition awareness, it is time-consuming to directly apply independent eigendecomposition to the formulated SCM in each sliding window. Since most SCMs entries in the two adjacent sliding windows are exactly identical, the relationship between two adjacent SCMs can be utilized to derive a more efficient algorithm based on rank-one modification for recursive eigenvalue updating, thereby reducing the computational complexity.
In the process of recursive eigenvalue updating, a firstin first-out (FIFO) sample queue is maintained with a fixed number of samples in the w-width sliding window. For the k th and (k + 1) th normalized sample matrices in the two adjacent sliding windows The associated normalized SCMsM k andM k+1 are represented by where D k and D k+1 are eigenvalue diagonal matrices; V k and V k+1 are corresponding left eigenvector matrices, respectively. Hence, (9) whereỹ k+1 = V T kx k+1 . Using rank-one modification, Then, substituting (13) into (12) yields Nevertheless, by rank-one modification, we can obtain and eventually, we have Due to the special diagonal-plus-rank-one structures of (10) and (12), the associated eigenvalues {λ i } are the zeros of the secular function [36], [37] for a generalized eigenvalue problem of a normalized real symmetric matrix in the form of D+zz T , In the LAPACK subroutines DLAEDx, the algorithms for solving (14) with corresponding eigenvectors have been implemented, which can be easily called by real-time programs such as Intel MKL (called by sstedc function) [38].

3) ASYMPTOTIC CHARACTERISTICS ANALYSIS
Since efficient eigenpairs updating of formulated slidingwindowed SCM has been achieved with the recursive algorithms, further RMT-based analysis can be performed statistically with the following meaningful asymptotic spectral laws.
On the one hand, for the normalized non-Hermitian random matrix X k (the associated SCM is denoted as X = XX T n t ∈ R m×m ), if the entries satisfy the i.i.d condition, then the eigenvalue distribution (empirical spectral distribution, ESD) of the singular value equivalentX u ∈ C m×m almost surely converges to the limit given by Ring Law [16], [17] as m, n t → ∞, and c = m n t ∈ (0, 1], whereX u = U n t √ m; and U ∈ C m×m is a standardized Haar unitary matrix. Benefiting from the Haar orthogonal VOLUME 8, 2020 mapping of the SCM eigenvalues from the real to the complex plane, the circular bounds of the ESD in the complex plane in i.i.d condition are also naturally given. This offers a more intuitive visualization approach for the systematic condition assessment based on the eigenvalue distribution of streaming SCMs. Based on the Ring Law, an averaged linear statistics of eigenvalues (ALES) based on information entropy (IE) is proposed to indicate the levitation condition of the maglev train, (16) where λ iX u is the eigenvalue set ofX u . The IE-based nonlinear test function is defined as Combined with the lower bound given by (15), the control limit (CL) of the Ring Law-based ALES (R-ALES) can be derived as, On the other hand, the ESD of X almost surely converges to the limit given by the Marchenko-Pastur Law (M-P law) as m, n t → ∞, and c = m n t ∈ (0, 1], On this basis, [39] proved that N • m [φ] is defined as Converging in distribution to a Gaussian random variable with zero mean and the variance where {λ i } is the eigenvalue set of X ; , and a m = c + 1. Based on (24), the theoretical bounds of the M-P law-based linear statistics of eigenvalues (MP-LES) indicator can be derived as Therefore, a combination of (16), (17), (18) and (24) can be utilized to provide the warning/alarm bounds of the levitation condition for the maglev train.

IV. RESULTS AND DISCUSSION
In the following two subsections, the implementation of the NARX based levitation control model and the levitation condition awareness evaluation based on the proposed architecture are illustrated, respectively.

A. LEVITATION CONTROL MODEL IMPLEMENTATION
As the control target of each LCM is to keep the associated levitation gap nearly constant, the resulting levitation gap series in the case of No. 5 LCM has the form shown in Fig. 8, which appears stationary. The ACF and PACF of the levitation gap series are presented in Fig. 7, while the cross-correlation between the levitation gap series and the current and vibration-acceleration series are depicted in Fig. 8.
From the correlograms shown in Figs. 7 and 8, it is obvious that the first two delayed samples (at least) of both the levitation gap series and the current/vibration-acceleration series have some influence on the present levitation gap value, since the line segments corresponding to these delays are above the predefined confidence level indicated by the blue line. Accordingly, the minimum TDL is identified as 2 for both the output and input signals of the LCM model. Furthermore, the proposed combinatorial GA algorithm is implemented to simultaneously optimize the structure and weight-bias parameters of the NARX-NN-based LCM model. The population size N is set as 50, and the SCG backpropagation epoch is set as 20 to avoid over-fitting. The probabilities of crossover and mutation operators are set to 0.9 and 0.1, respectively. λ is set to 0.004, and ε 1 = 0.05. In addition, the structure parameters are predefined based on a little expert  knowledge. That is, the upper bounds of {m, p, r} are set as {N m , N p , N r }={9, 9, 30}. Modeling identification is carried out based on the test dataset under normal operating condition, as shown in Fig. 10, and 5-fold cross-validation is utilized.
In Fig. 9, the convergence curve of the proposed combinatorial GA algorithm for No. 5 LCM model identification is presented. To facilitate 3D display of the convergence process, the structure parameters {m, p, r} are divided into two parts, r is separately labeled as ''Number of Neurons'', and the weighted TDL is calculated as N m · m + p. The proposed algorithm is proven to be fast-convergent, and the optimal object function value can be derived within dozens of iterations. The optimal NARX-NN-based LCM model can be obtained, where {m, p, r}={5, 6, 6}.
After that, to show the efficiency of the proposed GA based NARX-NN modeling methodology (MGAC-NARX) for LCM characterization, some other approaches, i.e., the ANFIS using FCM clustering (ANFIS-FCM), adaptive GA based fuzzy neural network (FNN-APTGA) and two traditional training algorithms based NARX-NN models (BR-NARX and LM-NARX) are selected for comparison in terms of modeling accuracy and experiment running time. It should be noted that, the model structure parameters cannot be optimized simultaneously in the two traditional training algorithms, and the structure parameters are the same as those obtained by the MGAC-NARX methodology. It can be observed from TABLE 1 that, the proposed method is superior to the other methods in terms of the number of mean rules/neurons and the average modeling RMSE, especially when compared to the two traditional algorithms using NARX-NN structure, although the identified delays for input and output signals are slight larger than the ANFIS-FCM and FNN-APTGA methods. It is concluded that the modified combinatorial GA algorithm is efficient to obtain good prediction accuracy, while retaining the compact structure which is benefited from the operator maintenance introduced in this work. However, the running time of the proposed method is longer than most of the other methods because the number of training data increases fivefold and more evolutionary generation is required during the optimization process. Nevertheless, as described in section III. A, the training procedure can be configured offline to achieve the model parameter identification, which is the most time-consuming step. Therefore, the longer running time has little influence on the efficiency of the online LCM model residual generation.
In Fig. 10, the levitation gap predictions for the No. 5 LCM are displayed and are based on the identified model structure. Table 2 lists the levitation gap residual results based on the samples provided by 20 levitation gap sensors located at the middle carriage M.
Benefited from the fast dynamics of the decentralized LCMs, the local levitation gaps can be regulated to fluctuate around 8.5 mm as expected. However, due to the comprehensive influences of the dynamic track condition, train-track coupling vibration and power receiving condition on the normal operation of maglev train, the levitation gaps generally fluctuate within a wide range. According to Figs. 7 and 8, the levitation gap dynamics are obviously temporally correlated. With the aid of the NARX-NN-based LCM modeling, the temporal correlation of the levitation gap during a normal train operating condition can be evaluated effectively. In addition, spatial correlation among adjacent LCMs is fully considered by introducing the vibration-acceleration time-   series. The RMSEs of levitation gap residuals are shown in Fig. 10 (for No. 5 LCM) and Table 1 (for the 20 LCMs in carriage M). Despite the insufficiently small number of predicted residuals of the levitation gaps, the statistical characteristics of the residuals tend to be normal distribution compared to that of the original levitation gap dataset, as indicated in Fig. 11. Further verification, the inherent correlation inside LCMs under a normal operating condition can be identified effectively with the proposed decentralized LCM model.

B. LEVITATION CONDITION AWARENESS EVALUATION
Based on the implementation of the NARX-NN-based levitation control modeling, as illustrated in subsection IV.A, the levitation condition awareness performance is evaluated in MATLAB R /SIMULINK R . As emphasized in subsection III. A, the training of an individual LCM model can be performed in an offline centralized manner while the validation can be carried out online with the trained model providing the levitation gap predictions. Therefore, the levitation gap residual generation from LCM models using the MGAC-NARX method are implemented using C programming in CMEX S-function format which provides a powerful mechanism for extending the capabilities of the Simulink R environment. On this basis, streaming random matrices can be provided to simulate the sampled dataflow uploaded from individual LCMs. Subsequently, the RMT-based systematic levitation condition awareness is also implemented in CMEX S-function format integrated with encapsulated Fortran subroutines. This implementation can take full advantage of the performance of Simulink engine. In addition, the verified architecture codes will be able to be complied and transplanted more conveniently in future studies.
Based on the above implementation, the curves of Ringlaw-based ALES (denoted as R-ALES) and M-P law-based LES (denoted as MP-LES) indicators obtained via the mathematical procedure of the data-driven random matrix analysis are illustrated in Fig. 12, as well as the associated control limits and the obtained abnormal dataset.
As shown in Fig. 12, in the way-in breaking and wayout acceleration sections (e.g., [150 s, 500 s] and [800 s, 1000 s] intervals), the velocity slope of the maglev train is steep. During these periods, the rapid change in the train tractive driving state is accompanied by the train-track vibration, which imposes a coupling effect on the levitation control stability. As a result, the control limits of the proposed R-ALES and MP-LES indicators are both exceeded in the short term (during the [800s, 1000s] interval). Nevertheless, due to the LCM configuration that has sufficient stability margins, the resulting disturbance can be quickly eliminated with no abrupt alteration of associated levitation control variables, and the levitation control condition will return to steady state in a gradual manner.
However, if the coupling effect of the train tractive driving system is greater than the equivalent stability margin of the configured LCM controllers, the synergy of the multi-node  LCM control process can be deteriorated, which can further result in operation anomaly of the LCS (during the [150s, 500s] interval). Based on the proposed levitation condition awareness architecture, a retrospective analysis of a typical anomalous process is conducted, as shown in the enlarged snapshot of Fig. 13. Stage I of the process is defined as incipient abnormality stage. In this stage, the occurrence of incipient fault in the acceleration of the No. 5 LCM resulting from the aforementioned coupling effect can be captured sensitively by the proposed condition awareness indicators. If the incipient fault cannot be perceived effectively, it propagates in the LCS through cooperative closed-loop control, as shown in Stage II defined as abnormality development stage. In this stage, the acceleration fluctuation of the No. 9 LCM is intensified abnormally, and the ''abnormal energy'' quantified by the proposed indicators further increases. Finally, in Stage III defined as the abnormality deterioration stage, the acceleration fluctuations of the No. 17 and No. 19 LCMs are also abnormally enlarged, with the range expanding from [−18, 0] to [−50, 30]. Furthermore, the ''abnormal energy'' increases sharply, thus seriously affecting the operational stability of the LCS.
Based on the above analysis, the qualitative relationships between the proposed indicators and the abnormal levels are summarized in Table 3. The ''abnormal energy'' indicated by the proposed statistics shows significant positive correlation with the severity of the abnormity. This also verifies the rationality of the proposed condition awareness indicators. In addition, LES seems to be more sensitive than ALES for abnormity detection.
In contrast with the above comprehensive results for the abnormal operation of the maglev train, curves of the LCS condition awareness under normal operating condition are illustrated in Fig. 14. Neither R-ALES nor MP-LES significantly exceeds respective control limits. In addition, the condition awareness indicators show similar increases during the train acceleration and braking processes. This also suggests that the tractive driving condition produces coupling effect on the LCS of the maglev train.
Therefore, the designed R-ALES and MP-LES indicators are proven to be effective in the timely recognition abnormal LCS operation conditions under the proposed condition awareness architecture, and can provide pre-alarm indication in order to prevent the levitation control system from encountering severe faults.

V. CONCLUSION
In this paper, a levitation condition awareness architecture for low-speed maglev train is proposed based on data-driven random matrix analysis. The proposed architecture is comprised of two parallel procedures: an engineering procedure based on decentralized data-driven NARX-NN modeling, and a mathematical procedure of centralized data-driven random matrix analysis. During the decentralized engineering procedure, the NARX-NN modeling for individual LCMs can be achieved by associated parameter identification using the modified combinatorial GA algorithm, and this is the data foundation for the subsequent mathematical procedure. The centralized mathematical procedure is the core algorithm of the proposed architecture, which is enlightened by the random matrix theory. Combined with recursive eigenvalue updating of SCM, the associated eigenvalues can be calculated efficiently to satisfy the real-time application requirement. On this basis, the high-dimensional analysis for the SCM and the associated singular equivalent of the LCS for the maglev train can be performed, and the control limits of the constructed R-ALES and MP-LES indicators can be well-designed based on the theoretical distributions. The results demonstrate that the time-varying behavior of individual LCMs can be well approximated by a model with compact NARX-NN structures. The correlation among the input-output signals can be extracted effectively while retaining the random characteristics of the associated prediction residuals. Furthermore, systematic condition awareness can be quantified based on the proposed high-dimensional statistics (R-ALES and MP-LES indicators) and associated control limits. Through comprehensive analysis based on the experimental datasets obtained under normal/abnormal train operating conditions, the rationality of proposed levitation condition awareness architecture for low-speed maglev train is verified. This proposed awareness architecture can detect the incipient train-guideway interaction vibration abnormality in a timely, and provides a framework and feasible research methods for the real-time diagnosis of the incipient/small state abnormalities of the LCS.