An Adaptive TimeStep Backward Differentiation Algorithm to Solve Stiff Ordinary Differential Equations : Application to Solve Activated Sludge Models

A backward differentiation formula (BDF) has been shown to be an effective way to solve a system of ordinary differential equations (ODEs) that have some degree of stiffness. However, sometimes, due to high-frequency variations in the external time series of boundary conditions, a small time-step is required to solve the ODE system throughout the entire simulation period, which can lead to a high computational cost, slower response, and need for more memory resources. One possible strategy to overcome this problem is to dynamically adjust the time-step with respect to the system’s stiffness. Therefore, small time-steps can be applied when needed, and larger time-steps can be used when allowable. This paper presents a new algorithm for adjusting the dynamic time-step based on a BDF discretization method. The parameters used to dynamically adjust the size of the time-step can be optimally specified to result in a minimum computation time and reasonable accuracy for a particular case of ODEs. The proposed algorithm was applied to solve the system of ODEs obtained from an activated sludge model (ASM) for biological wastewater treatment processes. The algorithm was tested for various solver parameters, and the optimum set of three adjustable parameters that represented minimum computation time was identified. In addition, the accuracy of the algorithm was evaluated for various sets of solver parameters.


Introduction
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) [1] [2] [3].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 [4] [5] [6] [7].Henze et al. [8] proposed the first ASM (ASM1), and since then, several other versions of ASMs have been proposed [9].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 [10], and the use of parameter estimation algorithms that require numerous simulations.For example, Alikhani et al. [11] 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 [12] 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 [13] [14] [15] [16].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. [17] 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 [18], Bayesian parameter estimation framework [19], and other application of numerical methods (e.g.[20] [21] [22]).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.

Adaptive Backward Differentiation Algorithm
The system of ODEs can be written in the following general form [23]: where, t is the time, ( ) is the state variables vector, and 0 C 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 [17] [24]: 1.5 1.5 0.5 0.5 0 where i indicates the time-step with non-constant step size i h , ( ) represents the discretized form of Equation (1) when ( ) ( ) . In Equations ( 2) and (3),

F
are the first and second order BDFs considering non-constant step size at each time-step, respectively, which will be henceforth denoted as F for the sake of generality.The conventional Newton-Raphson (NR) me- thod for solving the system of non-linear F to find

C
can be formulated as: where k is the iteration index beginning from 0, 1 l J − 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, l J in Equation ( 4) is a recycled Jacobian matrix evaluated at some previous time l t .
( ) , , where n is the size of the ODEs, ij δ is the Kronecker Delta, 1st l J and 2nd l J denote the Jacobian matrix for the first and second order F , respectively, where K l c is the K th element of vector l C , and similarly, I f is the I th element of vector f at l t .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 [25]: where The adaptive time interval scheme adjusts the size of time-steps based on monitoring the convergence of the NR method [25] [26].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 k max , and this inverse of Jacobian matrix ( 1 l J − ) will be used in Equation ( 4).
The iteration continues until the number of iterations to achieve convergence in the NR method exceeds k max .At this time, a new Jacobian matrix will be computed and/or the time-step size will be reduced by the factor of 1 depicts a de- tailed flowchart of this algorithm, in which ρ, γ, and k max 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.

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 [11].The total volume of the reactor was 17,500 m 3 , 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 , Hydrolysis of entrapped organic nitrogen

Governing ODEs
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 [9] [23]: where V (L 3 ), C (M•L −1 ), Q (L 3 •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, * C (M•L −1 ) is the saturated concentration, and m  (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, N m is the total number of tanks, and N r is the total number of reactions.
, m m Q ′ is the flow rate from tank m' to tank m, due either to sequential stages or through feedback or bypass (a positive , m m Q ′ 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 ( , j ret C ), a dynamic clarifier model [27] 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.

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 k max = 10.

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 (k max ).A small initial time-step of h 0 = 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 k max 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 , as expected, because the algorithm needs smaller timesteps to resolve the higher disturbance in the unfiltered ( ) . Nevertheless, the difference in the running time is not significant.

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 k max = 10, the time-step can reach a maximum value of 0.35 days, whereas, at k max = 4, the maximum value is 0.14 days.In other words, a higher k max 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 k max depends on 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 k max 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 k max , 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.

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  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 4 th 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 [26].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 2 O F = 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:  2 lists the relative maximum (e max ) and average (e avg ) of ( ) e t 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.

Conclusion
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 . 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 2 .
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 inTable1, respectively.

Figure 3 .
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.

Figure 4 .
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.

Figure 5 .
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 .
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).

Figure 7 .
Figure 7. Adaptive time-step (dt) variation for 120 days of applying adaptive BDF method on unfiltered U(t) for different adjustable parameters.
concentrations at Tank 1, obtained by the proposed algorithm and RK4, respectively.Table

Table 2 .
Maximum and the average error of the ABDF algorithm using various parameter selections as compared with RK4.