The use of suspended microorganisms to remove undesired components, including organic carbon, nitrogen, and phosphorus species, is referred to as activated sludge processes and is widely used around the world in municipal wastewater treatment plants (WWTPs)    . Mathematical models of activated sludge processes, referred to as Activate Sludge Models (ASMs), provide a cost-effective way to evaluate design, control, and optimization of the processes and have been used extensively in practical applications and academic research     . Henze et al.  proposed the first ASM (ASM1), and since then, several other versions of ASMs have been proposed  . ASM1 represents wastewater composition using 13 constituents, where each constituent may interact with each other through 8 reactions in which the general mass balance of each variable results in a system of ordinary differential equations (ODEs) describing the change in concentrations over time. Non-linear reaction rates based on Monod Kinetic cause a system of highly nonlinear ODEs. In some applications, the computation time of the ASM is critical due to a long simulation period, the need for online simulation  , and the use of parameter estimation algorithms that require numerous simulations. For example, Alikhani et al.  modeled a large scale WWTP with an ASM by using Markov-Chain Monte Carlo (MCMC) sampling and solved a large system of ODEs for 500,000 times to obtain the ASM’s parameter posterior probability distribution. Therefore, applying a fast algorithm to solve ASMs that can significantly decrease the computation time of such applications would be highly desired. In addition to the non-linearity of the system of ODEs comprising ASMs, they typically cover a wide range of biochemical reaction rate scales, ranging from seconds (for example, oxygen transfer rate) to days (for example, microbial growth rates) and result in a mixed (stiff/nonstiff) system of ODEs  that generally requires small time-steps when conventional ODE solver algorithms are used. Furthermore, due to the inherent fluctuating behavior of the external forcing vectors, including variation in the time series of the influent rate and the characteristics in wastewater treatment streams, the optimal time-step size can vary greatly during the course of a simulation.
Explicit and semi-explicit ODE solvers have been shown not to perform well for systems of ODEs with high degrees of stiffness, while backward differentiation formulas (BDFs) have been shown to be much more suitable     . In addition, because the degree of stiffness during the course of the simulation can vary, using a fixed but small time-step can translate to a heavy computational burden. A dynamically adaptive time-step algorithm can automatically adjust the time-step according to the degree of stiffness, resulting in a more optimal use of computational resources. For example, Celaya et al.  proposed a method to select or reject a proposed time-step in adaptive algorithms by evaluating the local truncation error at each time-step to maintain the error below a given threshold.
The goal of this study is to propose a fast, adaptive algorithm based on BDFs that can be implemented in any ASM solver package and in other large, non-linear systems of ODEs. The proposed adaptive algorithm can be optimized for specific problems by properly adjusting the solver parameters and can be used in other numerical method applications such as advection-dispersion-diffusion equations solution  , Bayesian parameter estimation framework  , and other application of numerical methods (e.g.    ). The proposed adaptive BDF algorithm in this study contains three adjustable parameters that can be manipulated to achieve the minimum computation time for a given system of ODEs. In addition, a sensitivity analysis is conducted to identify the optimally adjustable parameters of the adaptive algorithm for a system of ODEs obtained for a sample ASM problem.
2.1. Adaptive Backward Differentiation Algorithm
The system of ODEs can be written in the following general form  :
where, t is the time, is the state variables vector, and is the initial condition. is the boundary condition or the external forcing vector given for the entire simulation period. The first and second order BDFs can be represented considering a variable time-step in a given time-step (i) as   :
where i indicates the time-step with non-constant step size, represents the discretized form of Equation (1) when and. In Equations (2) and (3), and are the first and second order BDFs considering non-constant step size at each time-step, respectively, which will be henceforth denoted as for the sake of generality. The conventional Newton-Raphson (NR) method for solving the system of non-linear to find can be formulated as:
where k is the iteration index beginning from 0, is the inverse of Jacobian matrix (Equations (5) and (6)). Due to the fact that inverting the Jacobian matrix is computationally intensive, especially as the size of the matrix grows, it is beneficial if the inverse of the Jacobian matrix can be reused and is recalculated only when necessary. Therefore, in the proposed method, in Equation (4) is a recycled Jacobian matrix evaluated at some previous time.
where n is the size of the ODEs, is the Kronecker Delta, and denote the Jacobian matrix for the first and second order, respectively, where is the Kth element of vector, and similarly, is the Ith element of vector at. In addition, to further speed up the NR convergence, the following initial guess of the state variables is used at the first iteration (initial guess) of the NR  :
The adaptive time interval scheme adjusts the size of time-steps based on monitoring the convergence of the NR method   . The time-step size keeps growing by the factor of 1 + ρ (ρ ≥ 0) until the NR’s iteration k in Equation (4) remains smaller than a threshold kmax, and this inverse of Jacobian matrix () will be used in Equation (4). The iteration continues until the number of iterations to achieve convergence in the NR method exceeds kmax. At this time, a new Jacobian matrix will be computed and/or the
time-step size will be reduced by the factor of (). Figure 1 depicts a de-
tailed flowchart of this algorithm, in which ρ, γ, and kmax are, respectively, the time-step inflation factor, the time-step depression factor, and the desired maximum iteration number of NR method that can be manipulated for each particular case of ODEs to achieve the minimum computation time.
2.2. Activated Sludge Model
The proposed algorithm was applied to solve an ASM based on a real, full-scale WWTP. Influent flow and its composition data were collected from nitrification-denitrification reactors at the Blue Plains wastewater treatment plant in Washington, DC over a period of 120 days. This system was modeled as three reactors in a series based on mass balance equation (Equation (8)) with a modified ASM1 reaction network consisting of 11 reactions and 15 constituents, as shown in Table 1  . The total volume of the reactor was 17,500 m3, divided into three tanks in a series, as depicted in Figure 2. Figure 3 shows the influent flow and its ammonia concentration with return activated sludge (RAS) and waste activated sludge (WAS) flow. In addition, an alternative, hypothetical influent was created by applying a moving average filter with a seven-day mask to reduce noise and a high-order frequency to evaluate the effect of temporal fluctuations on the choice of the solver’s adjustable parameters. The time-series vector of flows and concentrations and the temperature (partly shown in Figure 3) represent the in Equation (1). Significant fluctuation can be observed in the input data, especially in the inflow rate, mainly because the Blue Plains WWTP receives combined sewer flow, and storm events can generate unpredictable fluctuations in both rate and composition of influent into the plant.
Figure 1. Adaptive time-step BDF algorithm with adjustable parameters of ρ, γ, and kmax.
Figure 2. Activated sludge process containing three tanks in series with aeration in tank 1 and external loading rate in tank 3. The mixed liquor inside the bio-reactors contains both soluble and particulate constituents which are shown with S and X sign in Table 1, respectively.
Table 1. Process reaction rates (R) and constituents stoichiometeries (
Figure 3. Influent flow and ammonia concentration, RAS and WAS over time: (Left) Unfiltered; (Right) Filtered using moving average technique over a 7-day span.
In ASMs, the activated sludge process is commonly represented as a series of attached, mixed reactors. These reactors hold a series of processes that can occur simultaneously, resulting in changes in the concentration of constituents. Figure 2 depicts a schematic of the activated sludge process of this study in the WWTP. A transient mass balance can be mathematically performed for each state’s variables, consisting of the concentration of each model constituent, resulting in a system of ODEs, as expressed in Equation (8). Therefore, Equation (8) expresses the changes in constituents’ concentrations over time as controlled by influents, effluents, reactions, and various mass-transfer processes   :
where V (L3), C (M∙L−1), Q (L3∙T−1), and ω indicate the volume, concentration, flow, and flow-fraction factor (determining how much of the inflow and return flow enters each tank in a step-feed system), respectively, H is the Heaviside (unit step) function, R (M∙T−1) is the reaction rate (Table 1), ∅ is the stoichiometric coefficient (Table 1), F (T−1) is the mass transfer rate, (M∙L−1) is the saturated concentration, and (M∙L−1) is the external mass flow rate. And as for the indices, j indicates a constituent, m represents a tank, “ret” indicates a return flow, “in” shows conditions in influent, and r indicates the reaction. In addition, Nm is the total number of tanks, and Nr is the total number of reactions. is the flow rate from tank m' to tank m, due either to sequential stages or through feedback or bypass (a positive value means to flow into tank m). In Equation (8), the term on the left side is the rate of change in the total mass of the constituents j in tank m. The first term on the right side of Equation (8) is the mass inflow due to return flow; the second term is the mass inflow of the influent; the third term represents inflow from the preceding stage or feedback flow if it is present; the fourth term is the outflow of constituents; the fifth term is the production or disappearance of constituents due to reactions; the sixth term is the effect of rate-limited mass transfer (for example, aeration); and the last term is the direct addition of the constituents (for example, the addition of a carbon source for denitrification). To obtain solid-particle concentrations in the return flow (), a dynamic clarifier model  or a quasi-steady-state approximation (that is, performing mass balance while ignoring the solid storage changes in the clarifier) can be applied. The study described in this paper used the former approach.
3. Results and Discussion
The resulting system of ODEs, consisting of 15 state variables and 45 ODEs, was solved with a range of adjustable parameters for both unfiltered and filtered external vectors, and the running time for each set of solver’s parameters was recorded. The solver algorithm was implemented into the BioEst program, which was developed using the C++ language. The results presented in the next section were obtained by running the executive version of the C++ code on a desktop computer with a quad-core Intel® 2.67 GHz Xenon® CPU and 8 Gb RAM.
Figure 4 shows the results of the simulation for soluble chemical oxygen demand (sCOD), dissolved oxygen (DO), ammonia, and heterotroph biomass concentrations based on unfiltered influent obtained using the solver parameters ρ = 0.01, γ = 0.01, and kmax = 10.
3.1. Optimum Adjustable Parameters
To get the minimum computation time (running time) and acceptable accuracy, the user can assign the three adjustable parameters, including the inflation growth rate of the time-step (ρ), the depression rate of the time-step (γ), and the maximum NR’s iteration (kmax). A small initial time-step of h0 = 10−5 day is chosen for all the cases. First, the effect of changing ρ and γ is investigated by varying ρ between 0.0001 and 0.1 and γ between 0.0001 and 1 at the fixed kmax of 4. Figure 5 shows the results for effect of ρ and γ on running time. In both scenarios, filtered and unfiltered, the running time is very sensitive to ρ, and the minimum running time occurs at ρ = 0.01. Figure 5 also shows that the sensitivity of the algorithm to γ is lower than to ρ in both the filtered and
Figure 4. Results obtained for 120 days by using adaptive BDF method in a 3 tanks in series scheme with solver parameters of ρ = 0.01, γ = 0.01, and kmax = 10.
unfiltered scenarios. For example, at ρ = 0.01, in the case of unfiltered influent, γ = 0.01 gives a minimal running time and γ = 0.1 gives a similar running time, showing the low sensitivity of the algorithm to γ. The effect of the maximum number of iteration kmax is examined by assigning the values of 4, 10, 40, and 100. Figure 6 shows the effect of kmax on the running time for two sets of inflation-depression parameters. In both scenarios, the algorithm shows a reasonably high sensitivity to kmax, and a minimum running time occurs at kmax = 10. In conclusion, for the selected system of ODEs, the minimum computation time is achieved at ρ = 0.01, γ = 0.01, and kmax = 10.
For a given set of ρ, γ, and kmax, the run-time is always smaller for the filtered than for the unfiltered, as expected, because the algorithm needs smaller time- steps to resolve the higher disturbance in the unfiltered. Nevertheless, the difference in the running time is not significant.
Figure 5. Effect of ρ and γ in a constant kmax = 4 in run-time for solving the system of 45 ODEs during 120 days. (Left) Unfiltered (t); (Right) Filtered U(t).
Figure 6. Effect of kmax in run-time for solving the system of 45 ODEs during 120 days. (Left) Unfiltered (t); (Right) Filtered U(t).
3.2. Variable Time-Step Pattern
The adaptive time-step algorithm increases the size of the time-step in the case when the NR algorithm continues converging and decreases the size of the time-step when the NR fails to converge. This pattern can be seen in Figure 7, which shows the variation of the accepted time-steps during the course of the simulation for four sets of adjustable parameters. The results show that the pattern of variation in the time-step is substantially affected by different sets of adjustable parameters. For example, as Figure 7 shows, at kmax = 10, the time-step can reach a maximum value of 0.35 days, whereas, at kmax = 4, the maximum value is 0.14 days. In other words, a higher kmax allows larger time-steps and further reuse of the Jacobian matrix at the price of a higher number of NR iterations at each time-step. This indicates that the optimal value of kmax depends on
Figure 7. Adaptive time-step (dt) variation for 120 days of applying adaptive BDF method on unfiltered U(t) for different adjustable parameters.
the effort needed in the Jacobian inverting process relative to the computational cost of NR iterations, which are mainly matrix-vector multiplications. Therefore, the optimal kmax can also depend on the size of the Jacobian matrix (that is, the number of state variables) and how well-posed the Jacobian matrix is.
As Figure 7 shows, at any given kmax, when γ is varied from 0.1 to 0.01, the pattern of the time-step variation differs significantly and shows more fluctuation at a lower γ. This indicates that the solver parameters ρ and γ essentially control how much the time-step is increased or decreased. A smaller ρ and γ lead to more conservative increases or decreases of the time-step and a larger ρ and γ lead to less conservative increases or decreases.
3.3. Error Assessment
One way to evaluate the accuracy of the method is to compare the solution with the exact solution. However, the system of ODEs represented in Equation (1) cannot be solved using analytical mathematical methods to find the exact solution so that the
Table 2. Maximum and the average error of the ABDF algorithm using various parameter selections as compared with RK4.
numerical error of the proposed adaptive BDF method can be evaluated. For this reason, to evaluate the accuracy of the method, the well-known Runge-Kutta 4th order (RK4) method is used to solve the presented system of ODEs with a small constant time-step size. The selected time-step size is small enough that the truncation error remains below 10−8  . Then, the numerical results obtained from the RK4 method are used as a reference to assess the maximum and average error of the proposed adaptive BDF algorithm.
Since the oxygen mass transfer rate in Tank 1 is set at = 190 day−1 (Equation (8)), resulting in the highest order of stiffness in the system of ODEs, the dissolved oxygen concentration in Tank 1 is susceptible to showing higher sensitivity to the applied numerical method. Therefore, relative error is calculated based on the following equation:
where and show oxygen concentrations at Tank 1, obtained by the proposed algorithm and RK4, respectively. Table 2 lists the relative maximum (emax) and average (eavg) of for four sets of adjustable parameters. The accuracy assessment results show acceptable accuracy for the proposed algorithm. Specifically, for the four cases, the maximum and average errors contain less than 2 and 0.3 percent error, respectively.
In many applications, it is important to enhance the computational efficiency of ASM solvers. This paper presents a novel adaptive BDF algorithm to solve the system of ODEs arising from ASM models. The proposed method uses the partial evaluation of inverse Jacobian matrix in the BDFs and monitoring the convergence of the NR method to select the size of the time intervals. The method contains three adjustable parameters, including inflation and depression rates of the time-step and the maximum allowable number of iterations for the NR method, to control how time-step size varies through the course of the numerical solution. These three adjustable parameters can be optimized to reduce the computation time while maintaining acceptable accuracy. The algorithm is applied to solve a system of 45 ODEs representing an activated sludge model containing 15 constituents in each of 3 reactors. The sensitivity analysis is conducted, and the optimum adjustable parameters of the adaptive algorithm are identified. The accuracy analysis shows an acceptable level of accuracy for the proposed algorithm.
The authors would like to thank the Innovation section of the DC Water and Sewer Authority (DCWASA) for providing the data and financially support the research.
 Gernaey, K.V., van Loosdrecht, M.C.M., Henze, M., Lind, M. and Jørgensen, S.B. (2004) Activated Sludge Wastewater Treatment Plant Modelling and Simulation: State of the Art. Environmental Modelling & Software, 19, 763-783.
 Hauduc, H., Rieger, L., Ohtsuki, T., Shaw, A., Takacs, I., Winkler, S., Heduit, A., Vanrolleghem, P.A. and Gillot, S. (2011) Activated Sludge Modelling: Development and Potential Use of a Practical Applications Database. Water Science & Technology, 63, 2164-2182.
 Alikhani, J., Al-Omari, A., Clippeleir, H.D., Murthy, S., Takacs, I. and Massoudieh, A. (2016) Assessment of the Endogenous Respiration Rate and the Observed Biomass Yield for Methanol-Fed Denitrifying Bacteria under Anoxic and Aerobic Conditions. Water Science and Technology, In Press. http://dx.doi.org/10.2166/wst.2016.486
 Le Moullec, Y., Potier, O., Gentric, C. and Leclerc, J.P. (2011) Activated Sludge Pilot Plant: Comparison between Experimental and Predicted Concentration Profiles Using Three Different Modelling Approaches. Water Research, 45, 3085-3097.
 Busch, J., Elixmann, D., Kühl, P., Gerkens, C., Schlöder, J.P., Bock, H.G. and Marquardt, W. (2013) State Estimation for Large-Scale Wastewater Treatment Plants. Water Research, 47, 4774-4787. http://dx.doi.org/10.1016/j.watres.2013.04.007
 Hauduc, H., Rieger, L., Oehmen, A., van Loosdrecht, M.C., Comeau, Y., Heduit, A., Vanrolleghem, P.A. and Gillot, S. (2013) Critical Review of Activated Sludge Modeling: State of Process Knowledge, Modeling Concepts, and Limitations. Biotechnology and Bioengineering, 110, 24-46. http://dx.doi.org/10.1002/bit.24624
 Han, H.-G., Qian, H.-H. and Qiao, J.-F. (2014) Nonlinear Multiobjective Model-Predictive Control Scheme for Wastewater Treatment Process. Journal of Process Control, 24, 47-59.
 Henze, M., Leslie Grady Jr, C.P., Gujer, W., Marais, G.V.R., Matsuo, T., Henze, M., Leslie Grady Jr, C.P., Gujer, W., Marais, G.V.R. and Matsuo, T. (1987) A General Model for Single-Sludge Wastewater Treatment Systems. Water Research, 21, 505-515.
 Alikhani, J., Takacs, I., Al-Omari, A., Murthy, S., Mokhayeri, Y. and Massoudieh, A. (2015) Time Series Analysis for Uncertainty Evaluation of Performance of Biological Wastewater Treatment. Proceedings of the Water Environment Federation, 15, 3519-3526.
 Alikhani, J., Takacs, I., Omari, A.A., Murthy, S., Mokhayeri, Y. and Massoudieh, A. (2014) Inverse Modeling of Nitrification-Denitrification Processes: A Case Study on the Blue Plains Wastewater Treatment Plant in Washington, DC. Proceedings of the Water Environment Federation, 19, 4556-4567. http://dx.doi.org/10.2175/193864714815942620
 Lindberg, C.-F. and Carlsson, B. (1996) Adaptive Control of External Carbon Flow Rate in an Activated Sludge Process. Water Science and Technology, 34, 173-180.
 Cash, J.R. (2000) Modified Extended Backward Differentiation Formulae for the Numerical Solution of Stiff Initial Value Problems in ODEs and DAEs. Journal of Computational and Applied Mathematics, 125, 117-130. http://dx.doi.org/10.1016/S0377-0427(00)00463-5
 Hojjati, G., Rahimi Ardabili, M.Y. and Hosseini, S.M. (2004) A-EBDF: An Adaptive Method for Numerical Solution of Stiff Systems of ODEs. Mathematics and Computers in Simulation, 66, 33-41. http://dx.doi.org/10.1016/j.matcom.2004.02.019
 Jannelli, A. and Fazio, R. (2006) Adaptive Stiff Solvers at Low Accuracy and Complexity. Journal of Computational and Applied Mathematics, 191, 246-258.
 Akinfenwa, O.A., Jator, S.N. and Yao, N.M. (2013) Continuous Block Backward Differentiation Formula for Solving Stiff Ordinary Differential Equations. Computers & Mathematics with Applications, 65, 996-1005. http://dx.doi.org/10.1016/j.camwa.2012.03.111
 Celaya, E.A., Aguirrezabala, J.J.A. and Chatzipantelidis, P. (2014) Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s. Procedia Computer Science, 29, 1014-1026. http://dx.doi.org/10.1016/j.procs.2014.05.091
 Alikhani, J., Shayegan, J. and Akbari, A. (2015) Risk Assessment of Hydrocarbon Contaminant Transport in Vadose Zone as It Travels to Groundwater Table: A Case Study. Advances in Environmental Technology, 2, 77-84.
 Alikhani, J., Deinhart, A., Visser, A., Bibby, R., Purtschert, R., Moran, J., Massoudieh, A. and Esser, B. (2016) Nitrate Vulnerability Projections from Bayesian Inference of Multiple Groundwater Age Tracers. Journal of Hydrology, In Press.
 Shoghli, B. and Lim, Y.H. (2015) Predictive Scheme for Failures of Embankment Dams with High Overtopping Potential: Case Studies Using Remote Sensing Data. World Environmental and Water Resources Congress, Austin, 17-21 May 2015, 1332-1340.
 Shoghli, B., Lim, Y.H. and Alikhani, J. (2016) Evaluating the Effect of Climate Change on the Design Parameters of Embankment Dams: Case Studies Using Remote Sensing Data. World Environmental and Water Resources Congress, West Palm Beach, 22-26 May 2016, 575-585. http://dx.doi.org/10.1061/9780784479872.059
 Zhang, X., Scheving, B., Shoghli, B., Zygarlicke, C. and Wocken, C. (2015). Quantifying Gas Flaring CH4 Consumption Using VIIRS. Remote Sensing, 7, 9529-9541.
 Sharifi, S., Murthy, S., Takács, I. and Massoudieh, A. (2014) Probabilistic Parameter Estimation of Activated Sludge Processes Using Markov Chain Monte Carlo. Water Research, 50, 254-266. http://dx.doi.org/10.1016/j.watres.2013.12.010
 Metzger, M. (2000) A Comparative Evaluation of DRE Integration Algorithms for Real-Time Simulation of Biologically Activated Sludge Process. Simulation Practice and Theory, 7, 629-643. http://dx.doi.org/10.1016/S0928-4869(99)00015-4