Denitrification Control in a Recirculating Aquaculture System—A Simulation Study

Abstract

The recirculating aquaculture system (RAS) is a land-based water treatment technology, which allows for farming aquatic organisms, such as fish, by reusing the water in the production (often less than 5%). This technology is based on the use of filters, either mechanical or biological, and can, in principle, be used for any species grown in aquaculture. Due to the low recirculation rate, ammonia accumulates in the system and must be converted into nitrate using nitrification reactors. Although less toxic for fish, nitrate can also be further reduced into nitrogen gas by the use of denitrification biofilters which may create several issues, such as incomplete denitrification, resulting in toxic substances, such as nitrite and nitric oxide, or a waste of carbon source in excess. Control of the added quantity of carbon source in the denitrification biofilter is then mandatory to keep nitrate/nitrite concentrations under toxic levels for fish and in accordance with local effluent regulations, and to reduce costs related to wasted organic carbon sources. This study therefore investigates the application of different control methodologies to a denitrification reactor in a RAS. To this end, a numerical simulator is built to predict the RAS behavior and to allow for the comparison of different control approaches, in the presence of changes in the operating conditions, such as fish density and biofilter removal efficiency. First, a classical proportional-integral-derivative (PID) controller was designed, based on an SIMC tuning method depending on the amount of ammonia excreted by fish. Then, linearizing and cascade controllers were considered as possible alternatives.

1. Introduction

Recirculating Aquaculture Systems (RAS) have been increasingly used due to the growth of the aquaculture industry [1]. They can be defined as systems where less than 10% of the total water volume is replaced per day and respond to the need of more intensive practices, growing environmental constraints on water consumption and effluent quality, as well as the possibility to supply fish in places where it would be otherwise difficult [2,3]. The main drawback of these systems is the accumulation of some toxic compounds, such as ammonia, when there is no proper treatment. Ammonia can be introduced into the system by fish excretion, 60–70% of the nitrogen consumption by fish being excreted as ammonia through the gills, or the degradation of uneaten feed [4]. The nitrogen removal process reduces the ammonia level using microorganisms that transform it into nitrate (nitrification). The latter, although much less toxic for fish, also accumulates in the system and needs to be removed to avoid system sustainability and environmental issues. Thus, another unit is often considered in order to convert nitrate into nitrogen gas (denitrification). The use of anaerobic denitrification to remove nitrate is not yet widely applied to commercial RAS due to its level of efficiency, its complexity and cost [5]. The complexity of RAS, created by its intrinsic closed-loop and the interactions between water treatment and fish grow-out, implies to build dynamic models to analyze the process behavior and to optimize it (configuration, size, fish, feed, flows etc) with respect to cost, robustness and water quality [6]. Indeed, manual control of denitrifying biofilters may easily lead to an unstable process performance and high cost [7]. Nitrite, which is toxic for fish, can accumulate in transient phases [8]. It is an intermediate compound in the denitrification process, a set of four reactions converting nitrate $N{O}_3^{-}$ to nitrogen $N_2$. Particularly, the lack of carbon source for complete denitrification could make this problem more severe. The absence of nitrite/nitrate/oxygen, however, if large quantities of organic matter exist, causes the production of sulfides (such as hydrogen sulfide) which are also extremely toxic [9]. Although nitrate control has been addressed many times regarding municipal wastewater treatment, its application in the context of RAS is not widespread [7,10,11,12]. A multivariable PID control, manipulating recirculation flowrates and aeration to control nitrate levels is proposed in [12] while an online control strategy/algorithm optimizing a denitrification biofilter through carbon dosage and backwash is shown in [7]. In [10], a simple feedback control strategy by methanol addition is used to minimize nitrite production and [11] proposes an output-feedback control scheme for a wastewater treatment biofilters in order to regulate nitrate and nitrite concentrations.

In this study, a global assessment of nitrate control using acetic acid (pH is considered regulated and constant) as manipulated variable is first achieved and two different control methodologies are considered: a classical PID control whose parameters are estimated by a SIMC tuning method [13], and a simple model-based linearizing control. PID control is widely used in process industries [14], and have been extensively used to control WWTPs. Unfortunately, PID controllers may present robustness limitations when applied to nonlinear systems with modeling uncertainties as well as measurement noise. They might indeed require tuning each time the operating conditions significantly vary. Linear controllers, on the other hand, take advantage of the availability of accurate nonlinear models often based on reaction networks and mass balances. However, as bioprocess models generally present uncertain kinetic structures, adaptive and/or robust solutions are required [15,16,17,18]. Taking advantage of both controller structures, a cascade of PI and linearizing control (accounting for possible model inaccuracies) is further proposed and compared with an adaptive linearizing strategy. The controller performances are tested with regard to typical industrial plant variations: changes in fish density, biofilter nitrate conversion efficiency, effluent guidelines regarding nitrate and water/acid pump malfunction in the denitrification reactor.

2. Problem Description

2.1. Industrial Plant

The industrial-scale RAS under consideration was composed of fish tanks, a physical filtration system and two biofilters (one for nitrification and another for denitrification). The reactor sizes and flowrates are presented in Figure 1.

The fish were reared in the tanks from where water flowed into a physical filtration system in order to remove particles. This water flow was then bypassed by a denitrification system where a carbon source, acetic acid, was injected. The nitrification filter was composed of 6 compartments, each one filled with fixed plastic support. Oxygen was added to the nitrification filter up to its saturation point and fresh water could also have been introduced in the system at a rate $F_{fresh}$.

In the following analyses, the RAS system is assumed to operate at a specific setpoint and model-based control strategies are designed to allow the user to choose specific operating conditions (i.e., to keep the current setpoint in steady-state or to chose another setpoint) and the closed-loop system to compensate possible disturbances. Regarding the simulated industrial plant, the nominal steady-state corresponds to a $6.4$ L/h acetic acid flow rate ($\overline{F_{Acid}}$) added to keep a constant nitrate concentration in the denitrification biofilter of 26.6 gN/m³. Water leaves the fish basins (and enters the denitrification biofilter) with oxygen, ammonium and nitrate concentrations of, respectively, 6.8 gCOD/m³, 0.5 gN/m³ and 98 gN/m³.

2.2. Model Description

Each block of the industrial RAS is described by the following models whose components and nomenclature are related to the Activated Sludge Models (ASM) from [19].

• Fish Tanks:
  The fish excretion rate is considered with assumed constant fish size and population. It can be translated into ASM variables using the waste matrix reported in [6]. The fish respiration is calculated based on the values reported in [20,21] (see

  Table 1. Fish waste production rates.

  Corresponding	Description	${\mathit{\vartheta }}_{\mathbf{waste}}$	Units
  component $\xi$
  $S_I$	Inert soluble organic material	994	gCOD/h
  $S_S$	Readily biodegradable organic material	0	gCOD/h
  $X_I$	Inert particulate organic material	994	gCOD/h
  $X_S$	Slowly biodegradable substrate	$1.2\times {10}^4$	gCOD/h
  $X_{BH}$	Active heterotrophic biomass	0	gCOD/h
  $X_{NS}$	Active ammonia oxidizing bacteria	0	gCOD/h
  $X_{NB}$	Active nitrite oxidizing bacteria	0	gCOD/h
  $X_P$	Part. products from biomass decay	$4.0\times {10}^3$	gCOD/h
  $S_0$	Dissolved oxygen	$-6.1\times {10}^3$a	gCOD/h
  $S_{NO3}$	Nitrate nitrogen	0	gN/h
  $S_{NH}$	Ammonium and ammonia nitrogen	1.5$\times {10}^3$	gN/h
  $S_{ND}$	Soluble biodegradable organic nitrogen	355	gN/h
  $X_{ND}$	Particulate biodegradable organic nitrogen	355	gN/h
  $S_{Alk}$	Alkalinity (as HCO3− equivalents)	0	gCOD/h
  $S_{CO2}$	Dissolved carbon dioxide	8.4$\times {10}^3$	gCO2/h
  $S_P$	Phosphorus	497	gP/h
  $S_{NO2}$	Nitrite concentration	0	gN/h

  ^a Negative since fish consume oxygen rather than expel it.

  ).
  The Fish basins are modeled as perfectly mixed tank reactors where no biological reaction occurs, applying mass balances as in:

  $$\frac{d\xi }{dt}=-D({\xi }_{inlet}-\xi )+\frac{{\vartheta }_{waste}}{V}$$

  (1)

  where variable ξ may be used to represent both soluble (S) and particulate (X) waste compounds and V is the fish tank volume. The corresponding waste production rates ${\vartheta }_{waste}$, all quantified in Table 1, are estimated for an assumed constant fish density of 12 kg/m³ and body weight of 8.5 kg. These rates can also be translated into ASM variables using the waste matrix reported in [6] while the fish respiration is calculated based on the values reported in [20,21] (see Table 1). The fish basin may therefore be modeled by the replication of the generic Equation (1) for each waste component of Table 1.
• Physical particle filter:
  90% of the particulate components in the inlet are filtered, leading to:

  $$X_{k,outlet}=0.1X_{k,inlet}$$

  (2)

  where $X_{k,outlet}$ and $X_{k,inlet}$ represent, respectively, the outlet and inlet concentrations of particulate component k.
• Denitrification filter:
  The denitrification filter consists of a moving bed perfectly-mixed tank reactor. Solubles that exist in the reactor bulk (index b) dissolve into the biofilm (index c) where the main part of the biological reactions occurs. The same ASM model, which is applied for nitrification, is used with the exception that ammonia is not considered as substrate for heterotrophic growth ($\frac{S_{NH}}{(K_{NHH}+S_{NH})}$ is not considered in ${\rho }_1$, ${\rho }_2$ and ${\rho }_3$). Acetic acid is added to this bioreactor (4) as a carbon source for denitrification, which can therefore be modeled as:

  $$\frac{d{S}_{i,b}}{dt}=D(S_{i,inlet}-S_{i,b})+\frac{{\dot{m}}_{i,Added}}{V_b}-\frac{k_{diff}(S_{i,b}-S_{i,c})}{V_b}+{\varphi }_{S_{i,b}}$$

  (3)

  $${\dot{m}}_{SS,Added}=F_{Acid}\phi$$

  (4)

  $$\begin{array}{c}\hfill \frac{d{X}_{k,b}}{dt}=D(X_{k,inlet}-X_{k,b})+\frac{K_d}{V_b}{X}_{k,c}+{\varphi }_{X_{k,b}}\end{array}$$

  (5)

  $$\begin{array}{cc}\hfill \frac{d{S}_{i,c}}{dt}& =\frac{k_{diff}}{V_{bf}}(S_{i,b}-S_{i,c})+{\varphi }_{S_{i,c}}\hfill \end{array}$$

  (6)

  $$\begin{array}{cc}\hfill \frac{d{X}_{k,c}}{dt}& =-\frac{K_d}{V_{bf}}{X}_{k,c}+{\varphi }_{X_{k,c}}\hfill \end{array}$$

  (7)

  where $F_{Acid}$ is the acetic acid flow rate in $L/h$, $\phi$ is the coefficient of conversion of acid into easily biodegradable organic matter in $gCOD/L$ and ${\varphi }_{\xi }$ is the production/consumption rate of component ξ, which can be calculated using stoichiometry—that is, the corresponding coefficient ${\upsilon }_{\xi ,w}$ (provided in Table S1, in Supplementary Materials) and rate ${\rho }_w$ (provided in Table S2, in Supplementary Materials), as in:

  $${\varphi }_{\xi }=\sum _{w=1}^{11}{\upsilon }_{\xi ,w}\times {\rho }_w$$

  (8)

  where w is an index denoting the process referenced in Tables S1 and S2 in Supplementary Materials.
• Nitrification filter:
  The nitrification filter is hydraulically modeled assuming six sequential perfectly mixed tank reactors filled with solid media. Biologically, a modified ASM1 model with the inclusion of two-step nitrification/denitrification is used (provided in Tables S1 and S2, in Supplementary Materials). This model has been validated using the data from the COST Simulation Benchmark in [22]. It has to be noticed that oxygen is also added to this bioreactor. Particles present in the inlet remain in the outlet. However, the particles already existing in each compartment ($rs$) do not move to others. The resulting mass balances lead to:

  $$\frac{d{S}_{i,rs}}{dt}=\frac{F}{V_b}(S_{i,rs-1}-S_{i,rs})+\frac{{\dot{m}}_{O2,Added}}{V_b}\pm {\varphi }_{S_{i,rs}}\frac{V_{biofilm}}{V_b}$$

  (9)

  $$\frac{d{X}_{k,rs}}{dt}={\varphi }_{X_{k,rs}}$$

  (10)

  $${\dot{m}}_{O2,Added}=k_{La}(S_{O,sat}-S_O)$$

  (11)

  $$S_{i,rs-1}=S_{i,inlet},atrs=1$$

  (12)

  where D represents the dilution rate, $k_{diff}$ is the diffusion coefficient in m³/h, $V_b$ is the empty volume of the reactor in m³, ${\varphi }_{S_{i,rs}}$ and ${\varphi }_{X_{k,rs}}$ are the consumption/production rates of solubles i and particulates k, calculated from (8), F is the volumetric flowrate in m³/h, $\dot{m}$ is the mass flowrate in g/h, $K_d$ is the detachment coefficient in h⁻¹·m³ and $k_{La}$ is the oxygen mass transfer coefficient in h⁻¹.

2.3. Plant Dynamics: A Preliminary Study

A first hydraulic study of the plant is achieved in order to assess the RAS inherent dynamics. To this end, a simulation of the concentration variation of an inert soluble component in the denitrification reactor following a production rate unitary step increase (in gCOD/(m³·h)) is shown in Figure 2 where, obviously, the plant takes up to 4000 h to reach the new steady state. Indeed, the step response can be approximated using a first order transfer function $G_1$ with time constant ${\tau }_C=744$ h—i.e., $5{\tau }_C\approx 4000$ h. As already mentioned, the plant is a semi closed system with 1% of regenerated water flowrate. This small inlet of water, along with the high considered volume, causes the very slow drift towards the new steady-state. This preliminary numerical analysis aims at demonstrating that the current RAS is very delicate to control due to the recirculating loop dynamics. In the upcoming sections, dedicated to control design, these slow dynamics are taken into account using adapted controller settings.

3. Control Implementation: Classical Approach

This section aims at first designing a classical PID controller to regulate the nitrate level in the denitrification reactor. From the corresponding model of differential Equations (3) and (4), acetic acid is obviously the sole manipulated variable and the following control law is proposed:

$$\begin{array}{cc}\hfill F_{Acid}\left(t\right)& =\overline{F_{Acid}}+K_C\left[e\left(t\right)+\frac{1}{{\tau }_I}{\int }_0^{t}e\left(t^{*}\right)d{t}^{*}+{\tau }_D\frac{de\left(t\right)}{dt}\right]\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill e\left(t\right)& =S_{NO3}^{*}-S_{NO3}\left(t\right)\hfill \end{array}$$

(14)

where $\overline{F_{Acid}}$ represents the initial acetic acid flowrate ensuring the industrial plant nominal steady-state ($\overline{F_{Acid}}=6.4$ L/h), $e\left(t\right)$ represents the measurement error, $K_C$ is the proportional gain, ${\tau }_I$ and ${\tau }_D$, respectively, are the integral and derivative time constants, $S_{NO3}^{*}$ is the desired nitrate concentration setpoint and $S_{NO3}\left(t\right)$ is the measured nitrate concentration. Expression (13) therefore aims at compensating any deviation from the chosen nominal steady-state trajectory which is assumed to be a priori knowledge from heuristics. Although (13) only presents three parameters, finding optimal settings without resorting to a systematic procedure may be complex. One way to design the gains is to first focus on the Proportional-Integral (PI) (i.e., ${\tau }_D=0$) part in order to set a satisfactory steady-state response (no steady-state error) and then consider the derivative part to improve the transient response [13,23].

In this study, the SIMC design method, as presented in [13], is chosen since it involves the transfer function of the system which can be easily computed considering a specific operating point and the corresponding step response, for instance, using the MATLAB model identification toolbox.

As shown in Figure 3, a one pole model already provides a very good approximation, and reads:

$$\begin{array}{cc}\hfill G\left(s\right)& =\frac{k{e}^{-\theta s}}{{\tau }_{1}s+1}=\frac{-32}{8.9\times {10}^{2}s+1}\hfill \end{array}$$

(15)

The numerator of (15) shows that the system presents a steady-state gain ($G\left(0\right)$) equal to $-32$ and no delay ($\theta =0$) while the denominator provides the system time constant (${\tau }_1$) which, as expected, is particularly large (890 h). Such slow systems, according to [13], may be approximated by a simple integrator:

$$\frac{k{e}^{-\theta s}}{{\tau }_{1}s+1}\approx \frac{k^{\prime }{e}^{-\theta s}}{s}$$

(16)

using the ratio between the gain $G\left(0\right)$ and the time constant ${\tau }_1$, denominated $k^{\prime }$, to determine PI-settings in a much shorter time frame (see Figure 4). $k^{\prime }$ is therefore calculated as follows:

$$k^{\prime }=\frac{\Delta y}{\Delta t\Delta u}=\frac{22.09-25.59}{(240-24)(6.762-6.44)}=-0.05$$

(17)

where $\Delta y$ represents the variation in nitrate concentration during $\Delta t$ hours, for a variation in acid flowrate $\Delta u$.

The SIMC method then proposes, for a general first-order process with delay, the PID tuning rules presented in

Table 2. PI parameter calculation using the SIMC method.

	Process	$\mathit{G}\left(\mathit{s}\right)$	${\mathit{K}}_{\mathit{C}}$	${\mathit{\tau }}_{\mathit{I}}$	${\mathit{\tau }}_{\mathit{D}}$
SIMC PID settings	First-order	$k^{\prime }\frac{e^{-\theta s}}{s}$	$\frac{1}{k^{\prime }}\frac{1}{{\tau }_C+\theta }$	$min\left\{4({\tau }_C+\theta )\right\}$	-
Industrial-scale plant model	First-order	$k^{\prime }=-0.05$${\tau }_1=8.9\times {10}^2$$\theta =0$	$\frac{1}{0.05}\frac{1}{{\tau }_C}$	$min\left\{4{\tau }_C\right\}$	-

, where ${\tau }_C$ represents the time constant of the desired first-order response of the system. The second row shows how this method is applied to the considered industrial-scale plant.

Considering y, the current nitrate concentration in the denitrification filter, $y_S$ the desired setpoint, and $\theta$ the time delay, Equation (18) shows the relation between these variables and ${\tau }_C$:

$$\begin{array}{c}\hfill {\left(\frac{Y\left(s\right)}{Y_s\left(s\right)}\right)}_{desired}=\frac{1}{{\tau }_{C}s+1}{e}^{-\theta s}\end{array}$$

(18)

The rate at which nitrate is produced/consumed in the system depends on either the quantity of ammonia released by fish (and consequent quantity converted into nitrate) or the amount of nitrate the denitrification biofilter is able to convert into nitrogen (gas). Regarding nitrate consumption, an increase in available carbon source induces biomass growth and nitrate consumption increases (assuming that the biofilter size is properly designed). However, nitrate production is entirely driven by the quantity of ammonia produced in the system. ${\tau }_C$ is then limited by the “maximum” production rate of nitrate and can be estimated using the quantity of ammonia produced by the fish. This can be further explained by Equation (18), considering $\theta =0$, represented in Figure 5. It has to be noticed that ${\tau }_C$ should always be greater than ${\frac{d\left(\frac{y}{ys}\right)}{dt}}_{max}$, that is, ${\tau }_C\ge {\tau }_{C,min}$ as shown in Figure 5, depending on the amount of ammonia excreted by fish.

The controller parameter calculation conditions, related to the case-study, are presented in

Table 3. Parameter values for the industrial case-study controller.

Parameter	Value	Parameter	Value
$r_{N{H}_4}$	$1.5\times {10}^3$ gNH4−N/h	${\tau }_C$	$\frac{0.63}{0.36}=1.8$ h
$V_{case-study}$	$4.2\times {10}^3$ m3	$K_{C,max}$	$-\frac{1}{0.05}\frac{1}{1.8}=-11$ L·m3/(h·gNO3−N)
$r_{NO{3}_{increase}}$	0.36 gNO3−N/h/m3	${\tau }_I$	$min\{0.3,4{\tau }_C\}=min\{8.9\times {10}^2,7.2\}=7.2$ h

. All the ammonia introduced in the system is considered as completely converted into nitrate.

As shown in Figure 6, where setpoint changes of +15/−15 gN/m³ are simulated, an anti-windup configuration (using Equation (19)) is recommended to avoid detrimental effects due to input saturation and the control law reads:

$$\begin{array}{cc}\hfill F_{Acid}\left(t\right)& =\overline{F_{Acid}}+K_C\left[e\left(t\right)+\frac{1}{{\tau }_I}{\int }_0^t(e\left(t^{*}\right)-e_{aw})d{t}^{*}\right]\hfill \\ \hfill e\left(t\right)& =S_{NO3}^{*}-S_{NO3}\left(t\right)\hfill \end{array}$$

(19)

with $e_{aw}$ being the anti-windup term calculated by the following equation:

$$\begin{array}{c}\hfill e_{aw}=\left\{\begin{array}{cc}{K}_{aw}(F_{acid,t}-F_{acid,max})\hfill & if{F}_{acid,t}>F_{acid,max}\hfill \\ K_{aw}(F_{acid,t}-F_{acid,min})\hfill & if{F}_{acid,t}<F_{acid,min}\hfill \\ 0\hfill & F_{acid,min}<F_{acid,t}<F_{acid,max}\hfill \end{array}\right.\end{array}$$

(20)

$e_{aw}$ is only greater than 0 when the controller output is greater than the allowed maximum value ($F_{acid,max}=200$ L/h) or lower than its minimum value ($F_{acid,min}=0$ L/h). When this happens, the control action is limited and $e_{aw}$ is calculated accordingly ($K_{aw}$ is the anti-windup gain). Figure 6 shows the effect of the anti-windup when a $K_{aw}$ value of −1 is used. When the nitrate setpoint is increased, the controller decreases the acid flowrate accordingly. However, nitrate concentration will continue to increase due to the ammonia conversion reaction which causes the integral term of the PI controller to go on increasing. When the nitrate concentration finally reaches the desired setpoint, the accumulated integral action entails overshoot. $K_{aw}$ values chosen between −0.1 and −10 are sufficient to solve the issue. Nevertheless, using a high $K_{aw}$ value tends to slow down setpoint tracking. A $K_{aw}$ value of $-1$ is ideally chosen and used in future simulations.

The design of the derivative gain is achieved following the analysis of the controller response to setpoint changes. Comparative results using different derivative gain values (driven by the choice of ${\tau }_D$) can be found in Figure 7. The input calculated by the PID therefore reads:

$$\begin{array}{cc}\hfill F_{Acid}& =\overline{F_{Acid}}+K_C\left[e\left(t\right)+\frac{1}{{\tau }_I}{\int }_0^t(e\left(t^{*}\right)-e_{aw})d{t}^{*}+{\tau }_D(e\left(t\right)-\eta )\right]\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \dot{\eta }& =(e\left(t\right)-\eta ){\tau }_D\hfill \end{array}$$

(22)

where the derivative term is approximated by a continuous high-pass filter using $\eta$ as intermediate variable. Including the derivative action leads to the new results of Figure 7, where oscillations are attenuated and the tracking is improved. A ${\tau }_D$ value of 1 provides the most satisfactory result and is therefore selected.

The resulting PID performance is assessed in Figure 8 where, following several setpoint changes and step disturbances in waste production (${\vartheta }_{S_{NH}}$), robust and fast tracking is achieved. Values of $-11.6$ L·m³/(g·h), 6.9 h, 1 h and $-1$ are, respectively, used for $K_C$, ${\tau }_I$,${\tau }_D$ and $K_{aw}$. Robustness analysis with respect to other operating condition disturbances is also presented in Figure 9. Variations in the available biomass concentration $X_{BH}$ and the bypass flowrate $F_{denitri}$ are identified as severe system disturbances. A sudden increase in biomass decay rate ($b_H$ term in Table S2 is increased 3 fold, from $2.6\times {10}^2$ to $7.8\times {10}^{-2}$ h⁻¹) is therefore simulated and indeed causes response oscillations, as shown in Figure 9, while step variations on the bypass flowrate are compensated very fast.

For the exact same perturbation in the biomass decay rate as in the previous test, the situation gets worse when adding measurement noise, as shown in Figure 10, where a Gaussian noise with zero mean and $10\%$ relative standard deviation is added to the measured variable ($S_{NO3}$). A sampling time of $0.01$ h is also used to match the probe measurement sampling process. The controller obviously faces tracking problems, sometimes reaching unacceptable nitrate concentration values even when trying to retune the controller proportional gain (Figure 10 shows the results for two values of the gain). To avoid this undesired effect, two solutions could be considered—they are, (a) low-pass filtering of the output signal or (b) replacing the PID by a PI, which is, in essence, less sensitive to measurement noise.

4. Linearizing Control

A linearizing control strategy, as illustrated in Figure 11, may also be proposed as an alternative solution.

The control law is established based on a simplified ASM1 model [19] (in order to limit the quantity of measured variables) and the following assumptions:

• Only the denitrification reactor is considered;
• This reactor only contains heterotrophic bacteria;
• Since the oxygen concentration is low, no nitrification can occur;
• Biomass is retained in the tank $(X_{BH,out}=0)$;
• Only three other components besides biomass are considered: oxygen ($S_O$), easily biodegradable organic matter ($S_S$) and nitrate ($S_{NO3}$);
• Only two reactions are assumed to occur:
  • Aerobic growth of heterotrophs:

    $$S_S+S_O\to X_{BH};$$

    (23)
  • Anoxic growth of heterotrophs on nitrate:

    $$S_S+S_{NO3}\to X_{BH};$$

    (24)
• Acid is added as a carbon source and $\phi$ represents the conversion of acid flowrate ($F_{acid}$, $L/h$) into easily biodegradable carbon source flowrate ($gCOD/h)$).

Applying mass balances to (23) and (24), the following simplified differential equation system is obtained:

$$\begin{array}{cc}\hfill \frac{d{S}_S}{dt}& =\frac{F_{Denitri}}{V_{Denitri}}{S}_{S,in}-\frac{(F_{Denitri}+F_{acid})}{V_{Denitri}}{S}_S+{\varphi }_{SS}+\frac{F_{acid}\phi }{V_{Denitri}}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \frac{d{S}_{NO3}}{dt}& =\frac{F_{Denitri}}{V_{Denitri}}{S}_{NO3,in}-\frac{(F_{Denitri}+F_{acid})}{V_{Denitri}}{S}_{NO3}+{\varphi }_{NO3}\hfill \end{array}$$

(26)

where ${\varphi }_{SS}$ and ${\varphi }_{NO3}$ are the easily biodegradable organic matter and nitrate reaction rates given by:

$$\begin{array}{cc}\hfill {\varphi }_{SS}& ={\upsilon }_{SS,aerob}{\rho }_1+{\upsilon }_{SS,anox}{\rho }_2\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\varphi }_{NO3}& ={\upsilon }_{NO3}{\rho }_2\hfill \end{array}$$

(28)

where ${\upsilon }_{SS,aerob}$, ${\upsilon }_{SS,anox}$ and ${\upsilon }_{NO3}$ are stoichiometric coefficients depending on biomass yield $Y_H$ as in:

$$\begin{array}{cc}\hfill {\upsilon }_{SS,aerob}& ={\upsilon }_{SS,anox}={\upsilon }_{SS}=-\frac{1}{Y_H}\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\upsilon }_{NO3}& =-\frac{1-Y_H}{2.86Y_H}\hfill \end{array}$$

(30)

and ${\rho }_1$ and ${\rho }_2$ are the aerobic and anoxic growth rates of heterotrophs:

$$\begin{array}{cc}\hfill {\rho }_1& ={\mu }_H\frac{S_O}{S_O+K_{OH}}\frac{S_S}{S_S+K_S}{X}_{BH}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\rho }_2& ={\eta }_g{\mu }_H\frac{S_S}{S_S+K_S}\frac{K_{OH}}{S_O+K_{OH}}\frac{S_{NO3}}{S_{NO3}+K_{NO}}{X}_{BH}\hfill \end{array}$$

(32)

Considering a first-order linear reference as in (33):

$$\frac{d{S}_{NO3}}{dt}={\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})$$

(33)

where $S_{NO3}^{*}$ is the desired nitrate concentration setpoint. Replacing Equation (33) in Equation (26):

$${\rho }_2=\frac{2.86Y_H}{1-Y_H}\left[\frac{F_{Denitri}}{V_{Denitri}}(S_{NO3,in}-S_{NO3})-\frac{F_{acid}}{V_{Denitri}}{S}_{NO3}-{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})\right]$$

(34)

In steady-state, $S_S$ variation tends to zero:

$$\frac{d{S}_S}{dt}\approx 0$$

(35)

and injecting (34) into (25), assuming (35), yields:

$$\begin{array}{cc}\hfill & \frac{F_{acid}}{V_{Denitri}}{S}_S-\frac{F_{acid}}{V_{Denitri}}\phi =\hfill \\ \hfill & \frac{F_{Denitri}}{V_{Denitri}}(S_{S,in}-S_S)-\frac{1}{Y_H}{\rho }_1-\frac{2.86}{1-Y_H}\left[\frac{F_{Denitri}}{V_{Denitri}}(S_{NO3,in}-S_{NO3})\right.\hfill \\ \hfill & \left.-\frac{F_{acid}}{V_{Denitri}}{S}_{NO3}-{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})\right]\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \iff F_{acid}\left[\frac{S_S-\frac{2.86}{1-Y_H}{S}_{NO3}-\phi }{V_{Denitri}}\right]& =\frac{F_{Denitri}}{V_{Denitri}}\left[S_{S,in}-S_S-\frac{2.86}{1-Y_H}(S_{NO3,in}-S_{NO3})\right]\hfill \\ \hfill & -\frac{1}{Y_H}{\rho }_1+\frac{2.86}{1-Y_H}{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})\hfill \end{array}$$

(37)

$$\begin{array}{c}\hfill \iff F_{acid}=\frac{\frac{F_{Denitri}}{V_{Denitri}}\left[S_{S,in}-S_S-\frac{2.86}{1-Y_H}(S_{NO3,in}-S_{NO3})\right]}{\frac{S_S-\frac{2.86}{1-Y_H}{S}_{NO3}-\phi }{V_{Denitri}}}\\ \hfill \frac{-\frac{1}{Y_H}{\rho }_1+\frac{2.86}{1-Y_H}{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})}{\frac{S_S-\frac{2.86}{1-Y_H}{S}_{NO3}-\phi }{V_{Denitri}}}\end{array}$$

(38)

As previously mentioned, an anti-windup term should be added to cancel input saturation effects:

$$\begin{array}{cc}\hfill F_{acid}& =\frac{F_{Denitri}(S_{NO3,in}-S_{NO3}+\frac{1}{2.86}{S}_{O,in})-\frac{1}{Y_H}{\rho }_1+{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})}{\frac{1-Y_H}{2.86}\phi +S_{NO3}}-\hfill \\ \hfill & \frac{{\lambda }_{lin,I}\int (S_{NO3}^{*}-S_{NO3}+e_{aw,lin})}{\frac{1-Y_H}{2.86}\phi +S_{NO3}}\hfill \end{array}$$

(39)

$$\begin{array}{c}\hfill e_{aw,lin}=\left\{\begin{array}{cc}{K}_{aw,lin}(F_{acid}-F_{acid,max})\hfill & if{F}_{acid}>F_{acid,max}\hfill \\ K_{aw,lin}(F_{acid}-F_{acid,min})\hfill & if{F}_{acid}<F_{acid,min}\hfill \\ 0\hfill & F_{acid,min}<F_{acid}<F_{acid,max}\hfill \end{array}\right.\end{array}$$

(40)

As shown in Figure 12, the linearizing controller provides good performances even in the presence of measurement noise with zero mean and $10\%$ relative standard deviation, and biomass concentration decrease. However, in an attempt to decrease the undesired oscillations, two additional strategies, which aim at being more robust to model uncertainties, are proposed in the following section, one combining both control structures in cascade and another estimating the biomass unknown dynamics.

4.1. Cascade Control

The cascade control methodology takes advantage of linearizing control to impose a first-order closed-loop dynamic in the inner loop, and a PI controller is used in the outer loop (see Figure 13). This cascade structure provides robustification, for instance, to model uncertainties, as follows:

$$\begin{array}{cc}\hfill F_{acid}& =\frac{F_{Denitri}(S_{NO3,in}-S_{NO3}+\frac{1}{2.86}{S}_{O,in})-\frac{1}{Y_H}{\rho }_{Lin}+{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})}{\frac{1-Y_H}{2.86}\phi +S_{NO3}}\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\rho }_1& ={\mu }_H\frac{S_{O,b}}{S_O+K_{OH}}\frac{S_S}{S_S+K_S}{X}_{BH}\hfill \end{array}$$

(42)

where

$$\begin{array}{cc}\hfill S_{NO3}^{*}\left(t\right)& =\overline{S_{NO3}^{*}}+K_{C,cascade}\left[e\left(t\right)+\frac{1}{{\tau }_{I,cascade}}{\int }_0^{t}e\left(t^{*}\right)d{t}^{*}\right]\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill e\left(t\right)& =S_{NO3}^{*}\left(t\right)-S_{NO3}\left(t\right)\hfill \end{array}$$

(44)

4.2. Adaptive Linearizing Control

Lumping the kinetic uncertainties into one single parameter $\widehat{\theta }$, Equation (38) can be re-written:

$$F_{acid}=\frac{F_{Denitri}(S_{NO3,in}-S_{NO3}+\frac{1}{2.86}{S}_{O,in})+\widehat{\theta }+{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})}{\frac{1-Y_H}{2.86}\phi +S_{NO3}}$$

(45)

with

$$\begin{array}{cc}\hfill \theta & =-\frac{1}{Y_H}{\rho }_1\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill \tilde{\theta }& =\theta -\widehat{\theta }\hfill \end{array}$$

(47)

where $\tilde{\theta }$ represents parameter mismatch $\tilde{\theta }=\theta -\widehat{\theta }$. An adaptive scheme, as described in [16,24], can be designed to develop an asymptotically stable control law, considering the following positive Lyapunov candidate function:

$$V_L\left(t\right)=\frac{1}{2}\left({\tilde{S}}_{NO3}^2+\frac{{\tilde{\theta }}^2}{\gamma }\right)$$

(48)

and its time derivative:

$$\dot{V_L}=\frac{d{\tilde{S}}_{NO3}}{dt}{\tilde{S}}_{NO3}+\tilde{\theta }\frac{d\tilde{\theta }}{dt}\frac{1}{\gamma }$$

(49)

where ${\tilde{S}}_{NO3}=S_{NO3}^{*}-S_{NO3}$, $\tilde{\theta }=\theta -\widehat{\theta }$ and $\gamma$ is a striclty positive parameter. All the terms in ${\rho }_{Lin}$ are assumed to be slowly varying so that $\frac{d\theta }{dt}=0$ and we get:

$$\frac{d\tilde{\theta }}{dt}=-\frac{d\widehat{\theta }}{dt}$$

(50)

Considering possible parameter mismatch, the error dynamics may be written:

$$\frac{d{\tilde{S}}_{NO3}}{dt}=-\lambda {\tilde{S}}_{NO3}-\tilde{\theta }$$

(51)

Replacing Equations (50) and (51) in (49), we obtain:

$$\dot{V_L}=-\lambda {\tilde{S}}_{NO3}^2-\tilde{\theta }{\tilde{S}}_{NO3}-\tilde{\theta }\frac{d\widehat{\theta }}{dt}\frac{1}{\gamma }$$

(52)

Considering Lyapunov stability theory, $\dot{V_L}$ should be negative and since the first term of (52) is, as long as the second and third terms cancel each other, the stability of the closed-loop is verified, as in:

$$\begin{array}{c}\hfill -\tilde{\theta }{\tilde{S}}_{NO3}-\tilde{\theta }\frac{d\widehat{\theta }}{dt}\frac{1}{\gamma }=0\end{array}$$

(53)

$$\begin{array}{c}\hfill \frac{d\widehat{\theta }}{dt}=-{\tilde{S}}_{NO3}\gamma \end{array}$$

(54)

The adaptive control law reads:

$$F_{acid}=\frac{F_{Denitri}(S_{NO3,in}-S_{NO3}+\frac{1}{2.86}{S}_{O,in})+\widehat{\theta }+{\lambda }_{lin}(S_{NO3}^{*}-S_{NO3})}{\frac{1-Y_H}{2.86}\phi +S_{NO3}}$$

(55)

with

$$\frac{d\widehat{\theta }}{dt}=-{\tilde{S}}_{NO3}\gamma$$

(56)

It should be noticed that the adaptive linearizing controller is equivalent to (38) where an integral action is added.

4.3. Numerical Results

Simulation results aim at assessing the controller performance following an increase in biomass (death rate) considering measurement white noise with zero mean and $10\%$ standard deviation on all measured variables ($S_{NO3}$, $S_S$, $S_O$ and $X_{BH}$), as shown in Figure 12. A sampling time of $0.01$ h is used to simulate the corresponding probe behavior. Regarding the linearizing controller, the parameters 0.5 h⁻¹, 6 L/h² m³ and $-10$ are, respectively, used for ${\lambda }_{lin}$, ${\lambda }_{lin,I}$ and $K_{aw,lin}$. The PI used in cascade is parameterized with $K_{C,cascade}=0.001$ and ${\tau }_{I,cascade}=0.5$ h (keeping ${\lambda }_{lin}=0.5$ h⁻¹) and the adaptive controller uses $\gamma =5$ (keeping ${\lambda }_{lin}=0.5$ h⁻¹). All the controllers perform reasonably well when noise is added and both the adaptive and cascade controllers are capable of attenuating the previously described oscillations.

In addition, comparative results between the linearizing controller with the cascade and the adaptive controller for several system disturbances are shown in Figure 14, Figure 15, Figure 16 and Figure 17.

It is noticeable from Figure 9 (left) and Figure 12 that all controllers oscillate following a perturbation in the biomass concentration. This is related to the operating point of the RAS and, by setting a new higher nitrate setpoint, it is possible to obtain a smoother response, as shown in Figure 18.

5. Conclusions and Perspectives

In this study, a control study of a recirculation aquaculture system is undertaken. Different control methodologies are tested in simulations for an industrial scale system. Both PID and linearizing controllers are capable of keeping the system stable when used with an anti-windup action. The PID controller gains can be calculated based on fish waste production (ammonia). The addition of a derivative action leads to a noticeable decrease in controller output oscillations. Nonetheless, the PID presents robustness issues when considering measurement noise. The linearizing controller is, however, able to overcome this problem and can also be further improved by adding an adaptive law to face possible model uncertainties. A cascade control strategy (combining PI and linearizing controllers) may also improve the control robustness to model uncertainties, compromising the tracking performance. The system has an operational region around which all the considered feedback strategies present oscillating responses. Operational conditions should then be changed to accommodate this phenomenon. The controllers show good robustness with respect to disturbances, such as the increase in ammonia production due to variations in the fish population. The control strategies discussed in this article are versatile and could be applied to other RAS configurations at different scales. Future perspectives could include the analysis of nitrite accumulation and the possibility to limit its concentration.