Real-Time Optimization of Wastewater Treatment Plants via Constraint Adaptation

Abstract

An important requirement of wastewater treatment plants (WWTPs) is compliance with the local regulations on effluent discharge, which are going to become more stringent in the future. The operation of WWTPs exhibits a trade-off between operational cost and effluent quality, which provides a scope for optimization. Process optimization is usually done by optimizing a model of the process. However, due to inevitable plant–model mismatch, the computed optimal solution is usually not optimal for the plant. This study represents the first attempt to handle plant–model mismatch via constraint adaptation (CA) for the real-time optimization of WWTPs. In this simulation study, the “plant” is a model adopted from the BSM1 benchmark, while a reduced-order “model” is used for making predictions and computing the optimal inputs. A first implementation uses steady-state measurements of the plant constraints to adjust the model in the optimization framework. A fast CA technique is also proposed, which adjusts the model using transient measurements. It is observed that, even in the presence of significant plant–model mismatch, the two proposed techniques are able to meet the active plant constraints. These techniques are found to reduce the pumping and aeration energy by 20%, as compared to that adopted in BSM1.

1. Introduction

Wastewater treatment plants (WWTPs) are large-scale facilities for restoring pol-luted wastewater to a desirable quality. Wastewater discharge standards are enforced to ensure a reduced environmental impact, the compliance of which are getting stricter with time. These stricter standards impose an extra burden on WWTPs, thus leading to a significant increase in their operational costs [1,2,3]. In the operation of a WWTP, most of the energy is consumed during the secondary treatment, where the majority of nutrients and organic matter is removed with the help of an activated sludge process. This energy, which consists of pumping energy for the recycle flows and aeration energy to provide oxygen to the aerobic processes, consumes 50–90% of the total electricity required for the plant [4,5,6]. The trade-off between cost (energy consumption) and treatment performance (effluent quality) can be addressed via optimization. In other words, optimal operation of WWTPs implies minimizing the operational cost, while complying with the wastewater discharge standards.

A major issue associated with the implementation of optimization techniques is casting a real-world problem into a mathematical model. Typically, a mathematical model cannot describe a real process (labeled here plant) accurately due to assumptions, approximations and simplifications associated with the modeling process. It follows that optimal solutions obtained by solving a model-based optimization problem are optimal for the model, but usually not for the plant due to inevitable plant–model mismatch. Furthermore, if these model-based solutions are directly implemented on the plant, they may lead to violation of some key constraints. Unfortunately, this fact is too often overlooked by academic researchers. For example, there are several studies that propose to optimize the performance of WWTPs by solving a model-based optimization problem [6,7,8,9,10,11]. Additionally, some studies use model predictive controllers that rely on solving a model-based optimization problem [12,13,14,15,16]. Note, however, that the limitations of model-based approaches can be overcome by including plant measurements in the optimization framework to handle plant–model mismatch. These measurement-based optimization techniques are referred to as real-time optimization (RTO) in the literature [17,18].

Ideally, RTO techniques use plant measurements to drive the plant to optimality, while ensuring constraint satisfaction. Measurements can be used at three different levels, namely, the process model (to update the model parameters), the cost and constraints (to modify the cost and constraints of the optimization problem), and the constraints (to modify only the bias in the constraints of the optimization problem) as discussed next:

• Parametric uncertainties can be handled by updating the uncertain model parameters using a so-called “two-step” approach. This works well in the case of parametric plant-model mismatch, that is, when the plant-model mismatch can be corrected by parameter adaptation [19]. In the case of structural mismatch, parameter estimation will not be able to adapt the model correctly, and the two-step scheme may not work well [20].
• Structural uncertainty can be handled by modifying the cost and constraints of the optimization problem. The first investigation in this direction is the method labeled integrated system optimization and parameter estimation (ISOPE) [21]. However, the parameter estimation step makes its implementation difficult in the presence of measurement noise and insufficient excitation in the data. Hence, it may be beneficial to skip the parameter estimation step and modify only the cost and constraints in the optimization framework. The corrections added to the cost and constraints are referred to as modifier terms, and the resulting RTO scheme is named modifier adaptation (MA) [22]. The correction of bias and gradients in the predicted cost and constraints are performed via so-called zeroth- and first-order modifiers, respectively.
• When only zeroth-order modifiers are used to compensate for plant-model mismatch, the RTO method is known as constraint adaptation (CA) [23]. Note that gradient estimation is difficult in practice, especially with a large number of inputs. CA does not involve gradient estimation and is often preferred over MA. It ensures feasible operation of the plant by driving it iteratively to a point where all constraints are met. However, this way, optimal plant operation can only be ensured when the number of inputs (or decision variables) is equal to the number of active plant constraints [18].

This paper describes the application of CA towards the optimization of WWTPs. To the best of our knowledge, there is no work in the literature using CA for the real-time optimization of WWTPs. CA is an iterative scheme that ensures feasibility, but not necessarily optimality, upon convergence. However, in this study, since the optimal operation is characterized by two quality requirements (or constraints) being met, CA can lead to plant optimality by working with only two inputs.

Hence, CA is capable of optimizing the BSM1 plant operation at steady state using an adequate (but not necessarily accurate) model and appropriate either steady-state or transient measurements (There is a noticeable difference between model accuracy and model adequacy. In the context of MA and CA, good performance requires model adequacy, but not necessarily model accuracy. Model adequacy for CA simply says that the model is capable of predicting the correct set of active constraints at the plant optimum. See A.G. Marchetti, T. de Avila Ferreira, S. Costello and D. Bonvin, Modifier Adaptation as a Feedback Control Scheme, Ind. Eng. Chem. Res., 59(6), 2261–2274 (2020) for more details.). The reason for which plant optimality can be reached even with an inaccurate model is the fact that optimization is set up in a way that (i) the model is corrected using bias terms on the constraints so as to be able to predict the active plant constraints correctly, and (ii) all degrees of freedom in the optimization problem are used to meet the active constraints [23]. Note that, this way, a higher priority is given to feasibility (meeting the active plant constraints) than to optimality, which corresponds to industrial practice.

The paper is organized as follows. Section 2 on Materials and Methods describes the various WWTP models, namely, the 147th-order BSM1 model, the 47th-order model that is used as “plant”, and the 20th-order model that is used as “model”. Then, the RTO methodologies CA and fast CA are briefly presented, first in general terms and then with respect to their application to WWTP. Section 3 on Results and Discussion presents the main results obtained with both CA and fast CA for the case of constant influents. Conclusions are provided in Section 4. In addition, there are three appendices supporting the model simplifications introduced to reduce the 147th-order BSM1 model to the 20th-order “model”.

2. Materials and Methods

In silico approaches are often used to evaluate the performance of control and optimization schemes without expensive experimentation. For this, a simulated reality is typically selected as the “plant”, as was done in [15,16]. These studies have considered the Benchmark Simulation Model No. 1 (BSM1) [24] as “plant” and a reduced-order model as “model”. In our work, for simplicity, a 47th-order model adopted from the Benchmark Simulation Model No. 1 (BSM1) [24] is selected as “plant” and a reduced 20th-order model is used as “model” in CA schemes. The reason for this choice is twofold: (i) numerical optimization of the BSM1 model that includes the settler equations is a rather difficult task, and (ii) the actual performance of CA depends on the ability of a modified model to appropriately predicting the measured plant constraints; this is easily done by the use of bias terms that express the differences between the plant measurements and the values predicted by the model.

2.1. Wastewater Treatment Plant: Models and Controllers

Wastewater treatment is carried out in three main stages, namely, primary, secondary and tertiary treatments [25]. This study focuses on the secondary treatment that involves biological processes and secondary settling. The secondary treatment is also referred to as the activated sludge process (ASP) in the literature [26,27,28,29,30]. The standard benchmark model named BSM1 is widely used among the research community for simulation-based comparative studies.

2.1.1. Benchmark Simulation Model No. 1

The BSM1 model represents the activated sludge process in WWTPs, where the biological process is modeled using the Activated Sludge Model No.1 (ASM1) [31] and the settling process is modeled using the double-exponential settling velocity function [32]. BSM1, like any conventional secondary treatment in a WWTP, has a series of biological reactors followed by a settling tank. BSM1 includes 5 reactors, of which the first two are under anoxic conditions, and the remaining three are under aerobic conditions. The anoxic section, the aerobic section and the settler have capacities of 2000, 4000 and 6000 ${\mathrm{m}}^3$, respectively. There are two recycle streams, namely, one from the last tank (labeled $Q_a$), and another one from the underflow of the settler (labeled $Q_r$). Furthermore, a WWTP typically involves two built-in control loops to regulate (i) the nitrate level in the second reactor,$S_{NOp,2}$, by manipulating $Q_{ap}$, and (ii) the oxygen level in the fifth reactor,$S_{Op,5}$, by manipulating the oxygen transfer coefficient ${\left(K_{L}a\right)}_p$ in the aerobic section. This can be done using, for example, two proportional-integral (PI) controllers. The schematic representation of the BSM1 plant is shown in Figure 1.

2.1.2. 47th-Order Plant Adopted from BSM1

In the BSM1 model, every biological reactor is described using 13 ordinary differential equations (ODEs). Hence, the series of five biological reactors contribute 65 ODEs to BSM1. The secondary settler is a ten-layer system, where each layer is represented by 8 ODEs. Hence, the secondary settler contributes 80 ODEs to BSM1. Furthermore, each controller adds one differential equation to the model if it possesses an integral term. Consequently, BSM1 has a total of 147 ODEs.

The settling process in BSM1 is assumed to be biologically inactive, which means there is no biological phenomena taking place in the settling tank, and the only process prevalent in the settler is the settling of solids. Consequently, the soluble matter concentrations are not affected by the settling process and are the same at both the inlet and the outlet of the settler. Moreover, the settler overflow has a negligible amount of particulate matter in the treated effluent, as can be seen in the simulated results presented in [24]. It follows that one can assume that the settler generates effluent with no particulate matter, which allows eliminating 80 ODEs from the BSM1 model. Hence, the assumption of a perfect settler in BSM1 reduces the number of ODEs from 147 to 67. The perfect settler equations are shown in Appendix A.

BSM1 can further be simplified by removing 20 additional ODEs from the model as each reactor contains two inert states and two so-called reaction invariants. Reaction invariants are state variables that are not directly affected by the reaction [33]. $X_p$ and $S_{Alk}$ are reaction invariants in ASM1 as they do not appear in the process rates shown in the ASM1 Peterson matrix presented in [31]. ASM1 also includes two inert states, $S_I$ and $X_I$, which do not participate in the reactions and leave the system at their inlet concentrations. Removing the mass balance equations for these 4 states in each of the 5 biological reactor results in a 47th-order model. The representation of the reduced ASM1 model with 9 state variables in the form of Peterson matrix is shown in Appendix B. The 47th-order model presents two main advantages. Firstly, it is a good representation of BSM1 as its predictions are fairly close to those of the full model as shown in

Table A1. Steady-state values of variables of interest of the 147th-order and 47th-order models at three different set of inputs.

Variables	Units	147th-Order Model	47th-Order Model	147th-Order Model	47th-Order Model	147th-Order Model	47th-Order Model
$Q_{ap}$	${\mathrm{m}}^3/\mathrm{d}$	55,338	55,338	36,892	36,892	18,446	18,446
${\left(K_{L}a\right)}_p$	$1/\mathrm{d}$	240	240	200	200	160	160
$S_{NOp,2}$	$\mathrm{gN}/{\mathrm{m}}^3$	8.54	8.35	4.27	4.00	0.28	0.23
$S_{Op,2}$	$\mathrm{gCOD}/{\mathrm{m}}^3$	3.93	3.98	2.97	2.99	1.10	0.96
$S_{NH,2}$	$\mathrm{gN}/{\mathrm{m}}^3$	0.69	0.55	0.89	0.72	3.3	3.00

. Secondly, it can be easily optimized owing to the assumption of a perfect settler.

2.1.3. 20th-Order Model

The 47th-order “plant” model is still too large to be used for real-time computations. The model can be further simplified by combining the two anoxic reactors into a single reactor of volume 2000 ${\mathrm{m}}^3$ and the three aerobic reactors into a single reactor of volume 4000 ${\mathrm{m}}^3$. This way, a reduced model having only two reactors is obtained, with each reactor being represented by 9 ODEs. The reduced model includes 20 ODEs, 18 ODEs from the reactors and 2 ODEs from the PI controllers. This 20th-order model will be used as the “model” in the RTO scheme. The schematic representation of this controlled model is shown in Figure 2.

2.1.4. Design of Controllers

The general equation for a PI controller is given as [34]:

$$u\left(t\right)=K_c\left(e\left(t\right)+\frac{1}{{\tau }_I}{\int }_0^{t}e\left(t\right)dt\right)$$

(1)

where u is the manipulated variable, e is the difference between the setpoint (SP) and the process variable (PV), $K_c$ is the controller gain and ${\tau }_I$ is the integral time constant.

Differentiating Equation (1) with respect to t gives the following ODE:

$$\dot{u}\left(t\right)=K_c\left(\dot{e}\left(t\right)+\frac{1}{{\tau }_I}e\left(t\right)\right)$$

(2)

Based on this equation, the controller equations for the plant and the model are provided in Appendix C. Controller tuning can be done using, for example, the Ziegler–Nichols step response method [34].

2.2. RTO via Constraint Adaptation

2.2.1. Formulation of the Optimization Problem

Let the steady-state plant optimization problem be expressed mathematically as [35]:

$$\begin{array}{c}{\mathit{u}}_p^{*}≔arg \underset{u}{\min }{\mathsf{\Phi }}_p\left(\mathit{u}\right)≔{\varphi }_p\left(\mathit{u},{\mathit{y}}_p\left(\mathit{u}\right)\right)\\ s.t. {\mathit{G}}_p\left(\mathit{u}\right)≔{\mathit{g}}_p\left(\mathit{u},{\mathit{y}}_p\left(\mathit{u}\right)\right)\le \mathit{0}\\ {\mathit{u}}^L\le \mathit{u}\le {\mathit{u}}^U\end{array}$$

(3)

where $\mathit{u}$ is the input vector of dimension $n_u$, ${\mathit{y}}_p$ is the measured plant output vector of dimension $n_y$, ${\varphi }_p$ is the plant cost function, ${\mathit{g}}_p$ is the vector of plant constraints of dimension n_G, and subscript ${\left(·\right)}_p$ indicates a quantity related to the plant. Furthermore, ${\mathsf{\Phi }}_p$ and ${\mathit{G}}_p$ are the plant cost function and vector of plant constraints expressed in terms of $\mathit{u}$, respectively.

It is clear in the above definitions that the plant cost and constraints are known function of $\mathit{u}$ and ${\mathit{y}}_p$. However, in practice, the plant steady-state input–output mapping ${\mathit{y}}_p\left(\mathit{u}\right)$ is typically unknown, and it is estimated using an approximate steady-state model:

$$\mathit{f}\left(\mathit{x},\mathit{u},\mathit{\theta }\right)=\mathit{0}$$

(4)

$$\mathit{y}=\mathit{h}\left(\mathit{x},\mathit{u},\mathit{\theta }\right)$$

(5)

where $\mathit{x}$ is the state vector of dimension $n_x$, $\mathit{y}$ is the model predicted outputs vector of dimension $n_y$, and $\mathit{\theta }$ is the model parameters vector of dimension $n_{\theta }$.

Solution to $\mathit{f}\left(\mathit{x},\mathit{u}, \mathit{\theta }\right)=\mathit{0}$ can be expressed as $\mathit{x}=\mathit{\xi }\left(\mathit{u}, \mathit{\theta }\right)$ with $\mathit{\xi }$ the steady-state map between $\mathit{u}$ and $\mathit{x}$. Hence, the input–output map predicted by the model can be represented as:

$$\mathit{y}\left(\mathit{u},\mathit{\theta }\right)=\mathit{h}\left(\mathit{\xi }\left(\mathit{u},\mathit{\theta }\right),\mathit{u},\mathit{\theta }\right)=\mathit{H}\left(\mathit{u},\mathit{\theta }\right)$$

(6)

This way, a model-based optimization problem can be formulated as follows:

$$\begin{array}{c}{\mathit{u}}^{*}≔\mathit{arg}\underset{u}{\min }\varphi \left(\mathit{u},\mathit{\theta }\right)≔\varphi \left(\mathit{u},\mathit{y}\left(\mathit{u},\mathit{\theta }\right)\right)\\ s.t.\mathit{G}\left(\mathit{u},\mathit{\theta }\right)≔\mathit{g}\left(\mathit{u},\mathit{y}\left(\mathit{u},\mathit{\theta }\right)\right)\le \mathit{0}\\ {\mathit{u}}^L\le \mathit{u}\le {\mathit{u}}^U\end{array}$$

(7)

where $\varphi$ is the model cost function, $\mathit{g}$ is the vector of model constraints of dimension $n_G$, $\mathsf{\Phi }$ is the plant cost function expressed in terms of $\mathit{u}$, and $\mathit{G}$ is the vector of plant constraints expressed in terms of $\mathit{u}$.

The model-based optimization problem Equation (7) can be solved numerically. However, the model optimum ${\mathit{u}}^{*}$ does not coincide with plant optimum ${\mathit{u}}_p^{*}$ due to plant–model mismatch. Hence, plant measurements are incorporated into the optimization framework to tackle this plant–model mismatch. In CA, these measurements are used to construct zeroth-order modifiers. Zeroth-order modifiers are used for bias correction in the predicted constraints. They represent the difference between the measured plant constraints and the constraints predicted by the model for given inputs. There is no need to add zeroth-order modifier to the cost function since adding a constant value to the cost function does not change the location of the optimal solution. Hence, it is used only with the constraints and computed as follows:

$${\mathit{ϵ}}_k^G={\mathit{G}}_p\left({\mathit{u}}_k\right)-\mathit{G}\left({\mathit{u}}_k,\mathit{\theta }\right)$$

(8)

where ${\mathit{ϵ}}_k^G$ is a vector of dimension $n_G$, and $k$ is the iteration number.

2.2.2. CA Using Steady-State Measurements

Using the aforementioned zeroth-order modifiers, a modified optimization problem is solved iteratively as per Figure 3. At $k$th RTO iteration, ${\mathit{u}}_k$ is obtained from the numerical solution to the following modified optimization problem:

$$\begin{array}{c}{\mathit{u}}_{k+1}≔\arg \underset{u}{\min }\mathsf{\Phi }\left(\mathit{u}\right)\\ s.t.{\mathit{G}}_m\left(\mathit{u}\right)≔\mathit{G}\left(\mathit{u}\right)+{\mathit{ϵ}}_{f,k}^G\le \mathit{0}\end{array}$$

(9)

where ${\mathit{ϵ}}_{f,k}^G$ are filtered modifiers as discussed next, and subscript ${\left(·\right)}_m$ indicates a quantity that has been modified. Indeed, using modifiers Equation (8) often results in overcorrection, thereby causing the CA scheme to oscillate. This can be avoided by applying a first-order filter on the modifiers as follows:

$${\mathit{ϵ}}_{f,k}^G=\left(\mathit{I}-\mathit{K}\right){\mathit{ϵ}}_{f,k-1}^G+\mathit{K}{\mathit{ϵ}}_k^G$$

(10)

$\mathit{I}$ is the identity matrix, while $\mathit{K}$ is a filter matrix whose eigenvalues are between 0 (full filtering) and 1 (no filtering). The matrix $\mathit{K}$ is often chosen as a diagonal matrix, with $n_G$ identical eigenvalues, that is, $\mathit{K}=K\mathit{I}$. The more noise, the lower the K value. Values between 0.5 and 0.8 are often used. The effect of more filtering is to reduce the convergence speed, but it does not affect the converged values. The CA algorithm is given next:

• Initialize the filtered modifiers ${\mathit{ϵ}}_{f,0}^G=\mathit{0}$. Choose the filter matrix $\mathit{K}$ and the convergence threshold δ. Set $k≔0$.
• Solve the modified optimization problem Equation (9) to get the inputs ${\mathit{u}}_{k+1}$.
• If $\Vert {\mathit{u}}_{k+1}-{\mathit{u}}_k\Vert \le \delta$, set $k≔k+1$; otherwise stop.
• Apply the inputs to the plant, wait for steady state and take the measurements ${\mathit{G}}_p\left({\mathit{u}}_k\right)$.
• Evaluate the modifiers Equation (8), apply the filter Equation (10) and return to Step 2.

If the above RTO scheme converges, it will satisfy the plant constraints. For instance, if the RTO scheme converges for $k\to \infty$, then the modified quantities for $k\to \infty$ correspond to the plant quantities, and ${\mathit{G}}_m$ approaches ${\mathit{G}}_p$ as shown below:

$$\begin{array}{c}{\mathit{G}}_m\left(\mathit{u},\mathit{y}\left({\mathit{u}}_{\infty },\mathit{\theta }\right)\right)≔ \mathit{G}\left({\mathit{u}}_{\infty },\mathit{y}\left({\mathit{u}}_{\infty },\mathit{\theta }\right)\right)+{\mathit{ϵ}}_{\infty }^G\hfill \\ {\mathit{G}}_m\left(\mathit{u},\mathit{y}\left({\mathit{u}}_{\infty },\mathit{\theta }\right)\right)≔ \mathit{G}\left({\mathit{u}}_{\infty },\mathit{y}\left({\mathit{u}}_{\infty },\mathit{\theta }\right)\right)+{\mathit{G}}_p\left({\mathit{u}}_{\infty },{\mathit{y}}_p\left({\mathit{u}}_{\infty }\right)\right)-\mathit{G}\left({\mathit{u}}_{\infty },\mathit{y}\left({\mathit{u}}_{\infty },\mathit{\theta }\right)\right)\\ {\mathit{G}}_m\left(\mathit{u},\mathit{y}\left({\mathit{u}}_{\infty },\mathit{\theta }\right)\right)≔ {\mathit{G}}_p\left({\mathit{u}}_{\infty },{\mathit{y}}_p\left({\mathit{u}}_{\infty }\right)\right)\end{array}$$

2.2.3. Fast CA Using Transient Measurements

CA solves a modified optimization problem that uses steady-state plant measurements. When the plant requires much time to reach steady-state, this results in slow convergence to plant optimality. This issue can be mitigated by estimating the plant steady state using transient measurements and a dynamic model as shown in Figure 4 [36,37]. For this, at time $t_k$ corresponding to the $k$th RTO iteration during transient, the plant measurements $y_p^{dyn}$ are taken and complemented by the model prediction of the remaining part of the transient to steady state. This prediction $\mathsf{\Delta }$ is obtained via a dynamic model of the plant. Hence, an estimate of the plant steady-state $y_p^{ss}$ is obtained as:

$${\widehat{y}}_p^{SS}=y_p^{dyn}+\mathsf{\Delta }=y_p^{dyn}+\left(y^{SS}-y^{dyn}\right)$$

(11)

The estimated plant steady state is used in the CA scheme, which can then converge in a single settling time (i.e., the time it takes for the plant to reach steady state). CA using transient measurements is referred to as fast CA in the literature [38].

Both CA and fast CA ensure feasible operation of the plant by driving it iteratively to a point where all the constraints are met. In some cases, it also provides plant optimality, in particular when the number of active constraints is equal to the number of inputs [23]. Note that both CA and fast CA need to work iteratively because they compute an optimum locally (and not globally) at each iteration.

2.3. RTO Applied to WWTP

This section describes the application of RTO to a WWTP as per the scheme shown in Figure 3. The optimization problem is solved by incorporating the model described in Section 2.1.3. as equality constraints. This model is also used to provide the model predictions needed to compute the modifiers. The plant measurements required in the RTO scheme are generated by simulating the plant described in Section 2.1.2. In addition, a plant optimization problem is solved to determine the plant optimum, which is required to verify that the RTO scheme reaches plant optimality. As highlighted earlier, the mechanical pumping and the aeration process are the two most energy-intensive processes in WWTPs. Hence, the objective of optimization is to minimize the sum of pumping energy (PE) and aeration energy (AE) of the WWTP as presented below.

PE, which represents the energy used for pumping the internal recycle stream $Q_a$, is calculated as follows:

$$\mathrm{PE}=f_{PE}·Q_a\left[\mathrm{kWh}·{\mathrm{d}}^{-1}\right]$$

(12)

where $f_{PE}$ is the PE factor and has the value of 0.004 $\mathrm{kWh}·{\mathrm{m}}^{-3}$ as given in [24].

AE is expressed in terms of the oxygen transfer coefficient $K_{L}a$ and represents the energy needed to transfer oxygen into the reactor. $K_{L}a$ is directly related to energy-consuming factors such as airflow and mixing velocity [39]. It is assumed that 1800 ${\mathrm{gO}}_2$ can be transferred per $\mathrm{kWh}$ of energy at 15 °C and zero dissolved oxygen concentration [24]. This assumption leads to the following oxygen transfer rate (OTR) expression:

$$\mathrm{OTR}=V_{aerobic}·K_{L}a·S_{O,sat}\left[{\mathrm{gO}}_2·{\mathrm{d}}^{-1}\right]$$

(13)

where $V_{aerobic}$ is the volume of the reactor under aerobic conditions expressed in ${\mathrm{m}}^{-3}$, and $S_{O,sat}$ is the oxygen saturation concentration with the value of 8 ${\mathrm{gO}}_2·{\mathrm{m}}^{-3}$ at 15 °C. Hence, AE can be calculated as follows:

$$\mathrm{AE}=f_{AE}·\mathrm{OTR}\left[\mathrm{kWh}·{\mathrm{d}}^{-1}\right]$$

(14)

where $f_{AE}$ is the AE factor and has the value of $\frac{1}{1800}$ $\mathrm{kWh}·{\left({\mathrm{g}\mathrm{O}}_2\right)}^{-1}$ as provided in [24].

Finally, the overall cost index (OCI) is expressed as the sum of PE and AE:

$$\mathrm{OCI}=0.004Q_a+\frac{S_{O,sat}}{1800}·V_{aerobic}·K_{L}a\left[\mathrm{kWh}·{\mathrm{d}}^{-1}\right]$$

(15)

Substituting $S_{O,sat}=8$ and $V_{aerobic}=4000$ in the above expression gives:

$$\mathrm{OCI}=0.004Q_a+17.8K_{L}a\left[\mathrm{kWh}·{\mathrm{d}}^{-1}\right]$$

(16)

The constraints for the optimization problems include bounds on the manipulated variables $K_{L}a$ and $Q_a$ and on the setpoints of the controlled variables $sp{S}_{NO}$ and $sp{S}_O$. Additionally, a constraint is also placed on the ammonia concentration in the effluent. Constraint adaptation guarantees feasible operation upon convergence; therefore, the bounds are not going to be violated. Hence, the effluent quality (EQ) is automatically enforced by the upper bounds on the concentrations. The various optimization problems can be formulated as discussed next.

2.3.1. Plant Optimization Problem

The plant optimization problem is solved to have some prior knowledge of plant optimality in this simulation study. The plant optimum will be used to validate the optimum obtained using RTO schemes. The plant optimization problem reads:

$$\begin{array}{c}\underset{sp{S}_{NOp,2},sp{S}_{Op,5}}{\min }0.004Q_{ap}+17.8{(K_{L}a)}_p\\ 45SSplantequations\\ 2SScontrollerequations\\ 180\le Q_{ap}\le 92,000\\ 1.5\le {(K_{L}a)}_p\le 360\\ 0.5\le sp{S}_{NOp,2}\le 5\\ 0.5\le sp{S}_{Op,5}\le 5\\ S_{NHp,5}\le 4\end{array}$$

(17)

This optimization problem is formulated with the setpoints $sp{S}_{NOp,2}$ and $sp{S}_{Op,5}$ as decision variables and includes bounds on the inputs $Q_{ap}$ and ${\left(K_{L}a\right)}_p$. These bounds are adopted from [24]. In Equation (17), SS stands for steady state and the optimization problem includes 45 steady-state plant equations along with 2 controller equations as constraints.

2.3.2. Model Optimization Problem

The model optimization problem includes 18 steady-state model equations along with 2 steady-state controller equations as constraints as represented below:

$$\begin{array}{c}\underset{sp{S}_{NO,1},sp{S}_{O,2}}{\min }0.004Q_a+17.8K_{La}\\ 18SSmodelequations\\ 2SScontrollerequations\\ 180\le Q_a\le 92,000\\ 1.5\le K_{L}a\le 360\\ 0.5\le sp{S}_{NO,1}\le 5\\ 0.5\le sp{S}_{O,2}\le 5\\ S_{NH,2}\le 4\end{array}$$

(18)

The bounds on inputs and setpoints of the model optimization problem are taken similar to that of the plant as described in the Section 2.3.1. In the presence of plant-model mismatch, the optimal solution to the above model-based problem does not typically lead to plant optimality. Hence, the model-based problem will be modified and used in the CA and fast CA schemes to drive the plant towards optimality.

2.3.3. Modified Optimization Problem

The following modified optimization problem is solved at each RTO iteration $k$ to compute the optimal setpoints for the next iteration:

$$\begin{array}{c}\underset{sp{S}_{NO,1},sp{S}_{O,2}}{\min }0.004Q_a+17.8K_{L}a\\ 18SSmodelequations\\ 2SScontrollerequations\\ 180\le Q_a\le 92,000\\ 1.5\le K_{L}a\le 360\\ 0.5\le sp{S}_{NO,1}\le 5\\ 0.5\le sp{S}_{O,2}\le 5\\ S_{NH,2}+{ϵ}_k\le 4\end{array}$$

(19)

When some of the constrained quantities are either perfectly known or known functions of the inputs, they possess no uncertainty and thus they do not require modification [40]. The bounds on the inputs $Q_a$ and $K_{L}a$ represent saturations of the control actions. These saturation limits are perfectly known. In addition, the bounds on the setpoints are also perfectly known. Consequently, the bounds on inputs and setpoints need not be modified in the RTO scheme. Hence, $S_{NH,2}$ is the only quantity that needs to be modified via zeroth-order modifier as described next.

• Zeroth-order modifier for CA

The inequality constraint $G$ in a minimization problem formulation is typically provided by specifying an upper bound at zero, that is, $G\le 0$. Consequently, the upper bounds on the ammonia concentrations can be expressed as follows:

$$G_p≔S_{NHp,5}-4\le 0$$

(20)

$$G≔S_{NH,2}-4\le 0$$

(21)

The zeroth-order modifier is then expressed using steady-state measurements as:

$${ϵ}_k^G=G_p-G=S_{NHp,5}^{}-S_{NH,2}^{}$$

(22)

• Zeroth-order modifier for fast CA

Fast CA uses transient measurements denoted ${\left(·\right)}_{dyn}$, and the steady state of the plant is estimated using a dynamic model. For example, the steady-state value of $G_p$ can be estimated using Equation (11) as follows:

$${\widehat{G}}_p=G_p^{dyn}+\left(G-G_{}^{dyn}\right)$$

(23)

The zeroth-order modifier can then be evaluated as:

$${\widehat{ϵ}}_k^{}={\widehat{G}}_p-G_{}=G_p^{dyn}-G_{}^{dyn}=S_{NHp,5}^{dyn}-S_{NH,2}^{dyn}$$

(24)

3. Results and Discussion

This study assumes that a constant influent is treated by the WWTP. This influent is adopted from [24], and its parameters are given in

Table 1. Constant influent data.

Variables	${\mathit{S}}_{\mathit{S}}$	${\mathit{X}}_{\mathit{S}}$	${\mathit{X}}_{\mathit{B}\mathit{H}}$	${\mathit{S}}_{\mathit{N}\mathit{H}}$	${\mathit{S}}_{\mathit{N}\mathit{D}}$	${\mathit{X}}_{\mathit{N}\mathit{D}}$	${\mathit{Q}}_{\mathit{i}\mathit{n}}$
Values	69.5	202.32	28.17	31.56	6.95	10.59	18,446
Units	g COD·m−3	g COD·m−3	g COD·m−3	g N·m−3	g N·m−3	g N·m−3	m3/d

$S_O=0$ g COD·m⁻³, $X_{BA}=0$ g COD·m⁻³, $S_{NO}=0$ g N·m⁻³.

.

The control parameters for both the plant and the model are obtained using the Control Systems Toolbox in MATLAB©. First, single-input, single-output models are obtained in the form of transfer functions using step responses. This is followed by the design of two PI controllers. The control parameters for the plant and the model are given in

Table 2. Control parameters for the model and plant.

Control Parameters	${\mathit{K}}_{\mathit{C}1}$	${\mathit{\tau }}_{\mathit{I}1}$	${\mathit{K}}_{\mathit{C}2}$	${\mathit{\tau }}_{\mathit{I}2}$
Plant controllers	294	0.0069	4.7	0.066
Model controllers	122.8	0.0035	$2.6$	0.024

.

When the BSM1 control strategy presented in [24] is applied to the 47th-order plant for regulating the concentration of $S_{NOp,2}$ at 1 $\mathrm{g}\mathrm{N}/{\mathrm{m}}^3$ by manipulating $Q_a$ and that of $S_{Op,5}$ at 2 $\mathrm{g}\mathrm{COD}/{\mathrm{m}}^3$ by manipulating $K_{L}a$ in the fifth reactor, the plant OCI is 3674 $\mathrm{kWh}·{\mathrm{d}}^{-1}$. Note that $K_{L}a$ for the third and fourth reactors are fixed at 240 ${\mathrm{d}}^{-1}$, similar to those in BSM1.

This work uses the CasADi framework [41] to solve the numerical optimization problems. The variables of interest corresponding to the plant optimum are reported in

Table 3. Plant steady-state optimization results.

Variables	$\mathit{s}\mathit{p}{\mathit{S}}_{\mathit{N}\mathit{O}\mathit{p},2}^{ }$	$\mathit{s}\mathit{p}{\mathit{S}}_{\mathit{O}\mathit{p},5}^{ }$	${\mathit{S}}_{\mathit{N}\mathit{H}\mathit{p},5}^{ }$	${\mathit{Q}}_{\mathit{a}\mathit{p}}$	${({\mathit{K}}_{\mathit{L}}^{ }\mathit{a})}_{\mathit{p}}$	$\mathit{O}\mathit{C}{\mathit{I}}_{\mathit{p}}$
Plant SS values	0.5	0.81	4	34,460.9	153.81	2872.2
Units	g N·m−3	g COD·m−3	g N·m−3	m3/d	1/d	kWh/d

. These results will be used to validate the converged CA and fast CA solutions.

The variables of interest corresponding to the model optimum are reported in

Table 4. Model steady-state optimization results.

Variables	$\mathit{s}\mathit{p}{\mathit{S}}_{\mathit{N}\mathit{O},1}^{ }$	$\mathit{s}\mathit{p}{\mathit{S}}_{\mathit{O},2}^{ }$	${\mathit{S}}_{\mathit{N}\mathit{H},2}^{ }$	${\mathit{Q}}_{\mathit{a}}^{ }$	${\mathit{K}}_{\mathit{L}}^{ }\mathit{a}$	$\mathit{O}\mathit{C}{\mathit{I}}^{ }$
Model SS values	0.5	0.76	4	20,321.5	159.32	2913.6
Units	g N·m−3	g COD·m−3	g N·m−3	m3/d	1/d	kWh/d

. These setpoints are implemented on the plant to show the consequence of direct application of the model optimum to the plant.

As highlighted earlier, operating the plant using the model-based optimal setpoints results in suboptimal operation. In the worst case, it may even lead to infeasible plant operation due to plant–model mismatch. Such an instance is shown in

Table 5. Plant steady-state values at model optimum.

Variables	$\mathit{s}\mathit{p}{\mathit{S}}_{\mathit{N}\mathit{O}\mathit{p},2}^{ }$	$\mathit{s}\mathit{p}{\mathit{S}}_{\mathit{O}\mathit{p},5}^{ }$	${\mathit{S}}_{\mathit{N}\mathit{H}\mathit{p},5}^{ }$	${\mathit{Q}}_{\mathit{a}\mathit{p}}$	${({\mathit{K}}_{\mathit{L}}^{ }\mathit{a})}_{\mathit{p}}$	$\mathit{O}\mathit{C}{\mathit{I}}_{\mathit{p}}$
Plant SS values	0.5	0.76	4.65	38,636.3	151.53	2848.4
Units	g N·m−3	g COD·m−3	g N·m−3	m3/d	1/d	kWh/d

. It can be observed that $S_{NHp,5}>4$, thereby representing infeasible plant operation at the model optimum.

The characterization of the model optimal solution (Table 4) indicates that two constraints are active at the optimum, namely, $sp{S}_{NO,1}^{*}=0.5$ and $S_{NH,2}^{*}=4$. Since the optimization problem has two decision variables, optimality is completely determined by these two active constraints. Hence, updating $Q_a$ and $K_{L}a$ to meet these two optimality-determining quantities will drive the plant towards optimality.

CA scheme. Plant optimality can be verified by the CA results in Figure 5 and Figure 6. The red-dotted lines in Figure 5 represent the plant optimal values obtained by solving Equation (17), whereas in Figure 6 they represent the bounds on the nitrate, oxygen and ammonia concentrations. RTO is started from the steady-state values achieved by implementing the BSM1 control strategy to the plant. It can be observed that the CA scheme takes two iterations (around 120 days) to drive the plant close to optimality, where blue RTO lines overlap the red-dotted plant optimal lines. We notice here an important dichotomy, namely, two iterations are indicative of quick convergence, whereas 120 days is an unacceptable long time! This is due to the long time (60 days) taken by the plant to reach steady state following a system perturbation. This long settling time is intrinsic to WWTPs and cannot be changed with RTO schemes that rely on steady-state measurements. However, the practical message that can be taken from Figure 5 and Figure 6 is the following: CA updates the setpoints $sp{S}_{NO,1}$ and $sp{S}_{O,2}$ for the model (or $sp{S}_{NOp,2}$ and $sp{S}_{Op,5}$ for the plant and thus also $Q_{ap}$ and ${(K_L^{}a)}_p$ via the two PI controllers) in a direction that improves the OCI. Another message is that it would be quite useful in WWTP to be able to use transient measurements in order to avoid having to wait so long for steady state. This is precisely what fast CA does, as shown next.

Fast CA scheme. Plant optimality can be achieved in a single settling time using the fast CA approach as shown in Figure 7 and Figure 8. Blue-solid and red-dotted lines in these figures represent the same as they do in Figure 5 and Figure 6. The time between transient measurements, ${\tau }_{RTO}$, is chosen to be 1 day. It can be noticed that the plant reaches optimality within about 50 days, where the blue RTO lines overlap the red-dotted plant optimal lines. The value of the filter $K$ for the CA and the fast CA is chosen to be 0.6 and 0.5, respectively.

There are three constrained quantities in the optimization problem, namely, $sp{S}_{NO,1},$ $sp{S}_{O,2}$ and $S_{NH,2}$. The beauty of the adopted strategy is illustrated by the fact that, if two of the constraints are active at the optimum, the PI controllers can directly drive the plant to optimality in a single RTO iteration. Note, however, that feasibility is not ensured prior to convergence. This is clearly seen with CA in Figure 6, with the ammonia concentration exceeding its upper limit in the time period 10 < t < 80 d. In contrast, fast CA, which relies on most frequent measurements, is less aggressive and thus does not exhibit overshoot and constraint violation at all (see Figure 8).

4. Conclusions

The work has considered the real-time operation of WWTPs at minimal energy cost, while still satisfying the effluent discharge norms. It has been proposed to incorporate plant measurements in a model-based optimization framework in order to drive the plant to plant optimality. The two methods investigated here are CA and fast CA that both use measurements of the constrained quantities for bias correction in the optimization framework.

The procedure is iterative because the model correction at each iteration is only locally (and not globally) correct. By design, CA ensures feasibility upon convergence. In addition, it can drive the plant to optimality when the number of active constraints is equal to that of decision variables. This is precisely the case in this study of WWTPs, and it is observed that the pumping and aeration energy can be reduced by about 20% as compared to the BSM1 control strategy, while still satisfying the requirement on ammonia concentration in the effluent.

In this work, the plant is considered to treat an influent of constant flow and composition. Most of the WWTPs have a buffer system or equalization tank to feed a constant flow and composition influent to the WWTP. Nevertheless, the influent coming to the WWTP has time-varying flowrate and composition. For time-varying influent, the proposed methodology can, in principle, still be implemented to reach plant optimality. However, this would call for the online solution to a dynamic (rather than static) optimization problem, a much more challenging task, which could be the subject of further research.

Another topic of considerable interest is the possibility of using machine-learning techniques for the modeling, estimation and optimization of WWTPs. However, one has to keep in mind that these techniques need to rely on “representative measurements” in order to perform well. In the case of WWTPs, representative measurements probably need to include the various concentrations, which unfortunately are rarely available as a whole. Hence, whether machine-learning techniques are suited to tackle real WWTPs is still an open question. A few preliminary research results are available, see for example [42,43,44].