Abstract

Online model predictive control (MPC) of water resource recovery facilities (WRRFs) requires simple and fast models to improve the operation of energy-demanding processes, such as aeration for nitrogen removal. Selected elements of the activated sludge model number 1 modelling framework for ammonium and nitrate removal were included in discretely observed stochastic differential equations in which online data are assimilated to update the model states. This allows us to produce model-based predictions including uncertainty in real time while it also reduces the number of parameters compared to many detailed models. It introduces only a small residual error when used to predict ammonium and nitrate concentrations in a small recirculating WRRF facility. The error when predicting 2 min ahead corresponds to the uncertainty from the sensors. When predicting 24 hours ahead the mean relative residual error increases to ∼10% and ∼20% for ammonium and nitrate concentrations respectively. Consequently this is considered a first step towards stochastic MPC of the aeration process. Ultimately this can reduce electricity demand and cost for water resource recovery, allowing the prioritization of aeration during periods of cheaper electricity.

INTRODUCTION

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.

CASE STUDY: NØRRE SNEDE WRRF

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.

Operation and design

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 m3. 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).

Figure 1

Overview of the process tank at Nørre Snede WRRF with important features labelled.

Figure 1

Overview of the process tank at Nørre Snede WRRF with important features labelled.

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

Data

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 m3/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.

Table 1

Overview of online data used in this study

SymbolDescriptionSample frequencyUnitUncertainty
Effluent Flow 2 min m3/L Unknown 
Actuator setting 2 min mgO2/L 
MsNH Measured ammonium 5 min mgN/L ±3% 
MsNO Measured nitrate 5 min mgN/L ±5% 
SymbolDescriptionSample frequencyUnitUncertainty
Effluent Flow 2 min m3/L Unknown 
Actuator setting 2 min mgO2/L 
MsNH Measured ammonium 5 min mgN/L ±3% 
MsNO Measured nitrate 5 min mgN/L ±5% 

The uncertainty is based on the information available from the sensor manuals (HACH Lange Aps 2013, 2014).

THEORY AND METHODOLOGY

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.

Stochastic grey-box model

A discretely observed stochastic grey-box model can be written on the state-space form as 
formula
 
formula
where the description of the dynamics of the states are divided into a (deterministic) drift term and a (stochastic) diffusion term . The system is observed only through which is linked to via the observation equation . The residual error is separated in to two terms. Diffusion, represents model approximations and undescribed noise disturbances, i.e. changes in biomass efficiency, unmodelled inflow, uncertainty of input variables , or true stochastic behaviour in the processes. Measurement noise, represents the noise occurring due to imperfect accuracy of the measuring equipment, i.e. measurement uncertainty in the ammonium and nitrate sensors.

Simplification of ASM1

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, ρ4, which governs ammonification of soluble organic nitrogen can be ignored. This is considered reasonable as the ammonification rate parameter kam is typically small compared to the process rates of nitrification and denitrification (Henze et al. 1987).

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

  • The state variables governing active heterotrophic and autotrophic biomasses XOHO and XANO 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, iN_COD, can be ignored. A stoichiometric calculation by Henze et al. (1987) (assuming a typical cell formation, C5H7O2N) indicates that iN_COD is 0.086. This is very small compared to 1/YANO which is approximately 4.2, and therefore it will be difficult to estimate.

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

Applying these assumptions, the new, simple model of the ASP can be identified. The shorthand notations αNH4, αNO3 and αO2 refer to the changes in concentration of ammonium, nitrate and relative oxygen respectively. 
formula
The half-saturation constant KcNH4,ANO is introduced because the state SO,MO is the Monod term indicating how quickly the process is running relative to the maximum rate 
formula
The seven new parameters to estimate online are therefore θi, KNO3,OHO, KNH4,ANO and KcNH4,ANO where θi relates to the original ASM1 parameters as 
formula
where C1 and C2 are correction factors that are introduced because SMO,O is the relative amount of oxygen.

Aeration control and inflow

For the purpose of using the model for MPC of N-removal, it is necessary to include the effect of aeration and incoming wastewater as external inputs in the model. The signal determining the intensity of aeration and measurements of incoming wastewater flow are available online and consequently the control should be a function of these, i.e. C_(O,Q) is a function that describes the effect of inflow and aeration control on the given state. These functions are here determined from the literature. More specifically, the two-films theory (Lewis & Whitman 1924) and diurnal variations in ammonium concentration and constant (low) nitrate concentrations in the incoming wastewater (Henze & Comeau 2008). This means that CNH4(O,Q) and CNO3(O,Q) are given as 
formula
where μNH4,in, μNO3,in, si and ci are parameters relating to the inflow. Note that μNO3,in is typically ∼0 (e.g. Henze & Comeau 2008). The parameters rc and ρ are related to the recirculation to and the volume of the aeration tank (see Figure 1). The control of the aeration is given as 
formula
where k1 is a transfer constant (Lewis & Whitman 1924) related to the size and efficiency of the aeration equipment. The maximum value of the relative oxygen state is 1, and hence SO,MOmax is set to 1 and should not be estimated.

Stochastic ASM

A three-state grey-box model governing the ammonium and nitrate concentrations in the aeration tank can be written as 
formula
 
formula
 
formula
where the deterministic terms α and C govern the ASP, the aeration and the inflow as described in previous sections. To avoid negative noise and to make estimation of small noise processes easier, the diffusion terms are estimated as exponential parameters (i.e. ). The system is discretely observed through ammonium and nitrate sensors in the aeration tank. The measurements (MsNH and MsNO) from these relate to the system as 
formula
where εNH4 and εNO3 are independent and identically distributed (i.i.d.) Ɲ(0, 1), i.e. the residuals of the measurements are normally distributed with zero mean and exp(s1,NH4), exp(s2,NO3) standard deviation. The changes in the states dSNH4, dSNO3 and dSMO,O are given as state variables where ɷ1, ɷ2 and ɷ3 are 1-dimensional standard Wiener processes and exp(s11), exp(s22) and exp(s33) represent the deviation of these processes.

Online parameter estimation in stochastic ASM

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 parameter estimates are based on a maximum likelihood method, by assuming Gaussian-distributed conditional probability densities. 
formula
where εt=yt – ŷt|t-1 (ŷt|t-1=E(yt|yt-1, θ)) and Rt|t-1=V(yt|yt-1, θ). εt and Rt|t-1 are computed by means of a version of the EKF (e.g. Jazwinski 1970). The likelihood function can be simplified to a simpler log-likelihood function by conditioning on y0 and taking the negative logarithm. However, this rewriting is omitted here. The parameter estimates are then obtained by minimizing this log-likelihood.
The EKF mentioned above is a continuous-discrete time version of the EKF. With initial conditions for the model values and variance estimate (1|0 and P1|0) the filter approximations of the output predictions are given as 
formula
which here translates to 
formula
Here C is the Jacobian of the observation equation, h. The innovation given by 
formula
The Kalman gain, Kk, for the EKF is then calculated as 
formula
and the system is updated 
formula
is the state estimates (SNH4, SNO3 and SO,MO). The state prediction is done numerically by 
formula
where A(t) is the Jacobian of the drift term fi(t,…). This Jacobian is calculated using a method based on Speelpenning (1980). In calculations of the Jacobian it is assumed that x = k|k-1, u = uk, t = tk and the parameters, θ, are known. The ordinary differential equations (ODEs) are solved by numerical integration schemes suggested by Hindmarsch (1983) (cited in Kristensen & Madsen 2003, p. 17). This is to ensure an intelligent re-evaluation of A and σ. From this construction we see that the approximation is only good when nonlinearities are not too strong.

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 SNH4,0, SNO3,0 and SO,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 , θ3 and θ4 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 [KNO3,OHO, θ3, θ4] = [3.0, 5.0E − 6, 1.0], thereby reducing the number of parameters to estimate.

RESULTS AND DISCUSSION

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 SNH4, SNO3 and SO,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, SNH4, is updated with present data and then predicted 2 hours ahead. Clearly, uncertainty increases with increasing forecast horizon.

Figure 2

An example of a 2-hour prediction of ammonium concentrations (which is run every 2 min in the online setup). The uncertainty increases further into the future, as estimated by the SDE.

Figure 2

An example of a 2-hour prediction of ammonium concentrations (which is run every 2 min in the online setup). The uncertainty increases further into the future, as estimated by the SDE.

Model dynamics

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.

Figure 3

2-hour predictions (60 timesteps of 2 mins) of the ammonium concentration in the aeration tank (SNH4) with the measured concentrations (top) and of the nitrate concentration in the aeration tank (SNO3) with the measured concentrations (bottom). Note that the y-axes differ because there is more variation in the nitrate observations.

Figure 3

2-hour predictions (60 timesteps of 2 mins) of the ammonium concentration in the aeration tank (SNH4) with the measured concentrations (top) and of the nitrate concentration in the aeration tank (SNO3) with the measured concentrations (bottom). Note that the y-axes differ because there is more variation in the nitrate observations.

Figure 4

Top: the input, O. Middle: the estimated Monod oxygen state, SO,MO. Bottom: the measured oxygen in the aeration tank (solid line) and the binary signal for aeration on/off (dotted line).

Figure 4

Top: the input, O. Middle: the estimated Monod oxygen state, SO,MO. Bottom: the measured oxygen in the aeration tank (solid line) and the binary signal for aeration on/off (dotted line).

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, NH4 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, SO,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 SO,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.

Parameter estimates

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

Table 2

Parameter estimates and standard deviations

ParameterEstimateStandard deviation
Simplified ASM 
KNH4,ANO [mgN/L] 4.80 × 10−1 1.80 × 10−1 
KCNH4,ANO [mgN/L] 1.05 × 10−3 6.99 × 10−4 
θ1 [mgN/L] 5.00 × 10−2 5.39 × 10−3 
θ2 [mgN/L] 2.47 × 10−1 1.59 × 10−2 
Aeration control and inflow 
s1 [mgN/L] 1.14 6.63 × 10−1 
s2 [mgN/L] −5.63 × 10−1 3.73 × 10−1 
c1[mgN/L] 4.25 × 10−1 3.24 × 10−1 
c2 [mgN/L] −5.82 × 10−2 1.75 × 10−1 
μNH4,in [mgN/L] 1.35 × 10−1 6.76 
μNO3,in [mgN/L] 1.10 × 10−2 2.01 × 10−2 
K1 [L/mgO] 1.12 × 10−1 6.48 × 10−3 
ρ[] 1.08 × 10−5 1.64 × 10−5 
rc [] 7.21 × 10−4 4.31 × 10−4 
Noise and diffusion terms 
s11 [log(mgN/L)] −4.79 7.65 × 10−2 
s22 [log(mgN/L)] −3.20 2.06 × 10−2 
s33 [] −2.86 2.16 × 10−1 
s1,NH [log(mgN/L)] −8.66 1.26 × 10−1 
s1,NO [log(mgN/L)] −1.80 × 10−1 2.56 × 10−1 
ParameterEstimateStandard deviation
Simplified ASM 
KNH4,ANO [mgN/L] 4.80 × 10−1 1.80 × 10−1 
KCNH4,ANO [mgN/L] 1.05 × 10−3 6.99 × 10−4 
θ1 [mgN/L] 5.00 × 10−2 5.39 × 10−3 
θ2 [mgN/L] 2.47 × 10−1 1.59 × 10−2 
Aeration control and inflow 
s1 [mgN/L] 1.14 6.63 × 10−1 
s2 [mgN/L] −5.63 × 10−1 3.73 × 10−1 
c1[mgN/L] 4.25 × 10−1 3.24 × 10−1 
c2 [mgN/L] −5.82 × 10−2 1.75 × 10−1 
μNH4,in [mgN/L] 1.35 × 10−1 6.76 
μNO3,in [mgN/L] 1.10 × 10−2 2.01 × 10−2 
K1 [L/mgO] 1.12 × 10−1 6.48 × 10−3 
ρ[] 1.08 × 10−5 1.64 × 10−5 
rc [] 7.21 × 10−4 4.31 × 10−4 
Noise and diffusion terms 
s11 [log(mgN/L)] −4.79 7.65 × 10−2 
s22 [log(mgN/L)] −3.20 2.06 × 10−2 
s33 [] −2.86 2.16 × 10−1 
s1,NH [log(mgN/L)] −8.66 1.26 × 10−1 
s1,NO [log(mgN/L)] −1.80 × 10−1 2.56 × 10−1 

Parameters are calculated using 3,000 timesteps of data corresponding to 4 days and 4 hours. The period is from the beginning of October 2016.

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, KNH4,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, s1, s2, c1 and c2 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 θ1 and θ2 are difficult to compare with the literature as they depend on many parameters and the states XANO and XOHO. However, following the typical parameter suggestions from Henze et al. (1987), these should be 3.33 XANO and 0.82 XOHO respectively.

  • The parameters related to the incoming water, rc and ρ, are estimated to be 7.21 × 10−4 and 1.08 × 10−5 respectively. Summing these and multiplying the mean flow with ρ we get 1.11 × 10−3. This is slightly more than the expected value of 3.43 × 10−4 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 k1 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.

Model performance

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.

Table 3

Summary statistics of the residuals of the model predictions, based on residuals obtained from using the model with data from Nørre Snede WRRF for a period of 1 month in September/October 2016

PrecipitationaPredictionMean absolute error [mgN/L]
Relative absolute error [–]
NoYesNoYes
Ammonium 2 min 0.0264 0.0278 0.0212 0.0223 
2 hours 0.0884 0.1065 0.0709 0.0855 
24 hours 0.0925 0.129 0.0855 0.104 
Nitrate 2 min 0.0804 0.081 0.0594 0.0598 
2 hours 0.212 0.2399 0.157 0.177 
24 hours 0.246 0.374 0.1815 0.276 
PrecipitationaPredictionMean absolute error [mgN/L]
Relative absolute error [–]
NoYesNoYes
Ammonium 2 min 0.0264 0.0278 0.0212 0.0223 
2 hours 0.0884 0.1065 0.0709 0.0855 
24 hours 0.0925 0.129 0.0855 0.104 
Nitrate 2 min 0.0804 0.081 0.0594 0.0598 
2 hours 0.212 0.2399 0.157 0.177 
24 hours 0.246 0.374 0.1815 0.276 

‘Precipitation’ indicates whether only dry weather periods (No) or all periods are considered (Yes).

aWet weather periods are as follows: 28 September to 2 October, 15–19 October, 22–23 October.

Figure 5

Mean absolute error for the 2-hour (60 timesteps of 2 mins) prediction of ammonium and nitrate. Plotted with daily precipitation data from DMI (2018).

Figure 5

Mean absolute error for the 2-hour (60 timesteps of 2 mins) prediction of ammonium and nitrate. Plotted with daily precipitation data from DMI (2018).

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.

Towards MPC – simplified and full ASM models

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.

CONCLUSION

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.

ACKNOWLEDGEMENTS

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.

REFERENCES

REFERENCES
Batstone
D. J.
,
Keller
J.
,
Angelidaki
R. I.
,
Kalyuzhnyi
S. V.
,
Pavlostathis
S. G.
,
Rozzi
A.
,
Sanders
W. T. M.
,
Siegrist
H.
&
Vavilin
V. A.
2002
Anaerobic Digestion Model No. 1. Scientific and Technical Report No. 13
.
IWA Publishing
,
London
,
UK
.
Bechmann
H.
,
Madsen
H.
,
Kjølstad Poulsen
N.
&
Nielsen
M. K.
2000
Grey box modelling of first flush and incoming wastewater at a wastewater treatment plant
.
Environmetrics
2000
(
11
),
1
12
.
Cadet
C.
2014
Simplifications of Activated Sludge Model with preservation of its dynamic accuracy
.
IFAC
47
(
3
),
7134
7139
.
Carstensen
J.
1994
Identification of Wastewater Treatment Processes
.
PhD Thesis
,
IMSOR, Technical University of Denmark
.
Carstensen
J.
,
Harremoës
P.
&
Madsen
H.
1995
Statistical identification of monod-kinetic parameters from on-line measurements
.
Water Science and Technology
31
(
2
),
125
133
.
Corominas
L. L.
,
Rieger
L.
,
Takács
I.
,
Ekama
G.
,
Hauduc
H.
,
Vanrolleghem
P. A.
,
Oehmen
A.
,
Gernaey
K. V.
,
van Loosdrecht
M. C. M.
&
Comeau
Y.
2010
New framework for standardized notation in wastewater treatment modelling
.
Water Science and Technology
61
(
4
),
841
857
.
CTSM-R
2018
CTSM-R – Continuous Time Stochastic Modelling for R
.
http://ctsm.info, visited 31 January 2018
.
Del Giudice
D.
,
Lowe
R.
,
Madsen
H.
,
Mikkelsen
P. S.
&
Rieckermann
J.
2015
Comparison of two stochastic techniques for reliable urban runoff prediction by modeling systematic errors
.
Water Resources Research
51
(
7
),
5004
5022
.
DMI
2018
Vejrarkiv
. .
Dürrenmatt
D.
&
Gujer
W.
2012
Data-driven modeling approaches to support wastewater treatment plant operation
.
Environmental Modelling and Software
30
,
47
56
.
http://dx.doi.org/10.1016/j.envsoft.2011.11.007
.
Gernaey
K. V.
,
van Loosdrecht
M. C. M.
,
Henze
M.
,
Lind
M.
&
Jorgensen
S. B.
2004
Activated sludge wastewater treatment plant modelling and simulation: state of the art
.
Environmental Modelling and Software
19
(
9
),
763
783
.
Gernaey
K. V.
,
Jeppsson
U.
,
Vanrolleghem
P. A.
&
Copp
J. B.
2014
Benchmarking of Control Strategies for Wastewater Treatment Plants, Scientific and Technical Report Series
.
IWA Publishing Company
,
London
,
UK
.
HACH Lange APS
2013
Amtax sc, Amtax indoor sc Brugsanvisning
, 8th edn.
HACH Lange APS
.
HACH Lange APS
2014
Nitratax sc – Betjeningsvejledninger
, 6th edn.
HACH Lange APS
.
Halvgaard
R. F.
,
Vezzaro
L.
,
Grum
M.
,
Munk-Nielsen
T.
,
Tychsen
P.
&
Madsen
H.
2017
Stochastic Greybox Modeling for Control of an Alternating Activated Sludge Process
.
(DTU Compute-Technical Report-2017, Vol. 08)
.
DTU Compute
.
Henze
M.
&
Comeau
Y.
2008
Biological Wastewater Treatment: Principles Modelling and Design – Chapter 3, Wastewater Characterization
.
IWA Publishing
,
London
.
Henze
M.
Jr
,
Grady
C. P. L.
,
Gujer
W.
,
Marais
G. v. R.
&
Matsuo
T.
1987
Activated Sludge Model no. 1
.
Technical report
,
IAWPRC
.
Henze
M.
,
Gujer
W.
,
Mino
T.
&
van Loosdrecht
M. C. M.
2000
Activated Sludge Models: ASM1, ASM2, ASM2d and ASM3
.
Scientific and Technical Report no. 9
.
IWA Publishing
,
London
,
UK
.
Hindmarsch
A. C.
1983
ODEPACK, A Systematized Collection of ODE Solvers
, 1st edn.
Scientific Computing (IMACS Transaction on Scientific Computation)
,
Amsterdam
.
Hu
Z.
,
Wentzel
M. C.
&
Ekama
G. A.
2007
A general kinetic model for biological nutrient removal activated sludge systems: model development
.
Biotechnology and Bioengineering
98
(
6
),
1242
1258
.
Isaacs
S. H.
&
Thornberg
D.
1998
A comparison between model and rule based control of a periodic activated sludge process
.
Water Science and Technology
37
(
12
),
343
351
.
Jazwinski
A. H.
1970
Stochastic Processes and Filtering Theory
.
Dover Publications, Inc., Mineola, New York
.
Juhl
R.
,
Møller
J. K.
,
Jørgensen
J. B.
,
Madsen
H.
2016
Modeling and prediction using stochastic differential equations
. In:
Prediction Methods for Blood Glucose Concentration. Lecture Notes in Bioengineering
(
Kirchsteiger
H.
,
Jørgensen
J.
,
Renard
E.
&
del Re
L.
, eds).
Springer
,
Cham
.
Kim
H.
,
Kim
Y.
,
Kim
M.
,
Piao
W.
,
Gee
J.
&
Kim
C.
2014
Performance evaluation of a full-scale advanced phase isolation ditch process by using real-time control strategies
.
Korean J. Chem. Eng.
31
(
4
),
611
618
.
Kristensen
N. R.
&
Madsen
H.
2003
Continuous Time Stochastic Modelling Mathematics Guide
.
Kristensen
N. R.
,
Madsen
H.
&
Jørgensen
S. B.
2004
Parameter estimation in stochastic grey-box models
.
Automatica
40
(
2
),
225
237
.
Lewis
W.
&
Whitman
W.
1924
Principles of gas absorption
.
Industrial and Engineering Chemistry
16
(
12
),
1215
1220
.
Löwe
R.
,
Vezzaro
L.
,
Mikkelsen
P. S.
,
Grum
M.
&
Madsen
H.
2016
Probabilistic runoff volume forecasting in risk-based optimization for RTC of urban drainage systems
.
Environmental Modelling and Software
80
,
143
158
.
Meijer
S. C. F.
2004
Theoretical and Practical Aspects of Modelling Activated Sludge Processes
.
PhD Thesis
,
Delft University of Technology
,
The Netherlands
.
Mulas
M.
,
Tronci
S.
&
Baratti
R.
2007
Development of a 4-Measureable States Activated Sludge Process Model Deduced from the ASM1
. In:
IFAC Proceedings vol. 1, June 4–6, Cancún, Mexico
.
Nielsen
M. K.
&
Önnerth
T. B.
1995
Improvement of a recirculating plant by introducing STAR control
.
Water Science and Technology
31
(
2
),
171
180
.
Petersen
B.
,
Vanrolleghem
P. A.
,
Gernaey
K.
&
Henze
M.
2002
Evaluation of an ASM1 model calibration procedure on a municipal-industrial wastewater treatment plant
.
Journal of Hydroinformatics
4
(
1
),
15
38
.
Rieger
L.
,
Alex
J.
,
Winkler
S.
,
Boehler
M.
,
Thomann
M.
&
Siegrist
H.
2003
Progress in sensor technology – progress in process control? Part 1: sensor property investigation and classification
.
Water Science and Technology
47
(
2
),
103
112
.
Sharifi
S.
,
Murthy
S.
,
Takács
I.
&
Massoudieh
A.
2014
Probabilistic parameter estimation of activated sludge processes using Markov Chain Monte Carlo
.
Water Research
50
,
254
266
.
Sin
G.
,
De Pauw
D. J. W.
,
Weijers
S.
&
Vanrolleghem
P.
2008
An efficient approach to automate the manual trial and error calibration of activated sludge models
.
Biotechnology and Bioengineering
100
(
3
),
516
528
.
Smets
I. Y.
,
Haegebaert
J. V.
,
Carrette
R.
&
Van Impe
J. F.
2003
Linearization of the activated sludge model ASM1 for fast and reliable predictions
.
Water Research
37
,
1831
1851
.
Sørensen
J.
,
Thornberg
D. E.
&
Nielsen
M. K.
1994
Optimization of a nitrogen-removing biological wastewater treatment plant using on-line measurements
.
Water Environment Research
66
(
3
),
236
242
.
Speelpenning
B.
1980
Compiling Fast Partial Derivatives of Functions Given by Algorithms. Technical Report UILU-ENG 80 1702, University of Illinois, Urbana, USA. doi: 10.2172/5254402
.
Stentoft
P. A.
,
Munk-Nielsen
T.
,
Mikkelsen
P. S.
&
Madsen
H.
2017
A Stochastic Method to Manage Delay and Missing Values for In-Situ Sensors in an Alternating Activated Sludge
. In:
Proceedings of ICA 2017 11-14 June 2017 in Québec City, Québec, Canada
.
Thordarson
F.
,
Breinholt
A.
,
Moller
J. K.
,
Mikkelsen
P. S.
,
Grum
M.
&
Madsen
H.
2012
Evaluation of probabilistic flow predictions in sewer systems using grey box models and a skill score criterion
.
Stochastic Environmental Research and Risk Assessment
26
(
8
),
1151
1166
.
Vanrolleghem
P. A.
,
Van Daele
M.
&
Dochain
D.
1995
Structural identifiability of biokinetic models of activated-sludge respiration
.
Water Research
29
(
11
),
2571
2578
.
https://doi.org/10.1016/0043-1354(95)00106-U
.
Wentzel
M. C.
,
Ekama
G. A.
&
Marais
G. v. R.
1992
Process and modelling of nitrification denitrification biological excess phosphorus removal systems – a review
.
Water Science and Technology
25
(
6
),
59
82
.
Zhao
H.
,
Freed
A. J.
,
Dimassimo
R. W.
,
Hong
S. N.
,
Bundgaard
E.
&
Thomsen
H. A.
2004
Demonstration of Phase Length Control of BioDenipho Process Using On-Line Ammonia and Nitrate Analyzers at Three Full-Scale Wastewater Treatment Plants
.
WEFTEC, Alexandria, VA, USA
.

Supplementary data