CHEMICAL ENGINEERING TRANSACTIONS
VOL. 88, 2021
A publication of
The Italian Association
of Chemical Engineering
Online at www.cetjournal.it
Guest Editors: Petar S. Varbanov, Yee Van Fan, Jiří J. Klemeš
Copyright © 2021, AIDIC Servizi S.r.l.
ISBN 978-88-95608-86-0; ISSN 2283-9216
Development of a Nonlinear Model Predictive
Controller of Flexible Structure for Greenhouse Gas
Emissions Minimisation in Pulp and Paper Wastewater
Treatment Processes
Marios Gkionisa, Athanasios I. Papadopoulosb, Wenhao Shenc, Panos Seferlisa,*
a
Department of Mechanical Engineering, Centre for Interdisciplinary Research and Innovation, Aristotle University of
Thessaloniki, 54124, Thessaloniki, Greece
b
Chemical Process and Energy Resources Institute, Centre for Research and Technology Hellas, 57001, Thermi-
Thessaloniki, Greece
c
State Key Laboratory of Pulp & Paper Engineering, South China University of Technology, Guangzhou, 510640,
PR China
seferlis@auth.gr
It is a well-established fact that activated sludge Wastewater Treatment Plants (WWTP) are placed among the
most significant contributors of Greenhouse Gas (GHG) emissions. The task of intelligent emission control of
such plants is a great challenge, mainly due to the highly nonlinear dynamics and associated uncertainty of its
constituent bioprocesses. The need for more efficient and environmentally conscious control is ever
increasing, due to the upsurge of industrial and urban water usage. In the current work, a nonlinear Model
Predictive Control (MPC) scheme that considers an activated sludge model is developed for pulp and paper
industry. The nonlinear MPC employs a flexible controller input structure through a process superstructure
formulation that enables the adaptive utilisation of the most suitable set of manipulated variables based on the
prevailing operating conditions. A generalised framework of this flexible controller type along with a nonlinear
dynamic model for the plant accounting for the simultaneous satisfaction of stringent control objectives on the
effluent quality and the GHG minimisation is developed. The proposed controller incorporates manipulated
variables for all admissible aeration intensities, internal flow recirculation streams, external carbon sources,
influent flow distribution in each bioprocessing unit, unit volume variations and recirculation flow percentage.
Accurate simulations indicate the achievement of a superior emission reduction performance compared to
state-of-the-art MPC implementations in WWTP. Every consideration of a manipulated variable set is
associated with respect to its contribution to the total GHG reduction for each mode of operation.
1. Introduction
Due to the constantly increasing demand in water usage for industrial and urban operations, WWTP are
targeted to operate under more stringent operating specifications. To successfully biometabolise the
wastewater and achieve acceptable effluent water quality, the WWTP system needs to be constantly
monitored and efficiently controlled. The biochemical reactions as contributors of GHG emissions have been
properly named as "On-Site" Emissions. Additionally, process equipment that are being deployed to control
the plant also contribute to a significant degree to the overall GHG emissions. Those are called the "Off-Site"
emissions, which result mainly from the aerator and pump operations.
The main challenges associated with the efficient control of WWTPs are:
• High variability of the operating conditions and the external inputs in the system (e.g., influent load
variation due to abrupt weather changes).
• Nonlinearity of the process system dynamics.
DOI: 10.3303/CET2188160
Paper Received: 6 June 2021; Revised: 7 August 2021; Accepted: 15 October 2021
Please cite this article as: Gkionis M., Papadopoulos A.I., Shen W., Seferlis P., 2021, Development of a Nonlinear Model Predictive Controller
of Flexible Structure for Greenhouse Gas Emissions Minimisation in Pulp and Paper Wastewater Treatment Processes, Chemical Engineering
Transactions, 88, 961-966 DOI:10.3303/CET2188160
961
• An abundance of unmeasured state variables that are important to describe the system behaviour (1).
• Uncertainty of model parameters, especially those related to the growth rate of microorganisms and
indicated in the work of Henze et al. (1987) and later enhanced by Alex et al. (2008).
Current state-of-the-art literature is limited to control approaches that use specific manipulated variables. A
common choice is the manipulation of the aeration intensity and recirculation flow in the aerator. Vega et al.
(2014) considered such a choice for the manipulated variables and employed a simplified plant configuration
consisted of two-units chamber. Zeng and Liu (2015) used a similar approach with the Benchmark Simulation
Model 1 (BSM1) plant configuration suggested by Alex et al. (2008). In Han and Qiao (2014) a similar input
controller structure is utilised in a MPC, which is complemented by a neural network structure that performs
online approximation of the system's states. In Stentoft et al. (2021), the on-off durations of aerator operation
are controlled. In Huang et al. (2020), the manipulated variables are the aeration intensities of chambers 3, 4
and 5 (Figure 1). In Sweetapple et al. (2014a), the input variables involved the internal recirculation flow,
wastage flow, aeration intensities of chambers 1 and 2, and the external carbon sources of chambers 1, 2 and
5. To the best of the authors’ knowledge, there is no published work that addresses WWTP control through a
flexible approach that allows the potential use of all available degrees of freedom in the process system.
There is, however, a study wherein the sensitivity of the available manipulated variables is investigated
(Sweetapple et al., 2014b), that identified the aeration intensities as the most influential factors with regards to
GHG emissions, a conclusion that agrees with observations made in the scope of this work. In the current
study, the deployment of all feasible manipulated variables is considered to enable the greatest possible
flexibility in the control scheme at every time instance. The main control objective is the minimisation of the
GHG emissions and operating costs, while maintaining the effluent quality within allowable limits.
Due to the constant increase in computing power, it is now possible to use complex optimisation algorithms in
industrial settings. Especially in the case of WWTPs, where time constants are in the order of magnitude of
approximately 15 min, it is realistic to efficiently use conventional non-linear optimisation algorithms in real-
time. In addition, high stiffness of the model differential equations, as WWTP systems encompass process
dynamics that are characterised by a large range of temporal scales, result in the need of efficient temporal
decomposition techniques. For instance, the control of Dissolved Oxygen (DO) requires a different time scale
than the control of ammonia and nitrate in the sludge (Weijers, 2000). There have been significant efforts
towards temporal decomposition of the WWTP dynamics over the past two decades. The proposed approach
aims also to provide a highly performing NMPC (Nonlinear MPC) leveraging the multiple control objectives.
2. Process system and manipulated variables’ configuration description
In this work, the objective function is consisted of a term for the generated GHG emissions, terms for the
manipulated variable variations, terms for state tracking errors between different levels of controller operation,
and a term for the Effluent Quality Index (EQI), which is a weighted average of effluent component loads that
have significant influence on the quality of effluent flow water and are usually considered to satisfy regional
regulations (Alex et al., 2008). As far as the relative importance between each optimisation term is concerned,
an extensive analysis in multi-objective optimal control is provided in Stentoft et al. (2021).
As an accelerator of nature’s processes, WWTP goal of operation is to achieve fast and effective elimination
of toxic nutrients in wastewater. Among the different biological treatment methods, Activated Sludge Process
(ASP) is usually selected. Inside the unit chambers of the WWTP (Figure 1), oxidation of the inserted biomass
and denitrification occur in aerated and anoxic conditions.
The plant configuration is described by two anoxic and three aerobic chambers, followed by a 10-layer non-
reactive secondary clarifier as shown in Figure 1. A generally preferred model that provides a good
compromise between dynamics accuracy and computational effort is the BSM1 model (Alex et al., 2008) and
Figure 1: Network diagram depicting the flow inputs, outputs and recirculations (arrows). Right-most
shape of Figure: notation for unit chamber k. Qk and Vask refer to the flow rate and volume of the k-th
chamber. Qin is the influent flow rate, which can split to the fQi fractions
962
is used in this work. The manipulated variables shown in Eq(1) involve the flow rates, aeration intensities and
unit chamber volumes.
uT = u
KLa
T , u
fQi
T , u
fQint
T , u
QEC
T , u
Vas
T , u
fQr
T , Qw (1)
where: u
vari
T = vari1,…,varinu
, var = {KLa, fQi, QEC, Vas, fQr }, ufQint
T = fQint1,1
,…, fQint1,nanox
,…, fQintnaer,nanox
KLa denotes the aerator-actuated aeration intensities or oxygen transfer (i.e., the oxygen flow rate that is
transferred to the microbial cultures to be biometabilised) applied to each unit chamber. Variables fQi, fQr , fQint
and Qw represent the influent flow rate, external recycle rate, internal recycle rates and wastage flow,
respectively (Figure 1). QEC,k represents the pump actuated supplementary biodegradable external carbon
source flow for the k-th unit chamber which facilitates denitrification. The time-dependent state vector (x(t) ≡ x
t
)
is defined as:
xT = x
1
T,…,x
nu
T , x
stll
T (Time index is omitted for simplicity)
The state vector (x
k
) that describes the k-th unit chamber is given below.
x
k
T = SIk , SSk , XIk , XSk , XBHk , XBAk , XPk , SOk , SNOk , SNHk , SNDk , XNDk , SALKk , 1 ≤ k ≤ nu (2)
In the above-given vector, the elements that start with “S” represent the soluble components, whereas those
with “X” represent the insoluble components. For example, SO represents the soluble oxygen concentration in
the flow of the k-th chamber. A similar structure is used for the influent flow x
in
(Figure 1). The state vector
(xstll) that describes the settler is given below.
x
stll
T = x
XTSS
T , x
SI
T , x
SS
T , x
SO
T , x
SNO
T , x
SNH
T , x
SND
T , x
SALK
T , x
solidsu
T
where: x
var
T = [var1,…,varnL ], var = {"XTSS", "SI", "SS", "SO", "SNO", "SNH", "SND", "SALK"}
and x
solidsu
T = XIu , XSu , XBHu , XBAu , XPu , XNDu , u: refers to underflow: Qu = Qf - Qe (Figure 1). XTSSl is the total
solid concentration in the i-th settler layer which is equal to 0.75·(XI+ XS+ XBH+ XBA+ XP).
BSM1 (Alex et al., 2008) is the basis for the formulations of the nonlinear differential state equations, derived
from mass flow balance equations for the five unit chambers, i.e. 13 equations for each chamber, and the 10-
layered non-reactive settler (the number of settler equations is equal to three related to the first three elements
of x
stll
T of Eq(2).The mass balance equations also include terms that represent the growth and decay of
biochemical components that take place inside the unit chamber due to reaction kinetics, as per Activated
Sludge Model 1 (ASM1) (Henze et al., 1987). Certain additions to the state equations of the BSM1 model are
considered, to account for the added flexibility of the controller, which include extra flow inputs and outputs to
unit chambers that arise from the potential flow recirculations (Figure 1). Each arrow in Figure 1 that points to
and from a chamber essentially represents a mass flow that has to be added in the state equations. No
modification is required for the settler equations that can be found in Alex et al. (2008). All the biochemical
components in the state vectors will be considered as known (i.e., perfectly estimated).
Table 1: Relationship between measurable variables and state variables in a unit chamber (table obtained
from Yin and Liu, 2018). T1 = XS+XI+XBA+XBH, conc.: concentration
Measurement State expression Measurement States’ expression
Oxygen Conc. SO Alkalinity Conc. S
Conc. of NH4
+ & NH3 S COD BOD+T1
Conc. of nitrate & nitrite oxygen S BOD SS+SI
Conc. of suspended solids T1+XP+XND
3. Methods
The nonlinear features of the BSM1 are fully included in the MPC, since linearisation of the model equations
results in significant performance deterioration (Vega et al., 2014). To calculate the GHG emissions a
combination of the BSM1, ASM1 and Bridle’s model is implemented, a choice that is also considered in Huang
et al. (2020). An analytic representation of the calculations can be found in Snip (2010).
963
The employed optimisation package is the open-source symbolic-numerical solver CasADi (Andersson et al.,
2019). The discretisation method used is the direct collocation method. The MPC related parameters are
given in Figure 2. The external disturbance is the influent flow load (i.e., flow rates) and component
concentrations of the components that correspond the chamber’s state vector x
k
in Eq(2). A premeasured
influent dataset in Alex et al. (2008) was utilised that corresponds to stormy weather conditions
(“STORMINFLUENT” dataset of Alex et al. (2008)). For this work, it is assumed that the influent flow signals
are correctly identified and forecasted. Similarly, to Huang et al. (2020), this work considers the following
sources of GHG emissions: endogenous decay, BOD removal, nitrification, denitrification, sludge disposal,
aeration, pumping and chemical usage. Other researchers such as Stentoft et al. (2021) consider as GHG
sources those attributed to nitrous oxide emissions and electricity consumption.
3.1 Temporal Decomposition
The temporal decomposition problem is not yet solved by any closed-loop method and has been tackled by
several authors in a heuristic manner. The timescale properties of the reaction kinetics model (ASM1) are not
very well understood (Weijers, 2000). A representation of the dynamics in different time scales is presented in
Weijers (2000). The proposed controller structure is a two-staged NMPC, as depicted in Figure 2.
The supervisory NMPC solves the state equations in steady state and minimises the GHG emissions. The
reactions can be assumed as static for large time scales and this module further utilises weather data
forecasts that are in essence generated by a sliding window of the initial premeasured raw weather data for a
time window of 0.5 d. For the scope of this work, it will be assumed that the influent signal is correctly
predicted, by applying the aforementioned smooth signal and incorporating measurement disturbances in the
future. The process level NMPC aims to implement the outcome of the supervisory level controller and further
satisfy the quality specification of the effluent stream.
3.2 Objective Functions
The generic form of the MPC objective function is given below.
Jt = xt+j,t - xreft+j,t T M xt+j,t - xreft+j,t + ut+j,tTQ ut+j,t+ fGHG xt+j,t, ut+j,t + fEQI xt+j,t, ut+j,t + feffl_Qual_viol(xt+j,t, ut+j,t)H-1
j=0
(3)
The supervisory NMPC does not contain a state-tracking error penalty, first term of Eq (3), whereas the
process level NMPC considers it, because it is dictated by the supervisory controller. Eq(4) represents the
penalty of exceeding the legislation-mandated effluent index (or discharge) limits, see Eq(5).
feffl_Qual_viol(x, u) =
0, x
eff i
≤ x
eff,lim i
weffi, xeff i
> x
eff,lim i
5
i=1
(4)
x
eff
T = Ntoteff ,CODeff,SNH,eff, XTSSeff ,BOD5eff < xeff,limT = [18, 120, 4, 30, 10] g/m3 (5)
4. Results and discussion
The simulations yielded GHG emissions (Figure 3b, wherein kg CO2 eq. represents the equivalent CO2 yield
from the various generated greenhouse gases) that are low for WWTP standards, given the specific discharge
limits in Eq(5)(4). The only manipulated variables that the supervisory MPC optimises are the unit chamber
volumes, since it is not realistic to variate them with the rate dictated by the process level MPC’s sampling rate
(Figure 2). By observing Figure 4 it becomes apparent that the controller indeed chooses to variate the
available manipulated variables at different time instances. For presentation simplicity, only the aeration
intensities for the chambers 4 and 5 are shown. The variables that experience the largest levels of utilisation
by the controller are the aeration intensities and the external carbon sources (Figure 4). For the influent flow
x
plant
Figure 2: Simplified outline of the controller - System architecture. The sampling time of each NMPC
controller is denoted as T and the prediction and control horizons as Np and Nc
x
ref
Supervisory
NMPC
Process Level NMPC Plant
T=15 min, Np = 5, Nc = 2 T = 0.5 days, Np = 4, Nc = 1
964
split it was assumed that fQi1
=1 and fQik
= 0 for k >1, since toxic influent compound concentrations such as SNH
need to be biodegraded so that the respective concentrations in the effluent flow remain below discharge
limits. By observing Figure 5d, it becomes apparent that the consideration of all aeration intensities as
manipulated variables leads to an increase in the aeration energy. This is expected since additional aeration
energy is required to variate the aeration intensities. The effluent quality related components are depicted in
Figure 5. Effluent indexes are mostly kept below limits and BOD5eff exceeds the limit at time periods that
correspond to the most severe influent load peaks. Compared to Huang et al. (2020), all values of effluent
quality indexes computed in the present work are lower, except for BOD5eff . The GHG emissions are reduced
to a factor of 2.6 in the present work compared to the ones in Huang et al. (2020) and to a factor of 5
compared to the ones in Kim et al. (2015). However, the effluent indexes in Kim et al. (2015) are maintained at
lower levels (specifically SNHeff , which is kept around 2 g/m
3). This observation indicates that minimisation of
GHG emissions can compromise the minimisation of the values in x
eff
(Eq.(5)). In Stentoft et al. (2021), the
lowest GHG emissions were around 65 kg/day, which are generated from N2O production and electricity
consumption and are approximately 15 times smaller than the corresponding ones computed in the present
work. However, direct comparison of GHG emissions in this case is difficult, since the plant configuration in
Stentoft et al. (2021) is very different in terms of number of unit chambers and corresponding volumes. As
stated in Fighir et al. (2019) expert user knowledge is required when comparing output data between different
WWTPs.
Figure 3: Optimal unit chamber volumes (a) and total GHG concentrations (in kg CO2-eq./h) (b)
Figure 4: (a) Optimal aeration intensities (units 4 and 5), (b) Optimal external carbon source flows (units 1 and
2)
Figure 5: Effluent Quality signals (a-c) and Aeration and Pumping Energies (d)
5. Conclusions
The consideration of a flexible WWTP model with regards to the manipulated variables provides superior
results compared to schemes that employ a limited number of manipulated variables and similar plant
configurations found in recent literature. The proposed approach is computationally efficient and allows for the
implementation of the control system under real-time conditions, through the utilisation of conventional non-
linear optimisation packages.
Nomenclature
BOD5eff – Effluent biochemical oxygen demand,
g/m3
CODeff – Effluent chemical oxygen demand, g/m
3
fQint – Internal recirculation fraction, m
3/d
fQr – External recirculation fraction
KLa – Aeration intensity, 1/d
naer – Number of aerated chambers
nanox – Number of anoxic chambers
NH3 – Ammonia
to
ta
l
G
H
G
[
k
g
/h
]
X
T
S
S
,e
ff
[
g
/m
3
]
N
to
t,
e
ff
[
g
/m
3
]
(a) (b)
(a) (b)
(a) (b)
(c) (d)
965
NH4
+ – Ammonium
Ntot – Total nitrogen concentration, g/m
3
nu – Number of unit chambers
QEC – External carbon source flow rate, m
3/d
Qw – Wastage flow rate, m
3/d
SALK – Alkalinity concentration, g/m
3
SI – Soluble inert organic matter, g/m3
SND – Soluble biodegradable organic nitrogen
concentration, g/m3
SNH – NH4
+ & NH3 nitrogen concentration, g/m
3
SNH,eff – Effluent SNH, g/m3
SON – Nitrate and nitrite nitrogen concentration,
g/m3
SO – Soluble oxygen concentration, g/m
3
SO,sat – Saturation SO
SS – Readily biodegradable substrate
concentration, g/m3
u – Vector of manipulated variables
Vas – Activated sludge unit chamber volume, m
3
XBA - Active autotrophic biomass concentration,
g/m3
XBH – Active heterotrophic biomass concentration,
g/m3
XI – Particulate inert solid matter concentration,
g/m3
XND – Particulate biodegradable organic nitrogen
concentration, g/m3
XP – Particulate products arising from biomass
decay concentration, g/m3
XS – Slowly biodegradable substrate
concentration, g/m3
XTSSeff – Total effluent solid concentration, g/m
3(∎)eff – Variable that refers to the effluent flow (∎) – Vector variable (∎)i – i-th element of vector
Acknowledgments
The work is supported by the European Union and Greek national funds through the Bilateral and Multi-
lateral Scientific and Technological Cooperation Greece – China (project code: T7DKI-00426).
References
Alex J., Benedetti L., Copp J., Gernaey K.V.,Jeppsson U., Nopens I., Pons M.N., Rieger L., Rosen C.,
Steyer J.P., Vanrolleghem P., Winkler S., 2008. Benchmark Simulation Model no. 1 (BSM1),
<https://www.iea.lth.se/publications/Reports/LTH-IEA-7229.pdf>, accessed 16/09/2021.
Andersson A.E.J., Gillis J., Horn G., Rawlings B.J., Diehl M., 2019, CasADi: a software framework for
nonlinear optimization and optimal control. Math. Prog. Comp., 11, 1–36.
Fighir D., Teodosiu C., Fiore S., 2019, Environmental and Energy Assessment of Municipal Wastewater
Treatment Plants in Italy and Romania: A Comparative Study, Water, 11(8), 1611.
Han H.G., Qiao J.F., 2014, Nonlinear Model-Predictive Control for Industrial Processes: An Application
to Wastewater Treatment Process, IEEE Transactions on Industrial Electronics, 61(4):1970-1982.
Henze M., Grady C.P.L. Jr., Gujer W., Marais G.v.R., Matsuo T., 1987, Activated Sludge Model No.1.
International Association on Water Pollution Research and Control - Scientific Technical Reports
No.1.
Huang F.N., Shen W.Η., Zhang X.W., Seferlis P., 2020, Impacts of dissolved oxygen control on different
greenhouse gas emission sources in wastewater treatment process, Journal of Cleaner Production,
274, 123233.
Kim D., Bowen J.D., Ozelkan E.C., 2015, Optimization of wastewater treatment plant operation for
greenhouse gas mitigation, Journal of Environmental Management, 163, 39-48.
Snip L., 2010, Quantifying the greenhouse gas emissions of wastewater treatment plants. PhD Thesis
Systems and Control, University of Wageningen, The Netherlands.
Stentoft P.A., Munk-Nielsen T., Møller J.K., Madsen H., Valverde-Pérez B., Mikkelsen P.S., Vezzaro L.,
2021, Prioritize effluent quality, operational costs or global warming? –Using predictive control of
wastewater aeration for flexible management of objectives in WRRFs. Water Research, 196,
116960.
Sweetapple C., Fu G.T., Butler D., 2014a, Identifying sensitive sources and key control handles for the
reduction of greenhouse gas emissions from wastewater treatment, Water Research, 62, 249-259.
Sweetapple C., Fu G.T., Butler D., 2014b, Multi-objective optimisation of wastewater treatment plant
control to reduce greenhouse gas emissions, Water Research, 55, 52-62.
Vega P., Revollar S., Francisco M., Martína J.M., 2014, Integration of set point optimization techniques
into nonlinear MPC for improving the operation of WWTPs, Computers & Chemical Eng, 68, 78-95.
Weijers, S., 2000, Modelling, Identification and Control of Activated Sludge Plants for Nitrogen Removal.
PhD Thesis, TU Eindhoven, The Netherlands.
Yin X.Y, Liu J.F., 2018, State estimation of wastewater treatment plants based on model approximation,
Computers and Chemical Engineering, 111, 79–91.
Zeng J., Liu J.F., 2015, Economic Model Predictive Control of Wastewater Treatment Processes,
Industrial & Engineering Chemistry Research, 54 (21), 5710–5721.
966