Simulation and Dynamic Properties Analysis of the Anaerobic–Anoxic–Oxic Process in a Wastewater Treatment PLANT Based on Koopman Operator and Deep Learning

: The accurate simulation of the dynamics of the anaerobic–anoxic–oxic (A2O) process in the biochemical reactions in wastewater treatment plants (WWTPs) is important for system prediction and optimization. Previous studies have used real-time monitoring data of WWTPs to develop data-driven predictive models, but these models cannot be used to provide mathematical analysis of A2O dynamic properties. In this study, we developed a new simulation and analysis method for determining A2O dynamics in biochemical reactions using deep learning and the Koopman operator to address the above problems. This method was validated through data from a real-world WWTP in east China and compared it with the traditional deep learning model. According to the results, the new method achieved high-accuracy prediction. Meanwhile, with the help of the Koopman operator, the new method was able to analyze the asymptotical stability and convergence behavior of the A2O process, which provides a brand-new perspective for the in-depth study of biochemical reactor dynamics.


Introduction
The anaerobic-anoxic-oxic (A2O) process in biochemical reactors of wastewater treatment (WWT) is important in urban water systems [1].Mechanism models of A2O have been developed for dynamic prediction and control [2][3][4][5].However, they are built on fixed parameters, equations, and model assumptions, which lack the flexibility and adaptability to respond to external environmental changes and stochastic perturbation in realworld conditions.With the widespread adoption of intelligent devices, data monitoring systems, and SCADA, a large amount of dynamic data for the A2O process has gradually accumulated.As a result, the construction of A2O dynamic models based on datadriven methods has gradually emerged, e.g., linear models [6], machine learning [7], artificial neural networks [8,9], and deep learning models [10].These methods establish simulation and prediction models based on measured data that records various random events and perturbation factors.Therefore, they are more generalizable and adaptable to real-world situations.
However, existing data-driven models, especially deep learning models, are black-box models that lack interpretability and cannot be used to analyze the dynamic properties of the A2O process in a biochemical reactor.This is consistent with existing research on the dynamic properties based on deep learning and some data-driven methods [11][12][13][14].Therefore, finding the intersection point between data-driven methods and dynamic analysis through appropriate modeling for A2O dynamics is necessary.
The Koopman operator provides a solution to the above problem.For a given nonlinear dynamic system, obtaining its Koopman operator can help us find its eigenvalues, thereby revealing its linear structure and analyzing its dynamic properties mathematically [15][16][17][18][19][20][21].The theory and correlated methods have been widely used in model prediction [22], data fusion [23], fluid dynamics [24,25], and system control [26][27][28].The methods used to find the Koopman operator of a given system include purely mathematical methods [29] and data-driven methods [30,31].In recent years, deep learning has also been used to find the Koopman operator of a given nonlinear dynamic system, providing a new direction for mathematical analysis and system simulation [32][33][34].
This study constructed a new data-driven simulation model of A2O process in WWTP to solve the above problem.This model provides accuracy prediction and mathematical analysis for A2O dynamics using Koopman operator and deep learning.Based on this, a deep learning model that takes into account both dynamic properties analysis and accurate predictions is built for A2O to provide new insight into the dynamics of A2O from a mathematical perspective.

Materials and Methods
Deep learning and the Koopman operator were applied to establish a new simulation method and a dynamic analysis method for the A2O process in this study.First, the data for an A2O process in a WWTP in eastern China were collected and cleaned.Based on the cleaned dataset, a traditional deep-learning method was used to construct a simulation model as a baseline model for A2O dynamic prediction.Then, a new simulation model based on the combination of deep learning and the Koopman operator was designed and used to provide a mathematical analysis method for A2O dynamic properties.The following sections provide details on the methods.

Case Study and Data
The case study is a real-world WWTP located in eastern China.The WWTP adopts the A2O process, with a design flow rate of 800,000 m 3 /d and 8 groups of bioreactors, each with a capacity of 100,000 m 3 /d and can operate independently.The A2O process uses microorganisms to remove organic pollutants, ammonia, and phosphorus in water through anaerobic, anoxic, and aerobic treatment processes, so as to meet the discharge standards [35].The aerobic section of the biochemical reactor uses blowers at the bottom of the reactor to provide oxygen for microorganisms to grow, consume organic pollutants, and prevent activated sludge from settling.This study only focuses on bioreactor no. 1, which has an online data monitoring system.It has three parts, anaerobic, anoxic, and oxic; their volumes are 55,696 m 3 , 133,670 m 3 , and 321,680 m 3 , respectively.
This study collected the data for the bioreactor no.1 in the WWTP from 1 May 2020, at 00:00 to 21 October 2020, at 14:40.The monitoring sample frequency was 5 min.Four types of data (state, control, inflow, and outflow) were collected.Their information is given in Table 1.The state data represent the water quality data inside the bioreactor, which included dissolved oxygen (DO) in the aerobic section, oxidation-reduction potential (ORP) of the anaerobic and anoxic sections, nitrate nitrogen (NO − 3 ) of the aerobic section, and mixed liquor suspended solids concentration (MLSS) of the anaerobic, anoxic, and aerobic sections.The control data represent the data related to aeration (aeration volume), internal recirculation flow rate (Qr), and sludge internal recycle flow rate (Qsr).The inflow data represent the data of inflow water, which includes chemical oxygen demand (COD), total nitrogen (TN), total phosphorus (TP), water temperature (T), suspended substance (SS), and flow rate (flow).The outflow data represent the data of outflow water, which includes chemical oxygen demand (COD), total nitrogen (TN), total phosphorus (TP), ammonia nitrogen (NH + 4 − N), nitrate nitrogen (NO − 3 ).The original COD, TP, and TN data were measured every two hours and their time series data comprised the same values every two hours, resulting in the original data being shown in a stepwise manner.Due to the presence of outliers in the original data, data cleaning was performed.Based on the PauTa Criterion, extreme outlier data were removed and data were supplemented using interpolation.Firstly, the time interval for data cleaning was determined based on the hydraulic retention time (HRT) of the bio-reactor, and the mean and variance of various state data of the bio-reactor within the time interval were calculated.Then, according to the hypothesis testing PauTa Criterion, extreme outlier data within the time interval was deleted and supplemented using interpolation data.The criterion can be described by Equation (1), and the interpolation is shown in Equation ( 2), where N represents the time steps in the time interval, data i represents the i-th data point in the time interval, and data i represents the supplementing data.
After deleting the abnormal data, the moving average method was used to eliminate the small disturbances caused by the anomalies (Equation (3), where data t represented the data that have been corrected, n was the time span of the moving average and was 100 in this study).After being processed, the data relieved the impact of abnormal and disruptive values, providing a foundation for subsequent analysis and modeling.In addition, all the data were smoothed through this process.The specific details and statistical results of the cleansed dataset are shown in the table below.

Deep Learning for A2O Process Simulation and Prediction
The above data contain information on the influent, anaerobic, anoxic, and aerobic reaction processes, and effluent quality, covering the elements of carbon, nitrogen, and phosphorus, all of which have an impact on the dynamic of A2O process.Therefore, this study used the above data to construct a deep learning-based simulation model that focused on the inflow-state-control-outflow of A2O process.The inflow, state data, and control data were used as inputs for the simulation model, and the outflow data were used as the outputs of the simulation model.All input data were time series data for a period of time with a time span of hydraulic retention time (HRT), which was 7 h in this study.The model predicted the outflow data for the next half hour.
Accordingly, the deep learning-based model (called DNN in this study) can be described as Equation ( 4), where State t , Control t , In f low t are the state, control, and inflow in previous HRT, Out f low t+1 is the outflow in next half hour and is the output of the model.The training process of DNN is to minimize the loss function based on prediction error Equation ( 5), where Out f low t+1 is the output of model.All the data for inputs and outputs were normalized to avoid the influence of differences in data magnitude, and the results of the simulation were obtained by renormalization.The formula is shown in Equations ( 6) and (7), where DATA, DATAmin, and DATAmax are the original data and its minimum and maximum, data is the normalized data, and DATA is the renormalized data.

Out f low t+1
DNN is trained using the standard Adam training algorithm [36], with 0.001 initial learning rate and 500 iterations.After training, the model performs prediction by inputting the corresponding state, control, and inflow variables.The structure and parameters of DNN are shown in Figure 1 and Table 2, respectively.

Koopman Operator and Deep Learning for A2O Process Simulation and Prediction
Due to the lack of interpretability in deep learning methods, it is difficult for DNN to perform in-depth analysis of the A2O process dynamics.In this study, a new A2O dynamic simulation model was constructed by combining Koopman operators and deep learning methods.

Dynamic System and the Koopman Operator
The definition of the Koopman operator and its approximation theory can be found in previous works [32,33].We only provide a step-by-step introduction in this section.The dynamic of A2O can be described by a nonlinear dynamic system Equation ( 8), which has a corresponding Koopman operator to describe the evolution of a given observable through a linear system (Equation ( 9)), where  is Koopman operator and  is observable represented by functions in a Hilbert space [31].

Dynamic System and the Koopman Operator
The definition of the Koopman operator and its approximation theory can be found in previous works [32,33].We only provide a step-by-step introduction in this section.The dynamic of A2O can be described by a nonlinear dynamic system Equation ( 8), which has a corresponding Koopman operator to describe the evolution of a given observable through a linear system (Equation ( 9)), where κ is Koopman operator and f is observable represented by functions in a Hilbert space [31].
However, this linear system has an infinite dimension [23,30] and cannot be used for simulation and prediction directly.Thus, its finite dimensional approximation is required.Given a set of dictionary functions Ψ = [ψ 1 , ψ 2 , . . . ,ψ L ], the observable can be written as f (x(n)) = ΞΨ(x(n)), where Ξ = [ξ 1 , ξ 2 , . . . ,ξ L ] T is corresponding weights.Then Equation (10) can be obtained, where K, called Koopman matrix, is the representation of the Koopman operator κ and r is the residual term.By minimizing r, we obtain the corresponding Ψ and KΞ to formulate a finite dimensional linear system [37].This system (Equation ( 11)) is a finite dimensional linear system that can be used for simulation and prediction.

Dictionary Learning-Based Extended Dynamic Mode Decomposition
An approximation method is required to obtain Equation (11).Regression-based approximation methods use the regression model [29,38], or even deep learning model [33,34] to find a linear system, and is mainly considered in this study.One of the methods based on deep learning, called dictionary learning extended dynamic mode decomposition (DLEDMD) [31,33] is used as the approximation algorithm in this study to identify the dictionary functions and Koopman matrix.First, N pairs of data  12), where w 1 is the trainable parameter, K is the Koopman matrix, and λ(K, w 1 ) is a suitable regularizer.The trained K, φ are used to obtain a model for simulation and prediction.This trained model is called DLEDMD in the rest of this paper.Its architecture and parameters are given in Figure 2 and Table 3. Considering the stable and convergence of the training process for general initializations, the emulator is trained through the iterating method [33] combined with the Adam algorithm [36,39] with 0.001 initial learning rate and 500 iterations.The K is computed by Equation ( 13) during applications.
Water 2023, 15, x FOR PEER REVIEW 7 of 12 Then, it goes through three separate deep neural networks and a combined layer for embedding and combination.This part is the  in Equations ( 12)-( 14).All the inputs and their combination are used to obtain the Koopman matrix.This is the only difference between DNN and DLEDMD.The final output or result  +1 is calculated by the final deep neural network for simulation and prediction purposes.Then, it goes through three separate deep neural networks and a combined layer for embedding and combination.This part is the φ in Equations ( 12)-( 14).All the inputs and their combination are used to obtain the Koopman matrix.This is the only difference between DNN and DLEDMD.The final output or result Out f low t+1 is calculated by the final deep neural network for simulation and prediction purposes.

Dynamic Properties Analysis Based on Koopman Operator
The linear structure of DLEDMD can be easily obtained through Equation ( 14), where W is the matrix of eigenvectors, Σ is the diagonal matrix containing the eigenvalues, and ψ(X, a, rain; w 1 ), i = 1, 2, . . ., L are the eigenfunctions of the Koopman emulator.
If the eigenvalues and eigenvectors of the Koopman matrix cannot be obtained, the linear structure can also be given through the singular values and singular vectors of the matrix K [40].Through the eigenvalues, the dynamic properties of A2O process can be directly analyzed and evaluated using the theorems of existing dynamic system and system, including the Lyapunov stability and convergence of the dynamics [16].This can determine the ability of a biochemical reactor to withstand external disturbances and its long-term behavioral characteristics during operation.

Results and Discussion
DNN and DLEDMD were trained using the method provided in Sections 2.2 and 2.3.Then, all the trained models were used to predict and simulate the dynamic of A2O in the case study area.After that, an analysis of A2O's dynamic properties was provided based on the Koopman matrix of DLEDMD.The results and discussion are given as follows.

Training Process of DNN and DLEDMD
The loss function curves of DNN and DLEDMD during training process are shown in Figure 3.The training data consisted of randomly selected 1000 input-output pairs in the given time spans, and the validation data consisted of randomly selected 500 input-output pairs in the given time span.Each pair consisted of State t , Control t , In f low t , and Out f low t+1 .The training method and parameters were those provided in Sections 2.2 and 2.3.From the training results, it is seen that both models progressively reduce their errors on the training set, demonstrating good training performance.They also show improvement on the testing set during the training process, indicating an increase in their generalization ability through training.However, DNN exhibits faster and better improvement on the testing set, suggesting that it has relatively better generalization ability compared to DLEDMD.

Simulation Performance
All the trained models were used to predict and simulate the outflow dynamic of A2O in the case study area.The predicted time ranged from 1 May 2020 at 0:00 to 21 October 2020 at 14:40.The models predicted the outflow in half an hour every 5 min, i.e., the time interval of prediction was half hourly and the frequency of prediction was 5 min, and there were a total of 49,901 predictions.The mean square error (MSE, Equation (15), where Y t is the output of models, Y t+1 is the true output values, Num is the number of prediction, and is 49,901 in this study) and Nash-Sutcliffe efficiency coefficient (NSE, Equation (16), where E Y t+1 is the mean value of true output values) were used to measure the performance of the two models.
The comparison between the predicted results of the two models and the actual values, as well as the MSE and NSE results, are shown in Figure 4.It can be seen from the figure that both data-driven models capture minor system disturbances and achieve good predictions with NSE reaching 0.99.At the same time, the comparison shows that the prediction performance of DLEDMD is slightly worse than that of DNN, which is consistent with the generalization performance in the training process of the two models.

Results and Discussion
DNN and DLEDMD were trained using the method provided in Sections 2.2 and 2.3 Then, all the trained models were used to predict and simulate the dynamic of A2O in the case study area.After that, an analysis of A2O's dynamic properties was provided based on the Koopman matrix of DLEDMD.The results and discussion are given as follows.

Simulation Performance
All the trained models were used to predict and simulate the outflow dynamic o A2O in the case study area.The predicted time ranged from 1 May 2020 at 0:00 to 21 Oc tober 2020 at 14:40.The models predicted the outflow in half an hour every 5 min, i.e., the time interval of prediction was half hourly and the frequency of prediction was 5 min, and there were a total of 49,901 predictions.The mean square error (MSE, Equation (15), where   is the output of models,  +1 ̅̅̅̅̅ is the true output values,  is the number of predic tion, and is 49,901 in this study) and Nash-Sutcliffe efficiency coefficient (NSE, Equation (16), where [ +1 ̅̅̅̅̅ ] is the mean value of true output values) were used to measure the performance of the two models.

Eigenvalues of Koopman Matrix and Asymptotically Stablility of A2O Dynamic
Based on DLEDMD, we can directly obtain its corresponding Koopman matrix and use its eigenvalues to obtain the linearized description of the A2O nonlinear dynamical system through Equation ( 14).The eigenvalues of DLEDMD's Koopman operator are shown in Figure 5. From the figure, it can be seen that the eigenvalues of the simulated dynamic of A2O process in this case study are all less than 1.
According to the Lyapunov stability theorem, a discrete time system such as DLEDMD in Equation ( 11) is asymptotically stable if and only if all its eigenvalues are inside the unit circle [16].Therefore, the dynamic of the A2O in this case is asymptotically stable, and its output variables will gradually converge as the system runs without external influences.This means that the biochemical reactor in this case is capable of resisting fluctuations in influent water quality and control variables within a certain range.That is, the water quality of the effluent of the system will converge to a certain value and maintain stability when there is no continuous change in the influent water quality.Such a mathematical analysis can never be provided by DNN because the eigenvalues cannot be obtained from it.
The comparison between the predicted results of the two models and the actual values, as well as the MSE and NSE results, are shown in Figure 4.It can be seen from the figure that both data-driven models capture minor system disturbances and achieve good predictions with NSE reaching 0.99.At the same time, the comparison shows that the prediction performance of DLEDMD is slightly worse than that of DNN, which is consistent with the generalization performance in the training process of the two models.

Eigenvalues of Koopman Matrix and Asymptotically Stablility of A2O Dynamic
Based on DLEDMD, we can directly obtain its corresponding Koopman matrix and use its eigenvalues to obtain the linearized description of the A2O nonlinear dynamical system through Equation ( 14).The eigenvalues of DLEDMD's Koopman operator are shown in Figure 5. From the figure, it can be seen that the eigenvalues of the simulated dynamic of A2O process in this case study are all less than 1.According to the Lyapunov stability theorem, a discrete time system such as DLEDMD in Equation ( 11) is asymptotically stable if and only if all its eigenvalues are inside the unit circle [16].Therefore, the dynamic of the A2O in this case is asymptotically stable, and its output variables will gradually converge as the system runs without external influences.This means that the biochemical reactor in this case is capable of resisting

Conclusions
Although data-driven methods have shown promise in the prediction of A2O dynamic in WWTPs, they cannot be helpful in the analysis of A2O dynamic properties.This study designed a data-driven prediction model using Koopman operators and deep learning methods to solve this problem.The method was validated using data from a real WWTP and compared to traditional deep learning methods.On the one hand, with the help of abundant data, this model accurately predicted the dynamics of the A2O process and captures small disturbances in practical situations.On the other hand, the Koopman operator was used to provide an analysis method for the dynamic properties of the A2O process from the perspective of mathematics.
Both DNN and DLEDMD learned the dynamics of the data through training and made relatively accurate predictions with NSE exceeding 0.99.However, according to the testing set error during the training process and the MSE and NSE values in the prediction, DNN had better generalization ability than DLEDMD because of its higher non-linearity.DLEDMD can be used to analyze the dynamic properties of the A2O process by utilizing its Koopman operator.In the case of this study, the trained DLEDMD directly obtained the Koopman matrix and its corresponding eigenvalues.The obtained eigenvalues indicates that nonlinear dynamic system of A2O process in this study was Lyapunov asymptotic stable.This proves that the state and the effluent of the system will not undergo significant changes within a certain range of influent water quality disturbance in this case study.Moreover, when there is no sustained change in the influent water quality, the effluent quality of the A2O process will converge to a certain value and remain stable.

12 Figure 1 .
Figure 1.Architecture of DNN.The input of the DNN consists of   ,   ,   .Then, it goes through three separate deep neural networks and a combined layer for embedding and combination.The final output or result  +1 is calculated by the final deep neural network for simulation and prediction purposes.

Figure 1 .
Figure 1.Architecture of DNN.The input of the DNN consists of State t , Control t , In f low t .Then, it goes through three separate deep neural networks and a combined layer for embedding and combination.The final output or result Out f low t+1 is calculated by the final deep neural network for simulation and prediction purposes.

Figure 2 .
Figure 2. Architecture of DLEDMD.The input of the DLEDMD consists of   ,   ,   .Then, it goes through three separate deep neural networks and a combined layer for embedding and combination.This part is the  in Equations (12)-(14).All the inputs and their combination are used to obtain the Koopman matrix.This is the only difference between DNN and DLEDMD.The final output or result  +1 is calculated by the final deep neural network for simulation and prediction purposes.

Figure 2 .
Figure 2. Architecture of DLEDMD.The input of the DLEDMD consists of State t , Control t , In f low t .Then, it goes through three separate deep neural networks and a combined layer for embedding and combination.This part is the φ in Equations (12)-(14).All the inputs and their combination are used to obtain the Koopman matrix.This is the only difference between DNN and DLEDMD.The final output or result Out f low t+1 is calculated by the final deep neural network for simulation and prediction purposes.
The loss function curves of DNN and DLEDMD during training process are shown in Figure3.The training data consisted of randomly selected 1000 input-output pairs in the given time spans, and the validation data consisted of randomly selected 500 inputoutput pairs in the given time span.Each pair consisted of   ,   ,   , and  +1 .The training method and parameters were those provided in Sections 2.2 and 2.3.From the training results, it is seen that both models progressively reduce their errors on the training set, demonstrating good training performance.They also show improve ment on the testing set during the training process, indicating an increase in their gener alization ability through training.However, DNN exhibits faster and better improvemen on the testing set, suggesting that it has relatively better generalization ability compared to DLEDMD.

Figure 5 .
Figure 5. Eigenvalues of DLEDMD.Each blue point represents an eigenvalue, which is a complex number with both imaginary and real parts.

Figure 5 .
Figure 5. Eigenvalues of DLEDMD.Each blue point represents an eigenvalue, which is a complex number with both imaginary and real parts.

Table 1 .
Statistic information of A2O process data.

Table 2 .
Architecture and parameters of deep neural network in DNN.

Table 2 .
Architecture and parameters of deep neural network in DNN.
[State t 1 , State t 2 , . . ., State t N ], [Control t 1 , Control t 2 , . . ., Control t N ], [In f low t 1 , In f low t 2 , . . ., In f low t N ], and [Out f low t 1 , Out f low t 2 , . . ., Out f low t N ] are collected.Then, a deep neural network, φ, is defined to represent the dictionary functions.It is trained through Equation (

Table 3 .
Architecture and parameters of deep neural network in DLEDMD.

Table 3 .
Architecture and parameters of deep neural network in DLEDMD.