Implementing an Extended Kalman Filter for estimating nutrient composition in a sequential batch MBBR pilot plant

There has been a rapid increase in the implementation of advanced control strategies in process industries, including water resource recovery facilities (WRRF). These control strategies are essential for the optimal operation of process plants (O'Brien et al. 2011). Advanced control strategies such as model predictive control (MPC) would require continuous, real-time information of various wastewater compositions (Liukkonen et al. 2013).

A number of online instruments for measuring nutrient composition are available in the market today. However, their use is often limited to large-scale urban treatment facilities (Häck & Wiese 2006). Despite considerable developments in online instrumentation in the past decade, physical sensors for real-time measurement of some of the wastewater quality parameters are either extremely expensive or do not exist. Moreover, sensors available in today's market are vulnerable to fouling, drift or other inaccuracies resulting in unreliable measurement (Olsson & Jeppsson 2006). A lot of effort is put into fault detecting systems (either manual or real-time) to remove errors and correct sensors readings. Most online sensors (both ion selective and optical) require a lot of attention in terms of data quality checks, through regular inspection and adjustment (Olsson 2012). On the other hand, composition analyzers provide an accurate measurement of parameters such as ammonia, nitrate, chemical oxygen demand (COD), and phosphates, but are extremely expensive. Due to these problems, most WRRF depend on standardized laboratory measurements rather than online monitoring. However, parameters such as pH, oxidation-reduction potential (ORP), dissolved oxygen (DO) and conductivity can be measured with relative ease and higher accuracy. Therefore, these sensors, which are often inexpensive, are installed in almost all treatment plants including small to medium WRRF.

Soft sensors are an alternative to expensive composition analyzers and unreliable in-situ sensors (Thürlimann et al. 2018). Soft sensors, also denoted as state estimators, virtual sensors, etc., use system knowledge enclosed in the process model together with available measurements. These software sensors can be used not only as alternative to expensive analyzers but also as a support system for the existing composition analyzers by acting as a reliable backup in case of an instrument malfunction (Haimi et al. 2013). A detailed compilation of statistical models such as artificial neural networks (ANN), partial least-square regression (PLSR) or auto-regressive models (ARX) used for real-time estimation is provided in Haimi et al. (2013). Mechanistic models describing biological processes such as activated sludge model (ASM) can also be used with state estimators such as Extended Kalman Filter (EKF) or Moving Horizon Estimator (MHE) for developing soft sensors (Busch et al. 2013). However, most of these state estimators are tested and validated on the BSM 1 simulator. We have not found reports of implementation of state estimators with mechanistic models for soft sensing in a real system.

The aim of this work is to develop a dynamic state-estimator to provide real-time estimations of PO₄-P, volatile fatty acids (VFA), NH₄-N, and NO₃-N in a sequential batch reactor, by using data obtained from online sensors such as pH and DO.

The sequential batch moving bed Bio-P (SB-MBBR Bio-P) pilot plant, located in the wastewater laboratory at Norwegian University of Science and Technology (NTNU) was used to test and validate the soft-sensor algorithm. Figure 1 depicts the pilot plant and the data acquisition system. The sequential batch reactor (SBR) contains polyurethane carriers (with 60% filling) which are exposed to raw wastewater to ensure biofilm growth and attachment to the carriers (see Supplementary Figure 1, available with the online version of this paper). The reactor has a total working volume of 13 L. It operates with a fixed cycle time of 8 h and a temperature of 16 ± 0.5 °C. The carrier media have a cylindrical shape (10 mm in diameter and length) with a cross inside the cylinder and longitudinal fins outside, providing a large surface area for biofilm growth (Ødegaard et al. 1994). The SBR is fed with wastewater from a storage tank, which receives steady supply of raw wastewater from a nearby municipal sewer. The storage tank has a total volume of 3.5 m³ and a residence time of 24 h, which ensures an uninterrupted supply of wastewater while maintaining sufficient variations in influent wastewater quality.

After ensuring the attachment of active biofilm to the carrier, the reactor operates in an alternating anaerobic–aerobic cyclic process with a hydraulic retention time of 8 h (3 h anaerobic, 5 h aerobic, 3 min filling and 3 min emptying). A sieve at the outlet line retains the carriers in the SBR, allowing just the treated water and detached biomass to exit the reactor. The reactor operates as a biofilm process and all active biomass is attached to the carrier media. Therefore, the SBR can be operated with a 100% volume exchange ratio, resulting in a cycle time equal to the hydraulic retention time. The system is mechanically agitated using a custom-made mixer for the entire cycle period. During the aerobic phase, air is supplied with a constant flow rate of 4.5 L min⁻¹. The algorithm to execute the cycle time and the operational sequence is stored in the Supervisory Control and Data Acquisition (SCADA) system of the plant.

During the anaerobic-stage, the active biomass (primarily PAOs) utilizes the acetate from raw wastewater and converts it into polyhydroxy-alkanoates (X_PHA). As a result, PO₄-P is released due to Poly-P degradation, thus increasing its concentration in bulk solution. In anaerobic conditions, the readily biodegradable organic matter (BOD) is also converted to VFA by fermentation, and the PAOs utilize it for further PO₄-P release. When aeration begins, the PAOs utilize the stored X_PHA as an energy source to grow and take up PO₄-P from the water phase and synthesize new storage products within the biomass (X_PP). The biomass also contains significant proportions of nitrifying organisms, which convert NH₄-N to NO₃-N. The biomass also contains heterotrophs, which consume BOD during the aerobic stage.

When the biofilm layer on the carriers reach a steady thickness, the system exhibits more than 80% PO₄-P removal. In this period, a constant amount of PO₄-P is added at the beginning of every cycle to avoid low influent PO₄-P concentration due to dilution of wastewater by storm-water (9 ± 3 mg P L⁻¹). However, the soluble COD (127 ± 30 mg COD L⁻¹), NH₄-N (28 ± 16 mg N L⁻¹) is not adjusted.

After the biomass activity shows steady conditions, kinetic tests are conducted in the reactor. Samples are collected for the entire cycle every 30 min and immediately filtered through 0.45 μm filters to determine the concentration of soluble COD (sCOD), orthophosphate (PO₄-P), ammonium nitrogen (NH₄-N) and nitrate nitrogen (NO₃-N). All these compositions are measured using Dr. Lange cuvette and HACH LANGE GmbH (DR 1900, China) spectrophotometer. All analyses are performed immediately during the cycle test. The DO and pH sensor installed in the pilot plant are connected to the plants SCADA system. The data from the online sensors are logged every minute. The data received from the laboratory analyses and the online sensors are used to calibrate the model.

A number of mechanistic models explaining the process of removal of phosphorus, nitrogen, and carbon in sequential batch reactors are available in the literature. However, using a complex model such as ASM2d (Henze et al. 1999) for estimation and control would be challenging, especially when the number of measurements available in the reactors are limited (Jeppsson & Olsson 1994). The usual approach is to reduce the number of state variables by either lumping two or more states or assuming some of the states with relatively slower dynamics as constants (Steffens et al. 1997). We find similar model-reduction strategies implemented on ASM1 model, for use in state estimation and control (Julien et al. 1999; Gómez-Quintero et al. 2000; García-Usach et al. 2006).

The original ASM2d model (Henze et al. 1999) consists of 19 states and 21 processes. It describes the growth and decay of three types of biomass; the phosphate accumulating organisms (X_PAO), autotrophs (X_AUT) and heterotrophs (X_H), followed by the subsequent uptake and release of various substrates and storage products in the biomass. The following assumptions, similar to the simplification strategies presented in Zhang et al. (2010), Cadet (2014) and Li et al. (2018) were made to reduce the model.

• The effect of particulate hydrolysis reactions appears to be insignificant in the process and has been ignored.

• The non-reacting states in the model: particulate inert (X_I), soluble inert (S_I) and alkalinity (S_ALK) are eliminated.

• Since there is no anoxic condition in the SBR's operational sequence, all kinetics pertaining to denitrifiers in the biomass are removed. However, the denitrification reactions by X_PAO has been retained.

• The effect of growth and decay kinetics of all three biomass types can be excluded. Hence, the model only considers the kinetics of substrate degradation processes.

This reduced order model explains the dynamics of storage and consumption of soluble substrates such as readily biodegradable substrate (S_F), volatile fatty acids (S_A), phosphates (S_PO), ammonia (S_NH) and nitrates (S_NO). The model also includes the kinetics of storage products such as polyhydroxy-alkanoates (X_PHA) and polyphosphates (X_PP). The description of substrates and storage products used in the model is presented in Table 1. The list of processes included in the reduced order ASM2d model is presented in Table 2.

The following modifications are implemented to reduce the number of parameters in the process rate kinetics.

• All the biological processes involved in the aerobic uptake of (S_F) are written as a Monod equation with a rate constant and saturation coefficients K_S and K_O.

• The maximum fermentation rate is lumped with X_H and represented as a single term⁠.

• The rate constant for X_PHA storage is lumped with X_PAO to form a term ⁠.

• The rate constant for X_PP storage is lumped with X_PAO to represent a single term ⁠.

• The maximum growth rate of autotrophs and is lumped with X_AUT to represent a single term ⁠. This rate equation assigned for aerobic growth of autotroph is used to explain the kinetics of ammonia nitrification.

• All the reactions pertaining to denitrification by PAO are lumped with X_PAO to represent a single term ⁠.

• The Monod saturation coefficients associated with X_PHA (K_PHA) and X_PP (K_PP) and (K_MAX) was lumped with X_PAO to obtain the terms K^’_PP, K^’_PHA and K^’_MAX, respectively.

After considering the aforementioned modifications, the rate kinetics can be transformed to the equations mentioned in Table 2. All the kinetic parameters mentioned in this model will be estimated in the model calibration stage.

‘p’ is a vector consisting of the model parameters that are to be estimated; the upper and lower bounds of the parameter vector is represented as ub and lb, respectively. The objective function f_OBJ is defined as the least-square error between the experimental values y_i = [sCOD NH₄-N NO₃-N PO₄-P] and its corresponding model predicted values ⁠. The half-saturation constants and stoichiometric parameters were constrained within −100% and +150% of the original ASM2d model parameters, but a wider constraints range is provided to the rate constants. The initial concentrations of S_F and X_PP are also taken as parameters, which would be estimated by the optimizer. The non-linear optimization solver function in MATLAB (fmincon) is used to solve the optimization problem and obtain the parameters.

The discrete form of state-space equation presented in Equation (8) is derived from the continuous-time system function f_c (in Equation (1)) by an explicit Euler discretization method with a constant sampling time of T_S. f_c and h are the nonlinear state transition and measurement functions. In the equations above, x_k is the state variable at time k, z_k is the measurement, Q and R are the covariance matrices of the process and measurement noise, respectively. F is the state transition matrix, and H is the measurement matrix. I is the identity matrix, is the a priori estimate of the state, the a posteriori estimate of the state, the Kalman gain, the covariance of a priori estimation error, and the covariance matrix of the a posteriori estimation error.

Before implementing the estimator to the pilot plant a series of simulations are conducted to identify the appropriate tuning parameters that would result in a reasonably good correction of the states by the EKF. For this purpose, the mathematical model and the EKF equation described in the previous section were implemented in MATLAB. The model equations are written symbolically and the Jacobian function in MATLAB's symbolic toolbox is used to linearize the process model and obtain the F and H matrix at every time step. The implementation schematic is presented in Figure 2. The block mentioned as ‘simulator’ consists of the mathematical model and ‘plant’ represents the real pilot plant. The EKF equations are written as a separate section in MATLAB, which can accept the measurement z_k and the control variable u_k from either a simulator or from the plant's SCADA system. This shift from simulator based testing to a pilot plant implementation is done with the switch S1 and S2, as indicated in Figure 2. During the tuning of the EKF parameters, the simulated state vector, x_k generated by the model is compared to the estimated state predicted by the EKF to test the estimator's ability to provide a reliable estimation.

The SBR pilot plant SCADA provides an Open Platform Communication (OPC) server, which enables communication with other devices through the standard OPC protocol. A communication layer is written in MATLAB using its OPC toolbox to pull real-time data from the OPC server into the MATLAB workspace. This method enables scripts written in MATLAB to be used directly in any industrial PLC without the hassle of re-writing the code in a different programming language. While implementing the EKF in the pilot plant, a few overriding constraints are included in the estimator algorithm to prevent EKF from overshooting during the emptying and filling stages of the SBR. The states and the auto-covariance matrix are reinitialized at the beginning of every cycle and the estimator is kept on hold until the anaerobic conditions begin (DO = 0) and the pH sensor value stabilizes. Once the EKF is started, the algorithm continues to run until the end of the aerobic stage. When the emptying stage begins, the estimated states are frozen until the next filling stage begins.

The results of the model calibration step, obtained by the proposed optimization procedure, are presented in Table 3. The graphs presented in Figure 3 show the comparison between the model predicted values and the data obtained from the kinetic studies conducted in the pilot plant. The figure also provides the normalized root-mean-square error (NRSME) and R² values obtained during the model calibration. The R² values presented in the graph show that the values predicted by the model are close to the measured data collected from the pilot plant. Hence, the model can be considered to provide a good representation of the SBR process.

The optimization algorithm used for estimating the model parameters (Equations (2)–(4)) is also included in the MATLAB code. The model parameters presented in Table 3 can be updated by running the optimizer code with a new data set obtained from laboratory analysis. This strategy provides and effective platform for model calibration and effortless update of parameters whenever a drift is observed between the laboratory results and estimated values. It also generates a possibility of adapting model parameters such that the soft-sensor algorithm can be implemented in any SBR plant with similar operational sequence.

The values of⁠, ⁠, Q and R described in Equations (8)–(15) can be used as the tuning parameters of the Kalman Filter. In practical applications, these parameters are not known precisely, and therefore, trial-and-error tuning must be expected. We have tuned the EKF as explained below, mainly following the guidelines in Haugen et al. (2014).

Tuning of⁠: It is commonly set as a diagonal matrix. The values are adjusted based on how close our initial estimates are to the real value of states. The knowledge on raw wastewater composition can help provide a good guess for the initial composition of some of the states in the process.

• The raw wastewater does not contain any nitrates, therefore the value of S_NO is 0 g-N m⁻³.

• The stored PHA in the biomass usually depletes before the cycle ends, therefore we can assume that the initial concentration of X_PHA0 = 0 mg COD L⁻¹.

• S_O = 0 mg L⁻¹ since there is no dissolved oxygen in the anaerobic stage.

• The states S_F and X_PP do not affect the measurements (pH) in the anaerobic zone. These states exist in the model merely for the purpose of maintaining the model dynamics. A fixed initial value of X_PP = 65 mg P L⁻¹ and S_F = 45 mg COD L⁻¹ is provided to make the model fit to the kinetic data. These values should not be confused with the actual P content in the biomass or the available soluble BOD.

Tuning of R: It is estimated as the variance of the measurement time series collected from the SCADA system.

The EKF is updated with the tuning parameters presented in Equations (16)–(20). The results presented in Figure 4(a) and 4(b) show a comparison between the state variable obtained from the simulator and the states estimated by the tuned EKF during an 8 h cycle. The state correction by EKF is not clearly visible in Figure 4(a) and 4(b) since the estimated states usually approach the simulated states within the first 30 min of the SBR cycle. Therefore, the simulation result for the first 30 min is presented in Figure 4(c), where the state correction by the EKF is noticeable. The convergence of the estimated value to a randomly selected initial composition confirms that the EKF parameters are correctly tuned.

The tuned EKF algorithm is implemented in the pilot plant and the estimator results are recorded for a period of two weeks. Figure 5 presents a subset consisting of estimated states for a duration of 3 days (nine cycles). The variations in nutrient composition during the cycle follow similar trends to those observed in a typical SBR operation, as presented in literature such as Marsili-Libelli et al. (2001) and Sin & Vanrolleghem (2006). Due to the high holdup in the buffer tank, the diurnal variations in influent concentration is barely noticeable. However, a significant change in influent concentration can be observed during rain or snowmelt events, causing storm-water infiltration at the raw wastewater source. These variations can be observed in Figure 5(e), which shows the data from the online turbidity sensor installed in the raw wastewater tank. The estimation results presented in Figure 5(a)–5(d) show the increase in the estimated value of S_NH, S_PO, and S_A on 31 October 2018 around 12.00 h. This event corresponds to the influx of wastewater with a rather higher nutrient concentration at 9:00 AM on the same day. This provides a qualitative validation for the soft-sensor's ability to detect variations in influent wastewater quality (see Supplementary Table 2, available with the online version of this paper).

Further validation tests were conducted on 4 November 2018. Samples were collected every 30 min during the complete 8-h cycle and the composition of NH₄-N, PO₄-P, and NO₃-N are measured using standardized laboratory tests. These values are compared to the states S_NH, S_PO, and S_NO estimated by the EKF algorithm during the cycle. The comparison plots along with the NRSME and R² values are presented in Figure 6. This demonstrates the EKF's potential to estimate the states in the pilot plant.

The estimated value of the effluent water quality is also recorded for each treatment cycle. This is executed by logging the estimated value of S_NH, S_PO and S_NO at the end of every SBR cycle. The pilot plant achieves complete nitrification, resulting in a very low effluent ammonia concentration close to zero. However, the effluent concentrations of S_NO and S_PO are rather susceptible to variations in raw wastewater quality. These variations can be observed in Figure 7(a) and 7(b).

The soft-sensor is also used to estimate the ammonia concentration in the raw wastewater. It is observed that the concentration of NH₄-N stays nearly constant during the anaerobic stage. However, it can also be observed from Figure 4(b) that the estimator converges to the real state variables within the first 30 min of the cycle time. Therefore, it would be safe to assume the value of S_NH at t = 30 min to be very close to the concentration of ammonia in raw wastewater. The estimated value of S_NH at t = 30 min for the entire testing period is presented in the third plot in Figure 7(c). These estimated values of the influent and the effluent wastewater composition is compared to the periodic laboratory measurements. Except for the erroneous estimates received on 3-5 November 2018 due to an unexpected power failure, the plots presented in Figure 7(c) demonstrate an acceptable match between the experimental data and the values estimated by the EKF algorithm.

The grey box model discussed in this article explains the SBR with sufficient level of accuracy to constitute the basis of a state estimator for the reactor. The model calibration strategy discussed in this work provides the flexibility in adapting the model to any plant with similar process configuration. The validation tests conducted in the pilot plant demonstrates the possibility to estimate the concentration of VFA, NO₃-N, NH₄-N and PO₄-P in a SBR operating in alternative anaerobic–aerobic mode for the biological phosphorus removal process, by using physical sensors such as pH and DO. The EKF based soft-sensor can also provide accurate estimations of the ammonia concentration in the influent wastewater. Since the EKF algorithm is written in MATLAB's OPC compatible framework this soft-sensor code can easily be implemented in any treatment plant which has an OPC enabled PLC, SCADA or a distributed control system (DCS), without additional hardware requirement or modifications to the MATLAB code.

The change in biomass concentration over a longer operational period could lead to minor drift in state estimation. However, periodic model calibration and parameter re-tuning can resolve this problem. The optimization code estimating the model parameters is also included as a part of the overall estimator package. This allows the model to be recalibrated with relative ease. A cost-effective method of state-estimation using soft sensors can be a prelude to implementing advanced model based control strategies in wastewater treatment and recovery processes.