Performance Prediction and MgCl Setpoint Adjustment of a Utility-Scale Struvite Precipitation Reactor

Define Problem:

Metro Water Recovery (MWR) in Denver, CO uses enhanced biological phosphorus removal (EBPR) to remove phosphorus from their main treatment stream. During subsequent anaerobic digestion (AD) of the EBPR sludge, a portion of the biomass-bound phosphorus is re-solubilized and released into the sidestream of the plant. This high concentration of phosphorus leads to inadvertent struvite precipitation that causes scaling in pipes and fixtures (Figure 1) resulting in operational shutdowns and increased maintenance costs. Additionally, after recycling of the sidestream into the main stream, the phosphorus loading of the plant actually increases over time.

Figure 1. Images of struvite scaling on pipes at MWR’s facility

To address these issues and help meet effluent discharge limits, MWR uses a controlled struvite precipitation reactor. These reactors transform soluble orthophosphate (OP) in the anaerobic digester effluent into struvite (NH4MgPO4 6H2O) precipitate, enabling the removal and recovery of phosphorus as well as protecting critical infrastructure from scaling. In 2020, MWR installed a commercial struvite precipitation reactor to treat their AD effluent prior to sludge dewatering and sidestream recycling. The key control variables for the system are the MgCl dosage and the rate of aeration for pH control (Figure 2).

Figure 2. Diagram of the MagPrex® system installed at MWR. The key control variables are boxed in red.

Currently, MWR is adjusting their MgCl dosage based on a 1.3:1 Mg-P ratio with weekly orthophosphate (OP) measurements from the AD digestate. Due to the extended time it takes to obtain the OP laboratory measurements, MWR is unable to get an accurate, real-time understanding of the reactor’s performance and adjust set points accordingly. Furthermore, laboratory experiments can be expensive and time consuming. Costs, maintenance requirements, and inaccuracies in measuring OP concentrations in the complex digestate matrix have prevented adoption of commercially available sensors and online analyzers that could help with optimizing control processes. An investigation into leveraging a variety of statistical and machine learning (ML) methods is therefore undertaken to provide daily predictions for the influent OP, effluent OP, and OP removed from the reactor and to recommend MgCl setpoints (MgCl SP) that would lead to MWR’s ideal reactor effluent OP concentration of 30 mg P/L. 

Go/No-Go Decision: If your problem can be accurately modeled with a mechanistic model, improved by operator staff training, or if existing data collection infrastructure is not in place then data-driven modeling may not be the right solution.

Over 15 months, MWR collected a dataset spanning 95 variables and 439 days. 65 variables were laboratory measurements (e.g., digester effluent magnesium and OP concentrations) with observations generally occurring 1-3 times a week. 30 variables were online measurements (e.g., pH and MgCl feed) with a time weighted, average observation occurring every day. All of the laboratory and online measurements were pulled from MWR’s internal data collection system into separate MS Excel files. These Excel files were then imported into Python, merged together, and matched up by timestamp.

Prior to moving forward with prepping the data and building the models, it was important to start with some exploratory data analysis (EDA) to better understand the dataset that was collected. EDA typically involves plotting histograms, scatterplots (Figure 3) , and boxplots of individual variables in a dataset to help quickly identify diagnostic information and plan out next steps. We plotted key variables of interest to make sure we had enough variation in our data to justify modeling them. Of these key variables, we found that the aeration blower used to control pH was turned off for the majority of the data collection timeline. While initially we were interested in modeling a blower set point model, due to this lack of data, we instead focused on developing influent OP, effluent OP, and OP removal models.

Figure 3. Initial scatter plots of key operational data.

Go/No-Go Decision: ML may not be a viable solution to your problem if your dataset contains a majority of missing or constant data. If your dataset contains data with little variance then your ML models may not perform any better than taking the statistical mean of the data.

Prepare Data:

An important step in any data-driven modeling process, is to clean and process the raw dataset. Laboratory measurements are typically taken at various timestamps throughout the day, so we resampled the dataset to combine all observations within the same day into a single row of daily observations. We then removed any row in the dataset that did not include values for digester effluent OP (influent OP into the reactor), centrifuge feed OP (effluent OP out of the reactor), and MgCl discharge (indicator that the system was online) to ensure a baseline amount of knowledge on the system’s performance was included. Additional variables were removed that had the vast majority of their observations missing or constant.

After this initial data cleaning, we focused attention towards outlier removal and data imputation. Through EDA and consultation with MWR’s process engineers, we began by manually identifying known outliers in the dataset. We then used kernel density estimation (KDE) with a gaussian kernel to fit probability density curves to each individual variable’s data (Figure 5). We iteratively adjusted the bandwidth and significance level (cutoff for what is considered an outlier or not) of the KDE until the known outlier points were identified. Resultantly, an observation from each feature was identified as an outlier when it had a log probability less than the 10th percentile. All manually and statistically identified outlier observations were then set equal to N/A.

Due to the outlier removal process and the infrequency of laboratory measurements, several rows of timestamps included missing or N/A observations. Imputation was thus performed to fill in unavailable data. The last-observation-carried-forward technique, which assumes data on days with missing observations is equal to the value of the last measured observation, was used first. For any remaining missing observations, data was backfilled with the first available measurement. Once imputation was finished, OP removal values which equal the difference between the influent OP and the effluent OP concentrations were added to the dataset as an additional variable and any remaining outliers were removed using the IQR method. Our final dataset consisted of 87 variables and 152 observations.

Go/No-Go Decision: ML modeling performs best with a large amount of observations. If your dataset has more variables than observations or a small set of observations that do not capture the variability of your process then you may not want to proceed.

Six modeling methods were evaluated on predicting three different objective variables: influent OP, effluent OP, and OP removal. To prevent data leakage by giving the models extra information that could cause over-optimistic predictions, the objective variables and other OP-related data were dropped from the independent variable set. For this write up, the influent and effluent OP predictions are used as examples. To see results from OP removal predictions, visit the main project website.

Before moving on to more computationally complex ML modeling methods, it is often good practice to evaluate whether simpler and more easily interpretable statistical models can be effective at meeting your objective. Four different types of statistical regression modeling methods including linear regression (LR), partial least squares (PLS), ridge, and lasso regression were therefore evaluated.

To see if ML is a more effective modeling tool for this case, eXtreme Gradient Boosting (XGBoost) was evaluated as well. A baseline model utilizing the mean values of the OP removal and the MgCl SP data for predictions was also developed as a point of comparison for all other models. 

For the ridge, lasso, and XGBoost models there are a variety of model hyperparameters that can be tuned to improve prediction accuracies. For our ridge and lasso models, cross validation was used to identify ideal alpha values from the training data. For the XGBoost model, which has several more hyperparameters, a bayesian tuner was used instead.

To evaluate each model on its ability to predict a variety of unseen data, a repeated k-fold cross validation with five folds and one repeat was used to split up the processed dataset (Figure 6). The root mean squared error (RMSE) between the actual measured Y variable data and the predicted Y variable data from each of the cross validation splits was calculated and averaged to provide a score for all models. The standard error (SEM) of the average RMSE values was also calculated.

Given the large number of variables in our dataset compared to the number of observations, a variety of dimension reduction methods were investigated to reduce model complexity and filter noise. A comparison between three different dimension reduction techniques, Spearman’s correlation (ρ), Principal Component Analysis (PCA), and Maximum Relevance-Minimum Redundancy (mRMR), was undertaken. Additional techniques used in the full project can be found in the write-up on the main project website. 

ρ is a nonparametric variation measuring the strength and direction of monotonic relationships. The correlations of all variables to the Y variables were found using ρ and the lowest correlated variables were removed.

PCA involves creating linear combinations of existing variables to sequentially capture the most variation in the data. Across these linear combinations (i.e, principal components), those variables with the lowest contributions to the variations in the data were removed. 

For mRMR, a ranking based on F-scores and correlation coefficients was created and the highest scoring variables were kept.

To identify the best performing dimension reduction technique and ideal number of variables to keep, the six modeling approaches were trained and tested on different numbers of variables. Training occurred on the first five variables selected and then was repeated for five additional variables until 30 were selected, this process was repeated for each of our dimension reduction techniques.

Figure 6. Visual representation of a repeated k-fold cross validation process. Adapted from Sontakke et al. 2019. Classification of Cardiotocography Signals Using Machine Learning: Proceedings of the 2018 Intelligent Systems Conference.

The results from each dimension reduction technique and model evaluation were averaged and plotted for comparison (Figure 7).

For the influent OP predictions, selecting 20 variables with the mRMR method and using the linear regression model resulted in a minimum RMSE of 22.4 mg-P/L representing a 32.3% improvement over the mean model. For the effluent OP predictions, selecting 30 variables with PCA and using the XGBoost model resulted in a minimum RMSE of 12.2 mg-P/L representing a 53.6% improvement over the mean model.

These improvements over the mean help to justify the use of data driven modeling for predicting the objective variables and indicate MWR can reduce the number of monitored variables while maintaining similar or lower errors. Across the dimension reduction techniques and when considering all modeling types, there is little difference in RMSE values after the selection of 5 or more variables. This suggests a few key variables are influencing the majority of the predictive abilities of the models and that the type of dimension reduction technique used will be most influential when selecting a relatively small number of variables to keep.

Figure 8. Plots of the predicted values versus the actual measured values. The left plot shows the Influent OP model results using mRMR dimension reduction and a linear regression model. The right plot shows the Effluent OP model results using PCA dimension reduction and an XGBoost model.

Deployment involved operationalizing the model code and integrating it with MWR’s internal data historian system. Python scripts were created to pull new data every day from the historian, integrate it with the historical raw data, and process it following the same procedures used during development and training. Code for predicting the three modeling objectives and MgCl SP recommendations were also created. MgCl SP recommendations were obtained using the predicted OP concentration data multiplied by MWR’s 1.3:1 Mg-P ratio. To sufficiently capture changes in data while conserving computing power, it was determined that model retraining on any new data would occur every two weeks. Prediction outputs were set up to publish results back to the historian and be shown to MWR engineers and operators. Laboratory OP tests continued to be run to help validate the deployed model predictions.

Go/No-Go Decision: If water quality changes not captured in the original dataset lead to worsening model results or if deployed code generally performs worse than expected, consider retraining the models, reevaluating your data collection process, or changing your modeling objectives.