Towards model predictive control: online predictions of ammonium and nitrate removal by using a stochastic ASM

Mathematical modelling of water resource recovery facilities (WRRFs) is a widely established discipline for research, plant design, optimization, simulation of process control strategies, etc. For these purposes there are many models to choose between, such as the activated sludge models (ASM), the anaerobic digestion models (ADM), the University of Cape Town (UCT) model or the TU Delft phosphorous removal model (Wentzel et al. 1992; Henze et al. 2000; Batstone et al. 2002; Gernaey et al. 2004; Meijer 2004; Hu et al. 2007). These complex models have differences in focus and as a result, in their structure. Hence choosing a model structure is, as with all modelling tasks, crucial to the outcome of the project. One important thing to include in the choice of a suitable model is the number of states and parameters. On the one hand more states and parameters leads to a more detailed model. However, on the other hand more details introduce more inputs that need to be distinguished and therefore estimated, measured or, if this is not possible, guessed at (Vanrolleghem et al. 1995). Furthermore, numerically solving large models with many states leads to long simulation times which can be demanding for data-driven optimizations, which need to be run at short time intervals (seconds to minutes). Although not yet used in the online operation of WWRFs, models can also be used to forecast future variables of interest for use in model predictive control (MPC), which means they should be fast and adaptable to online data.

The activated sludge model number 1 (ASM1) (Henze et al. 1987) describes organic matter degradation, nitrification and denitrification in activated sludge bioreactors. The model contains 13 states variables and 19 parameters. One of the most important challenges in using ASM1 is arguably attributing the many stoichiometric and kinetic parameters (Gernaey et al. 2004). The information needed for the characterization of these can come from three sources (Petersen et al. 2002): (1) default values from the literature, (2) full-scale plant data such as those collected by online sensors, and (3) information obtained from laboratory-scale experiments. The type of data and calibration framework to use is highly dependent on the intended use (e.g. Petersen et al. 2002). While (1) might be good for educational purposes or comparison of control strategies (e.g. Gernaey et al. 2014), optimization of processes with respect to a specific plant requires (2) and/or (3) (Petersen et al. 2002).

MPC aims to predict processes as a function of potential control actions and then choose the best control scenario based on the optimization of some objective function. In WRRFs this can translate to real-time modelling and forecasting of plant performance based on aeration control, optimizing electricity costs and effluent. When it comes to the selection of a suitable model for WRRF MPC strategies, the structure of states and parameters becomes particularly important. This is for the following two reasons: firstly, because parameters should be statistically identifiable from online data to take proper advantage of the real-time setting and secondly, because the computational requirements should be sufficiently low to allow for real-time, recursive simulation of several control scenarios. This means that a good online model should not have strong correlations between parameters, which is the case for parameters of ASM1 (Sharifi et al. 2014). Furthermore, the calibration should only depend on online data to avoid delays in updating the model, and hence the recalibration routine should not depend on information obtained from laboratory-scale experiments.

ASM1 has already been simplified to a linearized version to provide faster and yet reliable predictions (Smets et al. 2003). Furthermore, ASM1 has been reduced in an effort to make more parsimonious models (e.g. Mulas et al. 2007; Cadet 2014). However, the focus of these models is not online operations, i.e. to be updated with online data only. Online model applications are here managed by the use of stochastic, data-driven modelling (DDM) techniques. Many DDM methods exist, depending on the purpose and data availability (e.g. Dürrenmatt & Gujer 2012) and they are generally good alternatives when mechanistic models are not available or not valid (Gernaey et al. 2004). Since the detailed mechanistic understanding of the activated sludge process (ASP) already exists, the use of DDMs would ignore all the existing empirical process knowledge of nonlinearities and correlations.

Discretely observed stochastic differential equation (SDE) based models are often referred to as stochastic ‘grey-box’ (GB) models because the structure of the models represents both the physical/chemical/biological, deterministic (‘white-box’) understanding of the processes and the statistical, stochastic (‘black-box’) information indicated by data. Parameter calibration can be managed in the SDE-GB model by e.g. combining extended Kalman filtering (EKF) techniques and maximum likelihood estimation. This can be done statistically directly from online data by using e.g. the frameworks suggested by Kristensen et al. (2004), Tullekin (1993) or Jazwinski (1970). Furthermore, the EKF allows for optimal state estimation and handles additive noise effectively. Thordarson et al. (2012), Del Giudice et al. (2015) and Carstensen (1994) concluded that in terms of process prediction and control, SDE-GB models of the wastewater processes perform significantly better than traditional black-box models, such as ARMAX models, and also used them to statistically identify Monod-kinetic parameters from online measurements in an ASP (Carstensen et al. 1995). SDE-GB models have also been used to model incoming ammonium loads and first-flush phenomena (Bechmann et al. 2000; Halvgaard et al. 2017), and to forecast rainfall-runoff flow and volume in sewer systems for use in real-time optimization (Thordarson et al. 2012; Löwe et al. 2016).

In this paper we show that ASM1 can be rewritten as a simpler SDE-GB model that is applicable to online MPC purposes by treating state variables that show only slow and minor changes over short time horizons as model parameters that are kept fixed or intermittently re-estimated using online data. Thus, changes that occur slowly over weeks or months, such as changes in biomass, temperature, maintenance, wastewater composition etc., will be included in the parameters that are re-estimated intermittently with data from the previous few days. The small error introduced by this simplification is estimated by a stochastic diffusion process and consequently it can be managed in the control setup. Following this methodology it is possible to create a stochastic ASM with only three states representing ammonium concentration, nitrate concentration and available oxygen. This model can then be used to optimize the ammonium and nitrate removal process within an MPC approach.

This article presents a simple ASM based on SDEs, which uses flow and aeration data as input and assimilates online measurements of ammonium and nitrate to update model states and thus prepare for providing the best possible forecasts at each time step. The model gives reliable online forecasts of the ammonium and nitrate removal process from a few minutes to up to 24 hours ahead and considers measurement errors. It was developed and tested with data from a small recirculation WRRF with alternating operation. The simplicity of the model makes it a general tool that can be useful in recirculation facilities with different configurations without changing the model setup.

The model is developed and tested with data from Nørre Snede WRRF, which is located in central Jutland, Denmark. The plant is designed to handle a maximum capacity of 9,700 population equivalent (PE) and the current load is approximately 4,000 PE.

The WRRF includes several typical treatment processes that the wastewater goes through before discharge. Listed in order from when the wastewater enters the process, these are pretreatment, grit removal and grease trap, chemical dosage, nitrification/denitrification and secondary treatment. The nitrification/denitrification in the Nørre Snede plant happens in a process tank with a total volume of 3,500 m³. The tank is divided into three smaller chambers operating under different conditions. This is illustrated in Figure 1, which also shows that the aeration tank is equipped with nutrient sensors, aeration equipment, a recirculation pump and rotors that control the flow direction/velocity (direction shown with arrows in the figure; rotors are located at the bridge).

The facility is currently operated with a rule-based control strategy, as described by Isaacs & Thornberg (1998), Zhao et al. (2004) and Kim et al. (2014), for example. In this case the control switches aeration on/off as a function of online ammonium and nitrate concentration measurements. Therefore the conditions switch between anoxic and aerobic and the cycle lengths depend on the conditions in the process tanks. This is managed in the control platform STAR Utility Solutions™ (Sørensen et al. 1994; Nielsen & Önnerth 1995).

The current control of the plant (i.e. actuator settings controlling aeration and inlet flow) is updated every 2 min and as result, aeration and inflow data are available every 2 min. The control rules are based on ammonium and nitrate signals, which are only sampled every 5 min directly in the aeration tank, meaning that observations are sampled at a different frequency to the control sampling. Calibration of sensors happens automatically two to four times per day, resulting in 30–60 min with no new observations. There is a response time from when aeration starts/stops until this is observed by the sensors. This is due to hydraulics in the tank and processing time in the sensors (Rieger et al. 2003). This response time is estimated using the method suggested by Stentoft et al. (2017): an estimate from the point at which conditions are shifting to the point at which a change in the trend in measurements is observed. Flow data are available at the outlet (after the settler) and change between 0 and ∼45 m³/h because of a pumping scheme. To account for this scheme, flow data are filtered by a second-order Fourier series. The available data used in this study are summarized in Table 1.

This section includes a description of the SDE-GB model, a simplification of the ASM1 model with noise terms added, the inclusion of aeration control and inflow as input. Finally, the parameter estimation is briefly described.

Following the notation proposed by Corominas et al. (2010) the ordinary differential equations that govern ammonium, nitrate and oxygen in ASM1 can be written in a Gujer matrix, as presented in the Supplement (available with the online version of this paper), together with a brief description of ASM1 parameters and state variables. The complexity of these equations is considered an obstacle for use in a real-time setting since many of the variables are unmeasured and consequently constants would be difficult to distinguish. We therefore make simplifications to obtain a more suitable model. The main assumption in this simplification is that the model parameters will be re-estimated frequently, and therefore several state variables of ASM1 will become constant and some parameters will become unimportant.

• The rate, ρ₄, which governs ammonification of soluble organic nitrogen can be ignored. This is considered reasonable as the ammonification rate parameter k_am is typically small compared to the process rates of nitrification and denitrification (Henze et al. 1987).

• The state variable S_B is constant throughout the day. In practice it will follow a diurnal pattern similar to that of S_NH4, but since S_B is not measured these will be difficult to distinguish.

• The state variables governing active heterotrophic and autotrophic biomasses X_OHO and X_ANO can be considered constant on a daily basis, and hence can be treated as parameters. This is considered reasonable as the biomass is known to only change over longer periods of time.

• The parameter for the relative amounts of N/COD in biomass, i_{N_COD}, can be ignored. A stoichiometric calculation by Henze et al. (1987) (assuming a typical cell formation, C₅H₇O₂N) indicates that i_{N_COD} is 0.086. This is very small compared to 1/Y_ANO which is approximately 4.2, and therefore it will be difficult to estimate.

• The half-velocity parameters K_O2,OHO and K_O2,ANO for oxygen use are considered to be equal (K_O2,OHO = K_O2,ANO = K_O2), as Henze et al. (1987) argue that they are not, quantitatively, that different.

The stochastic model presented fits the general model structure for continuous-discrete stochastic state space models, i.e. a model of the state variables in continuous time and measurements of some of the states at discrete times. The R-package CTSM-R (Kristensen et al. 2004; Juhl et al. 2016) can manage just this kind of system, and is therefore used to estimate parameters and predict the effect of control. This paper provides only a brief summary of how the package works and how it is used here. For further information on this, see CTSM-R (2018).

The estimation setup implies that initial state and parameter estimates are necessary in the parameter estimation procedure. These can be supplied either as prior distributions or simply as estimates with some maximum and minimum boundaries. For most parameters these initial estimates are based on the literature (e.g. Henze et al. 1987; Henze & Comeau 2008). However, a few parameters are unnecessary or cannot be estimated. The initial state values S_NH4,0, S_NO3,0 and S_O,MO,0 can easily be estimated directly from data, as ammonium and nitrate are directly measured and the oxygen signal is known, and hence these are considered unnecessary to estimate. Furthermore, the parameters ⁠, θ₃ and θ₄ show very small deviations and notably correlation with other parameters. It is also argued by Henze et al. (1987) that does not need estimating. For these reasons these are kept constant here at [K_NO3,OHO, θ₃, θ₄] = [3.0, 5.0E − 6, 1.0], thereby reducing the number of parameters to estimate.

Firstly the model is qualitatively evaluated by comparing the model predictions with data and discussing parameter estimates. Secondly, the model is quantitatively evaluated by running it for 1 month and discussing the statistics of the residuals. We stress that the model is run ‘online’ in the sense that parameters are estimated only by minimizing the objective function described in the previous section. Furthermore, the states S_NH4, S_NO3 and S_O,MO, are updated using the EKF whenever a new measurement becomes available. Figure 2 shows an example of one prediction of ammonium at a given inlet flow and aeration signal. The state, S_NH4, is updated with present data and then predicted 2 hours ahead. Clearly, uncertainty increases with increasing forecast horizon.

Parameters are estimated with data from a period at the beginning of October 2016, chosen arbitrarily among periods without rain. The length of the parameter estimation period is 4 days and 4 hours (corresponding to 3,000 time steps of 2 min). These parameters are used to predict the concentrations of ammonium and nitrate in the aeration tank. Figures 3 and 4 show predictions of the ASP 60 time steps ahead corresponding to 2 hours, given the aeration signal. This is done for 24 hours, meaning that each time a new measurement becomes available a prediction similar to Figure 2 is made and compared with data. This is done during normal operation of the plant i.e. no rain and no (known) problems.

In Figure 3 it is evident that under normal operation, the modelled ammonium concentration follows the same dynamics as the data. During the aeration phase the ammonium concentration decreases, and when aeration is switched off, NH₄ increases. In periods when no new data are received (i.e. calibration of the ammonium sensor from 17:30 to 18:30), the model continues to provide reliable estimates. The nitrate concentrations estimated in Figure 3 also follow dynamics similar to those in the data. It is noted that when aeration is off, nitrate decreases and when it is on, it increases. However, a dynamic starting at 06:00 does not follow the behavior shown by the sensor measurements. This period contains a relatively long timespan without aeration that would normally mean denitrification, but in this case we see that nitrate increases. This could be due to some unmodelled dynamics, problems with a drifting nitrate sensor or an unusually large load of nitrate in the influent from industry, for example. Overall, the results show that the uncertainty of the nitrate predictions is greater than the uncertainty of the ammonium predictions, and hence larger deviations from the modelled concentrations are expected. Figure 4 shows the estimates of the unmeasured state, S_O,MO. It is plotted together with the measured dissolved oxygen (DO) concentrations and the setpoint of the actuator, O. It is clear that when the setpoint is lower, it takes a longer time for S_O,MO to reach the maximum level. The comparison between the aeration status and the measured oxygen concentrations highlight that short periods of aeration are not registered in the measurements. This is caused by the location of the sensor, which is opposite the aeration grid (see Figure 1). This supports the choice of not including the measured DO as input/state in the model. Therefore the actuator signal is considered superior as it reports all periods of added oxygen and furthermore, it does not have any response time from when air is added until it is observed in the tank.

The parameters that are estimated in the period mentioned previously are presented in Table 2.

Some parameters are difficult to evaluate as they represent some catchment/plant-specific information. Nonetheless, parameters are discussed here and in some cases compared with the literature:

• The residual errors (which are split into diffusion and measurement noise) are found to be similar to what is estimated for ammonium in the aeration tank by Halvgaard et al. (2017) in another plant with similar sensors and configuration. However, the uncertainty of the measuring equipment is much smaller compared to what is stated by the sensor manufacturer (see Table 1 footnote).

• The Monod kinetics parameter, K_NH4,ANO, is found to be similar to what is statistically estimated by Carstensen et al. (1995). Here it was found to be between 0.44 and 0.76 among different plants and catchments. However, it is lower than ∼1, which is what is suggested by Henze et al. (1987).

• The mean incoming concentrations of ammonium, μ_NH4,in and nitrate, μ_NO3,in seem reasonable as these are similar to what is typically found in a ‘catchment with little industrial activity’ according to Henze & Comeau (2008). The diurnal variations in ammonium, s₁, s₂, c₁ and c₂ are catchment specific. However, these are considered reasonable as they produce a variation that is comparable with the one presented in Henze & Comeau (2008), i.e. a similar shape with peaks in the morning and afternoon.

• The new parameters θ₁ and θ₂ are difficult to compare with the literature as they depend on many parameters and the states X_ANO and X_OHO. However, following the typical parameter suggestions from Henze et al. (1987), these should be 3.33 X_ANO and 0.82 X_OHO respectively.

• The parameters related to the incoming water, r_c and ρ, are estimated to be 7.21 × 10⁻⁴ and 1.08 × 10⁻⁵ respectively. Summing these and multiplying the mean flow with ρ we get 1.11 × 10⁻³. This is slightly more than the expected value of 3.43 × 10⁻⁴ which is found by dividing mean flow with the volume of the process tank. This difference can be due to the recirculation that happens between the nitrification and denitrification tanks.

• The relative oxygen transfer rate k₁ depends on many factors such as tank design (e.g. reactor geometry, aeration design), physico-chemical properties (e.g. liquid composition, viscosity, temperature) and the presence of biomass (e.g. Pittoors et al. 2014). Therefore it is difficult to determine empirically as it varies between facilities and over time.

The model's predictive ability is tested by re-estimating parameters every hour for a period of 1 month and 1 week, starting late September 2016. The model is then used to predict concentrations of ammonium and nitrate 1, 60 and 720 time steps ahead (corresponding to 2 min, 2 hours and 1 day, respectively). The predictions are compared with data and a 24-hour running mean absolute residual is calculated. Figure 5 illustrates how this changes over time for predictions 2 hours ahead (60 time steps). Table 3 shows the statistics of the mean absolute residual for all the different prediction horizons.

In Figure 5 it is evident that during some periods (i.e. October 2nd, October 15th and October 22nd) the uncertainty increases. Comparison with rain data supplied by the Danish Meteorological Institute (DMI 2018) shows that many of these periods are characterized by wet weather. This is also indicated in Table 3 where the general picture is that uncertainty increases during wet weather. In Table 3 it is also clear that the relative 2 min uncertainty is ∼2% for ammonium and ∼6% for nitrate. This is comparable with the sensor uncertainty listed in Table 1 and indicates that the 2 min predictions cannot be further improved even with a more detailed deterministic model. The 2-hour and 24-hour predictions perform worse than the sensor uncertainty. This can, on the one hand, indicate that there is room for improvement, but can, on the other hand, also mean that there is some stochastic behavior (e.g. incoming nutrients or biomass efficiency) which is more pronounced when predicting further ahead from the EKF state update. It should be added that the treatment requirements for Nørre Snede WRRF state that during a 24-hour period, the ammonium concentration should be <2 mgN/L and total N should be <8 mgN/L in the outlet. Consequently this means that accurately predicting ammonium is more important because nitrate only affects total N. The relative uncertainties of ammonium and nitrate of <10% and <20% respectively (in dry weather, 24 hours ahead) are considered sufficient for stochastic MPC.

The results in Table 3 are difficult to compare with full ASM models, as to our knowledge there is no framework for making the full ASM models online adaptive to data. However, our results can be compared with data-driven model parameter estimations of full ASMs (i.e. other methods that rely only on data from online ammonium/nitrate sensors). One example of such an approach is provided by Sin et al. (2008), where the parameters of an ASM2d model were estimated using only frequently sampled online ammonium, nitrate and oxygen measurements (sampled every 5 min, similar to this study) in the 50,000 PE Haaren WRRF in the Netherlands with alternating control of aeration. Parameters were calibrated using Monte Carlo simulations to minimize a weighted sum of squared errors (WSSE) based on a calibration period of 16,117 measurements (56 days). Model performance was compared with data in a validation period of 9,217 measurements (32 days). The results showed a mean absolute error (MAE) for ammonium of 1.39 mgN/L and 0.98 mgN/L in the calibration and validation periods, respectively. For nitrate concentrations an MAE of 2.56 mgN/L and 2.31 mgN/L were found. These error values are 10 times larger than what we found in this study, cf. Figure 3. Also, the framework presented in Sin et al. (2008) differs from this study in that the model states are not updated when new data become available and hence short time predictions in the validation period (up to 24 hours ahead) are not based on all the information available in an online situation. Additionally it is noted that the framework by Sin et al. (2008) required some computation time (45 min per simulation (in 2008) on a PC, and 500 Monte Carlo simulations where used to obtain a model), which makes it non-ideal for online applications.

The development of tools for online performance optimization of WRRFs using models is crucial for exploiting the full potential of digitalization. Hence, the development of robust approaches to online identifiable ASMs for improved short horizon predictions is needed. These models should also include additional processes, such as biological removal of COD and P, to achieve an overall improvement of all the removal processes in the plant. This paper provides a first step in this direction with online predictions of ammonium and nitrogen removal.

Grey-box models based on SDEs are efficient tools as they can estimate both processes and noise from real-time data. A stochastic model of an aeration tank is proposed here. The model contains a deterministic term consisting of both a simplified ASM and input functions determining the influence of control and inflow. The model is used to predict the nitrification/denitrification in Nørre Snede WRRF in Denmark as a function of aeration and inflow.

The results show that despite the simple structure of the proposed model, the dynamics of the nutrient concentrations are captured. Quantitative investigations show that the processes are predicted accurately, i.e. 24-hour predictions of the ammonium and nitrate concentrations in the aeration tank are predicted with relative errors of <10% and <20% respectively. Consequently this is considered to be a step towards stochastic MPC of water resource recovery processes.

This work was partly funded by the Innovation Fund Denmark (IFD) under File No. 7038-00097B – the first author's industrial PhD: Stochastic Predictive Control of Wastewater Treatment Processes. Data for this study were supplied by Ikast Brande Forsyning A/S.