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

Online monitoring of water quality parameters can provide better control over various operations in wastewater treatment plants. However, a lack of physical online sensors, the high price of the available online water-quality analyzers, and the need for regular maintenance and calibration prevent frequent use of online monitoring. Soft-sensors are viable alternatives, with advantages in terms of price and ﬂ exibility in operation. As an example, this work presents the development, tuning, implementation, and validation of an Extended Kalman Filter (EKF) on a grey-box model to estimate the concentration of volatile fatty acids (VFA), soluble phosphates (PO 4 -P), ammonia nitrogen (NH 4 -N) and nitrate nitrogen (NO 3 -N) using simple and inexpensive sensors such as pH and dissolved oxygen (DO). The EKF is implemented in a sequential batch moving bed bio ﬁ lm reactor (MBBR) pilot scale unit used for biological phosphorus removal from municipal wastewater. The grey-box model, used for soft sensing, was constructed by ﬁ tting the kinetic data from the pilot plant to a reduced order version of ASM2d model. The EKF is successfully validated against the standard laboratory measurements, which con ﬁ rms its ability to estimate various states during the continuous operation of the pilot plant.


INTRODUCTION
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. ). Advanced control strategies such as model predictive control (MPC) would require continuous, real-time information of various wastewater compositions (Liukkonen et al. ).
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 ). 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 ). 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 ). 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. ). 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. ). 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. (). 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. ). 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 stateestimator to provide real-time estimations of PO 4 -P, volatile fatty acids (VFA), NH 4 -N, and NO 3 -N in a sequential batch reactor, by using data obtained from online sensors such as pH and DO.

Pilot plant setup
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. ). 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 3 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 anaerobicaerobic 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 À1 . 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 4 -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 4 -P release. When aeration begins, the PAOs utilize the stored X PHA as an energy source to grow and take up PO 4 -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 4 -N to NO 3 -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 4 -P removal. In this period, a constant amount of PO 4 -P is added at the beginning of every cycle to avoid low influent PO 4 -P concentration due to dilution of wastewater by storm-water (9 ± 3 mg P L À1 ). However, the soluble COD (127 ± 30 mg COD L À1 ), NH 4 -N (28 ± 16 mg N L À1 ) 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 4 -P), ammonium nitrogen (NH 4 -N) and nitrate nitrogen (NO 3 -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.

Mathematical 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. ) for estimation and control would be challenging, especially when the number of measurements available in the reactors are limited ( Jeppsson & Olsson ). 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. ). We find similar model-reduction strategies implemented on ASM1 model, for use in state estimation and control (Julien et

Model reduction strategy
The original ASM2d model (Henze et al. ) 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. (), Cadet () and Li et al. () 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 r 1 and saturation coefficients K S and K O .
• The maximum fermentation rate q fe is lumped with X H and represented as a single term r 2 .
• The rate constant for X PHA storage q PHA is lumped with X PAO to form a term r 3 .
• The rate constant for X PP storage q PP is lumped with X PAO to represent a single term r 4 .
• The maximum growth rate of autotrophs μ AUT and Y A is lumped with X AUT to represent a single term r 5 . 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 r 6 .
• 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.

State-space equations
The state-space equation denotes the rate of changes in the state variables. Equation (1) are the state variables and u ¼ K L a is the manipulated variable. The process rate equations f c are presented as a Petersen matrix in Table 2.
The reduced order ASM2d model has to be calibrated before use in a specific case (García-Usach et al.
'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 leastsquare error between the experimental values y i ¼ [sCOD NH 4 -N NO 3 -N PO 4 -P] and its corresponding model pre- The halfsaturation 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.

Observation equations
The observation equation expresses the online measurements as functions of the state variables (see Supplementary  Table 1, available with the online version of this paper). The work presented in Serralta et al. () and Aguado et al. () demonstrates the possibility of expressing the pH as a function of various ionic species available in the wastewater treatment system (S A , S PO , S NH , S NO ). A generic equation, which contains a weighted combination of the linear terms, squared terms, and binary interaction terms, is used as the preliminary model structure. The mathematical form of the equation is presented in Equations (5) and (6). A similar model structure can be found in dosing control predictors, and are found to fit a wide range of data (Ratnaweera & Fettig ). The regression learner app in MATLAB is used to establish a correlation between the ionic species (θ ¼ [S A S NH S NO S PO ]) and the data received from pH sensors and obtain the coefficients α i β i and γ i,j . The second measurement, dissolved oxygen is equivalent to the state S O in the model. The measurement function h(x) can be expressed as Equation (7).

State estimator
The Kalman Filter is an optimal state estimator for linear systems with Gaussian noise. However, a number of versions of original Kalman Filter algorithm such as Unscented Kalman Filter (UKF), and particle filters (PF) can be used to estimate the states of a non-linear system. In this case, the Extended Kalman Filter is used as a state estimator for the SBR. The equations of the EKF are given in Equations (8)- (15). The equations were written in the same order in which they were implemented in the PLC.
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, x À k is the a priori estimate of the state, x þ k the a posteriori estimate of the state, K k the Kalman gain, P À k the covariance of a priori estimation error, and P þ k the covariance matrix of the a posteriori estimation error.

Implementing the EKF
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 x þ k 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.

Model calibration results
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-meansquare error (NRSME) and R 2 values obtained during the model calibration. The R 2 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.

Tuning parameters of the EKF
The values of x þ 0 , P þ 0 , 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. ().
Tuning of x þ 0 : A reasonable value of x þ 0 is a representative steady state, here set equal to the designed operating point of the process. For an SBR operation, we chose the design values of influent water quality as x þ 0 .
Tuning of P þ 0 : 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 À3 .
• 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 À1 .
• S O ¼ 0 mg L À1 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 À1 and S F ¼ 45 mg COD L À1 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.
Therefore, the only uncertain initial estimates are S A0, S NH0 , and S PO0 . The P þ 0 corresponding to the states S F0 , X PP0 , X PHA0 , S NO0 , and S O0 can be set to a very low value or to zero. Therefore, the P matrix at the beginning of the cycle can be expressed as the following. Downloaded from https://iwaponline.com/wst/article-pdf/80/2/317/621359/wst080020317.pdf" /><meta name="description" content="Abstract. Online monitoring of water quality parameters can provide better control over various o by guest  By trial and error, it was found that the following values provide a good estimation. P S A ¼ 316 ; P SNH ¼ 0:12; P S PO ¼ 0:10 (18) Tuning of R: It is estimated as the variance of the measurement time series collected from the SCADA system.
Tuning of Q: Typically, Q is used as the main, or final, tuning parameter. It is assumed a diagonal matrix with diagonal elements related to the initial guess of the pertinent state variable through factors l i : In the simulator based testing, proper values of trial-anderror were found as {l i } ¼ [1 1 1 1 2 1 1 1:5] Ã 10 À3 (20)

Simulator test results
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.

Pilot plant implementation results
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. () and Sin & Vanrolleghem (). 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 4 -N, PO 4 -P, and NO 3 -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 2 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 4 -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.

CONCLUSIONS
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 3 -N, NH 4 -N and PO 4 -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 softsensor 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 retuning 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.