Economic Model Predictive Control for optimal struvite recovery

Resource recovery from municipal wastewater has been a prime focus for a decade. Although several recovery processes already exist in the market today, the high cost of material, inherent disturbance in the influent quality, lack of real time monitoring of critical parameters, and lack of a robust automation system may result in suboptimal performance. This work attempts to construct a model based predictive control for optimal operation of a struvite recovery unit in a full scale WRRF. A multi-parameter based predictive control has been developed by implementing an Economic Model Predictive Controller (EMPC) for optimal dosing of magnesium hydroxide in a struvite recovery unit. The EMPC used customized objective function for real-time optimization of performance and economical parameters of the crystallization unit. The effectiveness of the proposed EMPC controller is verified through tests conducted on the Benchmark Simulation Model No. 2 (BSM2d.). The results obtained from the simulator-based evaluation of EMPC demonstrate a significant improvement in resource recovery at reduced operational costs. The economic advantages of implementing an EMPC compared to proportional and constant magnesium dosage has also been enumerated.


Introduction
The Wastewater Treatment Plants (WWTPs) are consistently upgrading their processes to include more recovery operations and conform to its new terminology as Water Resource Recovery Facilities (WRRFs) (Regmi et al., 2019). Innovative treatment technologies are being implemented to enable better processing and disposal of the wastewater sludge. Struvite precipitation is one such process that gained popularity over the past decade (Jensen et al., 2015). Struvite (magnesium ammonium phosphate hexahydrate) has been of special interest due to its potential applicability as a slow-release fertilizer. The use of slow-release fertilizers (such as struvite) can offset the environmental deterioration caused by the excessive use of mineral-based fertilizers and eventually play a vital role in the modern eco-friendly sustainable agricultural sector (Rahman et al., 2014). Moreover, depleting reserves of mineral phosphorus also encouraged to explore alternative renewable sources (Cordell and Bennett, 2011). Struvite can be precipitated by adding Magnesium Hydroxide (Mg(OH) 2 ) to a stream rich in ammonium (NH 4 + ) and phosphate (PO 4 − ) ions. In a typical WRRF, these nutrient-rich streams can be found in the supernatants of the anaerobic-digested sludge (Rahaman et al., 2008). Therefore, a struvite recovery unit is often installed after the dewatering unit to recover the Nitrogen (N) and Phosphorus (P) before recycling the supernatant back to the biological reactors.
Although the struvite precipitation process was designed with an aim of generating a commercially marketable slow-release fertilizer, in most WRRFs it is often used as a strategy to prevent the scale formation along the sludge train and eventually reduce the maintenance cost (,b). Several social, economic, and technological reasons can be attributed to their inability to produce a commercially marketable product. The prices of Mg(OH) 2 and the energy required to operate fluidized bed reactors result in higher production costs. The fluctuations in the nutrient concentrations of the influent supernatant and the resulting inability to maintain a stable product quality add to the problem. These disturbances can often result in suboptimal performance of the crystallizers designed for struvite production. Instrument Control and Automation (ICA) offers several control strategies for ensuring optimal operation of various treatment processes. Although several works have already presented the advantages of introducing a struvite crystallization unit in a WRRF (Mbamba et al., 2016), we could not find control strategies for optimal operation of struvite production process. Several experimental studies have been conducted for the purpose of identifying the optimal operating conditions for struvite precipitation (Forrest et al., 2008;Jia et al., 2017). However, a continuous full-scale struvite production facility with ubiquitous disturbance in the influent stream benefits from real-time optimization of operational parameters. A well-designed control strategy using feedback from online nutrient sensors can reduce operating costs, improve struvite production, and enhance the performance of struvite crystallization unit.
The conventional single-input-single-output (SISO) control strategies used in the operations of WWTP have to be significantly upgraded to ensure profitable operations in a WRRF (Vanrolleghem and Vaneeckhaute, 2014). Model Predictive Control (MPC) is an advanced control strategy used for optimal control of various operations in process industries (García et al., 1989). Several implementations of MPC in wastewater treatment processes are also reported in literature (Ostace et al., 2012;O'Brien et al., 2011;Hasanlou et al., 2019). MPC decides the control movies based on a cost minimizing control strategy over a finite time-horizon using a mathematical model (mechanistic or data-driven) of the process. A variant of MPC, that uses the economic cost parameters as their objective function, called as Economic Model Predictive Control (EMPC), has also been reported in literature (Durand et al., 2016). Multi-parameter based optimal control strategies such as EMPC are especially suitable for processes that do not have a constant optimal operating setpoint (Ellis et al., 2017). Processes such as struvite recovery which is subjected to a considerably high influent disturbance in terms of flowrate, nutrient concentration as well as the final price of products, if operated under constant Mg(OH) 2 or buffer dosing could result in suboptimal performance (Crutchik et al., 2017). Strategies such as EMPC which decide the control moves based on a real-time optimization of an economic cost function can be a potential strategy for optimal control of the struvite crystallization unit.
This work aims to develop a multi-parameter-based EMPC for determining the optimal dosage of magnesium hydroxide (Q Mg ) in a struvite crystallization unit. Simulations were carried out on the standard Benchmark Simulation Model 2d. (BSM2d.) to study the effect of the novel EMPC control strategy on operational cost and recovery of phosphorus of a struvite production unit installed at a full-scale WRRF.

BSM2d simulator
The Benchmark Simulation Model No. 2 (BSM2d) is a comprehensive plant-wide model describing several processes in a typical WRRF. BSM2d presents a realistic simulation environment for various operations in a full-scale WRRF. This stimulation standard is commonly used for evaluating operational sequences and control strategies in a WRRF. The implementation of BSM2d is available in most of the popularly used simulator platforms such as BioWin, GPS-X, Matlab/Simulink, Simba, STOAT, WEST, etc. (Gernaey et al., 2014). The plant layout, process model, influent data, test procedures, and evaluation criteria for simulator-based testing of control strategies are mentioned in the benchmarking standard manual (Nopens et al., 2010).

Struvite crystallizer in BSM2d
The original BSM2d model has been upgraded to include phosphorus transformation kinetics and unit operation for phosphorus recovery (Solon et al., 2017), (Kazadi-Mbamba et al., 2016). The process flow diagram and the location of the struvite crystallization unit in the updated BSM2d is presented in Fig. 1. In this new configuration, the supernatant obtained from dewatering the anaerobic digested sludge is supplied to a crystallization unit and the crystallizer overflow is recycled back to the anaerobic chambers of the biological process. The underlying reaction involved in struvite precipitation is presented in Equation (1).
Experiments have reported that the optimal pH range for struvite crystallization process is between 8.0 and 10.0 (Daneshgar et al., 2018). Therefore, the crystallizer is dosed with sodium hydroxide (Na(OH)) to maintain the pH values above 8.1. It should be noted that the process presented in Fig. 1 is one of several configurations capable of recovering phosphorus from wastewater sludge. The mathematical models describing the sludge treatment processes in the BSM2d might not be ideal, resulting in an overestimation of the phosphate concentrations in the digestates entering the struvite recovery unit. The current work focuses specifically on the struvite recovery process alone. Therefore, improvements in the unit processes preceding the struvite recovery unit, as well as modeling inadequacies associated with the integration of the struvite recovery unit in the updated BSM2d, are not addressed. Within the aforementioned limits, the steady-state operating conditions of the struvite unit and the mass balances for N and P occurring in the struvite recovery unit are presented as Fig. 4S in the supplementary material.

Economic Model Predictive Control
The Model Predictive Control (MPC) is a commonly used strategy for optimal control. In a conventional MPC, the outputs y k+i are predicted for the finite prediction-interval N P using a mathematical model and the control moves (u k , u k+1 , …u k+Nc− 1 ) are calculated for a control horizon N C to minimize an objective function J. The objective function for an error minimization control is calculated by a weighted square average of the control error and change in manipulated variable. In Equation (2) the term k and i are the time indices along the prediction horizon, r k+i is the reference value (set-point), w SP and w Δu are weights for the control error and change in manipulated variables respectively. Penalizing the control error by increasing the values of w SP in the objective function keeps the output variable close to the reference value. Increasing the value of w Δu suppresses rapid changes in the manipulated variables and makes the controller more sluggish.
The Economic model predictive control (EMPC) is a variant of the MPC where the cost function includes process performances, energy savings or overall economic profit rather than the quadratic error between the reference and measured variable. The successful implementation of EMPC for optimal operation of various wastewater treatment processes can be found in literature (Zeng and Liu, 2015;Zhang and Liu, 2019). The cost function used for EMPC control in the struvite recovery process is defined in Equation (3) Where PO 4,EFF is the phosphate concentration in the overflow from the crystallizer, PO 4,IN is the phosphate concentration in the influent, Q IN is the flowrate of supernatant to the crystallizer and P Recovery is the realtime estimate for the mass of phosphate recovered as struvite. Q Mg is the volumetric flowrate of magnesium hydroxide, MW Mg is the molecular weight of Mg and ρ Mg is the molar density of the magnesium hydroxide solution.φ STR is the market price of recovered phosphorus in the form of struvite and φ Mg the market price of magnesium hydroxide used in struvite production.

Prediction model
The control action of the MPC is taken based on the prediction made by the model. Therefore, adequate model describing the relation between input and output variables are imperative (Revollar et al., 2018). Struvite production in the crystallizer depends on several variables which introduce non-linear interdependencies in process chemistry. Physio-chemical models such as PHREEQC can adequately explain the process of struvite precipitation (Daneshgar et al., 2019). However, the lack of suitable sensors for measuring the concentration of every ionic species during the precipitation process poses a significant challenge in the use of mechanistic models for the purpose of automation and control. With the advent of data-driven models, various system identification techniques exist that can be used to establish a statistically significant correlation between the input and output variables of the process.

State-space representation
The discrete form of the linear state-space model is presented in Equations (6) and (7).
In Equations (6) and (7), x k is the state variable at time instance k, u k is the manipulated variable, d k is the measured disturbance and y k is the measured output. The list of input (u k = Q Mg and and output variables y k = PO 4,EFF are provided in Table 1.

Generation of subspace identification model
Several algorithms are mentioned in literature for the purpose of identifying a linear, time-invariant, state-space model from input-output data (Verhaegen, 1994;Van Overschee and De Moor, 2012). In our work, the canonical variate analysis (CVA) approach for system identification algorithm, mentioned in (Larimore, 1990;Ljung, 1999) was used for estimating the state-space matrices of the multiple input-single output (MISO) model. The data for system identification is obtained from simulations performed in the BSM2d simulator using the dynamic influent. The n4sid function, provided as a part of the System Identification Toolbox in MATLAB is used to train the model and obtain the A, B, C and D matrices.

Control strategy
Four different scenarios, each with a different Magnesium dosing control strategy were evaluated. The base control strategy (C0) uses constant dosage of Magnesium (Q Mg = 0.13 m 3 /day) and Sodium (Q Na = 0.10 m 3 /day). C0 is the default dosing strategy in the BSM2d simulator, which provides a reference point against which the basic and advanced control strategies are assessed. The second control strategy (C1) is a feed-forward controller, where the Q Mg is proportional to the flow of supernatant entering the crystallizer Q IN (Fig. 2a). A feedforward proportionality constant K P = 0.000675 was provided for the controller C1. Feed-forward flow-proportional control (C1) is a commonly used basic dosing control strategy adopted in most WRRFs (Ratnaweera H and Fettig, 2015). Assessing the performance of a basic control strategy helps provide an intermediate reference point to highlight the benefits of implementing an advanced optimal control. The controllers C2 and C3 are EMPC controllers described in Equations (3)-(5). In C2 the market price of struvite (φ STR ) and magnesium hydroxide (φ Mg ) are held constant during the evaluation period. In C3 the cost (φ STR ) and (φ Mg ) are chosen as time varying inputs to the EMPC controller. The control strategies for C1, C2 and C3 are shown in Fig. 2 and details regarding the measured and manipulated variables are presented as Table 1S in the supplementary material. In addition to control strategies for magnesium dose prediction, an on-off pH controller was provided in C1, C2 and C3 to ensure the pH values stay above 8.1. To study the difference between C2 and C3 control strategy (Fig. 2b.), a hypothetical scenario is considered where the cost of struvite (φ STR ) is changed once every 30 days. The new values of φ STR are randomly selected between the range φ STR,MIN = 7.5 and φ STR,MAX = 9.5 in the first day of every month. The Φ Mg values were held constant during the evaluation of control strategies to better understand the effects of monthly price variations on the control action of EMPC. It should be noted that the prices φ STR and φ Mg are mere representative values taken from literature (Solon et al., 2017) and do not reflect the exact price of struvite or magnesium hydroxide in the market.

Performance evaluation
Several standardized criteria for evaluating the performance of control strategies are reported in literature (Vanrolleghem et al., 1996). However, in this work we limit our evaluation to the parameters that are directly influenced by the struvite crystallization unit. The performance criteria used in our evaluations are explained in Equations (8)-(11).
M Str (kg /day) is the average per day value for the mass of struvite produced by the crystallizer unit during the evaluation period. Equation (8) describes the calculation of M STR .
M Mg is the average (kg/day) mass for Magnesium hydroxide consumed by the struvite crystallization unit.
Operational Cost Index (OCI) is a standard economic measure used to calculate the total cost (material and energy) incurred during the daily operation of a WRRF. Since this work focuses on optimizing struvite crystallization, a simpler version of the operational parameters Profit CRYST was also calculated using factors that has a direct influence on the cost of the struvite crystallization unit.
For the scenario with time-varying market price of struvite (φ STR (t)), the method for calculating Profit CRYST is presented in Equation (11).
EV TP is the number of times effluent total phosphorus limits (TP = 2 gP m − 3 ) are violated during the evaluation period.

Implementation in simulink platform
The simulink implementation files for the modified BSM2d simulator is available in literature (Solon et al., 2017). The basic (C1) as well as advanced (C2, C3) control strategies are constructed in the base open-loop BSM2d SIMULINK file (C0) with constant dose of magnesium. The nonlinear MPC block provided in simulink was used to construct the EMPC and configure the controller parameters for C2 and C3. Four separate files were created each with a different control strategy described as C0, C1, C2, and C3. The SIMULINK files for implementation of the EMPC control strategy in a BSM2d can be provided on request. The standard procedure executing the simulator are provided in the BSM2d simulator manual (Jeppsson et al., 2007). The following steps were applied for simulation and subsequent evaluation of the control strategies.
1. Initialize the BSM2d with the default values provided in the simulator. 2. Simulate from t = 0 to t = 300 days using the constant influent data. 3. Initialize the simulink model with the final values of the steady state simulation (using constant influent data). This allows the next simulation (with dynamic influent data) to begin at the exact position as where the steady state simulation had ended. 4. Simulate the model with the dynamic influent file for a period of 609 days (from t = 0 to t = 609 days). 5. The data from the simulations with dynamic influent file is stored in the MATLAB workspace. 6. Utilize the data recorded from t = 245 to t = 609 days to assess the performance parameters mentioned in Equations 8-11.

System identification results
The data received from the dynamic simulation of C1 is used to develop the state-space model. The timeseries data used in the subspace identification method contains Q Mg , Q IN , and PO 4,IN as model inputs and PO 4,EFF as model output, with a time span of 609 days (t = 0 to t = 609 days) sampled at an interval of 15 min. The data is split to training data 75% (t = 0 to t = 455), which provides adequate number of data points necessary to construct a reliable model. The remaining 25% of the data is chosen as validation data (t = 455 to t = 609) which was used to confirm the ability of model to predict PO 4,EFF in the overhead stream of struvite recovery unit. The training data is used to calibrate the model and obtain the A, B, C, and D matrices (Eqs (6) and (7)) defined in the state-space matrix. The model is then used to predict the output for the validation data. The validation plot, showing a comparison between the validation dataset and the value predicted by the state-space model is presented in Fig. 3.
A close match is observed between the measured data and the data predicted by the model. The adequate match between the model predicted and plant data as well as the R 2 value of 0.91 demonstrate the ability of the state-space model to predict the effluent PO 4 -P concentrations in the overhead streams from the crystallizer.

MPC settings
The choice of prediction and control horizon, limits of manipulated variables determine the performance of an MPC. Systematic procedures are explained to determine the optimal values of these tuning parameters (Lee and Yu 1994). In our case, the parameters were tuned based on experience gained from running the BSM2d simulator with steady-state as well as the dynamic weather data. The MPC parameters used in the simulations are presented in Table 1.
The time-step in the influent disturbance file, the sampling interval of data used in prediction model generation, and the logging rate of simulations results in the MATLAB workspace were set at 15 min. In order to maintain uniformity, the same value was also used as EMPC's time-step. The Q Mg values in dynamic simulation using C1 strategy varies within the range 0.10 m 3 /day to 0.20 m 3 /day. In the EMPC strategy, range of variations in Q Mg was increased by reducing the lower limit value of 0.1 m 3 /day by 50% (− 0.05 m 3 /day) and increasing upper limit value of 0.20 m 3 /day by 50% (+0.10 m 3 /day). Therefore, the parameters u min and u max (Table 1) were set as 0.05 m 3 /day and 0.3 m 3 / day respectively. The steady state-simulations were used to determine the appropriate values of N P and N C .
Step changes were provided in the values of Φ STR and the performance of the EMPC was assessed by recording its ability to reach a new minimum point. At each simulation run the values of prediction and control horizon were gradually increased from a lower value N P = 2 and N C = 1, and the controller performance was evaluated at every step. The controller performance showed no significant change when the values were increased above, N P = 3 and N P = 2. Therefore, the values of prediction and control horizon were maintained at 3 and 2 respectively.

Controller performance evaluation
A quantitative assessment is necessary to compare the of the performance of various control strategies presented in Fig. 2. The wastequality parameters of the influent and effluent streams of the struvite crystallization unit for all four strategies are utilized to generate comparison plots and calculate the performance evaluation indicators mentioned in Eqs. 8-11. The dynamic values of the soluble PO 4 -P in the overhead flow of the struvite crystallization unit for the constant dosing scenario (C0) is presented in Fig. 4.
The benefits of introducing a struvite recovery unit in terms of improving effluent water-quality and reducing maintenance costs in WRRF is already reported in literature (Mbamba et al., 2016). However, Fig. 4 indicates that maintaining a constant magnesium hydroxide flowrate (C0) in the crystallizer would not be the most optimal dosing strategy. Several underdosing points (marked in circles) are indicated where increasing the dosage could have resulted in higher phosphorus recovery in the form of struvite.   tallizer. The improvement in dosing strategy can be reaffirmed in Table 2., which indicates an increase (8.01%) in average phosphorus recovery while reducing (8.62%) the Magnesium consumption. However, the dose prediction is entirely based on the flowrate, and flowproportional control strategy (C1) does not consider the fluctuations in the PO 4 -P concentrations in the influent. Therefore, suboptimal dosing is observed in situations with higher influent PO 4 -P fluctuations. The EMPC predicts the dose based on the optimal value of cost function (Equation (3)). Since the dose prediction considers both flowrate as well as influent PO 4 -P concentrations, a better control over the recovery of phosphorus can be expected. Fig. 5 shows a more stable effluent PO 4 -P concentration in C2 compared to C0 and C1. The dynamic plots for Magnesium consumed, recovery percentage and mass of struvite production are presented in supplementary material. A comparison between the performance indices presented in Table 2 indicates a 12.5% increase in the average daily struvite production and an 8.5% drop in total magnesium consumption compared to the base dosing control strategy (C0).
A comparison between the control strategies in terms of overall profits for struvite production indicate a 33.01% increase for C1 and 44.03% increase for C2 when compared to the base control strategy C0. The increase in overall profits are primarily due to improvements in phosphorus recovery and a reduction in magnesium used in the crystallizer by avoiding overdosing. Apart from the increase in profits for struvite production, implementing EMPC also demonstrates fewer effluent violations in the treated effluent from WRRF (Table 2). Fewer effluent violation would imply fewer effluent penalties, which could further add to the savings.

Economic assessment of time varying cost function
The influence of introducing a time varying cost function on the effluent PO 4 -P concentration and the magnesium dose prediction is presented in Fig. 6.
It is observed that in situations with lower φ STR , the C3 strategy predicts a lower magnesium dosage and a subsequent reduction in phosphorus recovery. The reduction in struvite prices (reflected in φ STR ) moved the optimal dosing point to a lower value (resulting in lower phosphorus recovery percentage) in order to generate savings on the magnesium cost. When the costs are increased, higher recovery of phosphorus (in the form of struvite) is restored. A comparison between the average (kg/day) values of struvite produced, magnesium consumed and the profits for all three control strategies in the presence of a monthly variations in struvite costs are presented in Table 3. It is observed that using costs as inputs to the EMPC controller as opposed to a constant cost function results in an additional 8.1% increment in the overall profits incurred in struvite production.
A complete list of evaluation parameters for all the three control strategies are presented in the supplementary material. The EMPC strategy also provides a convenient alternative to adjust the controller parameters based on the price factors without the need for retuning the controller parameters. Although the work primarily defines the optimization problem based on two parameters a.) Cost of Struvite b.) Cost of Magnesium, EMPC also provides the possibly of including more optimization parameters (Shaddel et al., 2019a) as long as a reliable correlation between the manipulated variables and the measured parameters exist.

Conclusion
The work demonstrates the advantages of a multi-parameter-based control strategy for optimal dosing of magnesium in a struvite crystallization unit. A systematic procedure for developing a data-driven model for establishing a correlation between the input and output parameters of the struvite crystallization process has also been presented. Performance evaluation of the EMPC indicates a significant improvement in the overall profits when compared to both constant as well as flow proportional dosing strategy. The operational flexibility of the EMPC controller was demonstrated by its ability to conveniently switch between multiple operating conditions by using the market price of struvite and magnesium as their input variables. Although the work demonstrated the optimization strategy based merely on two economic parameters (magnesium dose flowrate and phosphorus recovery); the flexible nature of the EMPC allows the possibility of introducing   multiple evaluation criteria in the objective function. The multiparameter based optimal control approach using a data-driven model presents an opportunity to further improve the process operation and achieve better product quality by deploying the optimal dose strategy based on criteria such as struvite crystal dimensions, settling properties etc.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.