ON OPTIMIZATION OF SUBSTRATE REMOVAL IN A BIOREACTOR WITH WALL ATTACHED AND SUSPENDED BACTERIA

. We investigate the question of optimal substrate removal in a bioﬁlm reactor with concurrent suspended growth, both with respect to the amount of substrate removed and with respect to treatment process duration. The water to be treated is fed externally from a buﬀer vessel to the treatment reactor. In the two-objective optimal control problem, the ﬂow rate between the vessels is selected as the control variable. The treatment reactor is mod- elled by a system of three ordinary diﬀerential equations in which a two-point boundary value problem is embedded. The solution of the associated singular optimal control problem in the class of measurable functions is impractical to determine and infeasible to implement in real reactors. Instead, we solve the simpler problem to optimize reactor performance in the class of oﬀ-on functions, a choice that is motivated by the underlying biological process. These control functions start with an initial no-ﬂow period and then switch to a constant ﬂow rate until the buﬀer vessel is empty. We approximate the Pareto Front numerically and study the behaviour of the system and its dependence on reactor and initial data. Overall, the modest potential of control strategies to improve reactor performance is found to be primarily due to an initial transient period in which the bacteria have to adapt to the environmental conditions in the reactor, i.e. depends heavily on the initial state of the dynamic system. In applications, the initial state, however, is often unknown and therefore the eﬃciency of reactor optimization, compared to the uncontrolled system with constant ﬂow rate, is limited.


1.
Introduction. Biofilms are macrostructures of bacterial cells that are attached to an immersed surface and/or to each other in an aqueous environment. The first cells adhere to a surface, become sessile and start producing a gel-like matrix of extracellular polymeric substances [21]. Formation and growth of the macrostructure occurs through proliferation among the already sessile bacteria or through attachment of new bacteria from the surrounding aqueous phase. Eventually, bacteria may also leave the biofilm by means of various processes collectively known as detachment [19]. Bacteria exist in a suspended mode of growth in addition to biofilms [28].
Both biofilms and suspended bacteria are used in wastewater treatment for removal of unwanted compounds. In the activated sludge process bacteria are kept in suspension while substrates enter the reactor with the influent [34]. The treated water and bacterial sludge then leave the reactor with the effluent. Biofilm reactors provide a longer retention time for bacteria due to the surface attachment. These processes are particularly suitable for e.g. nitrogen and phosphorus removal, which are dominated by slow growing bacteria. Biofilm systems do not require sludge recirculation in order to maintain a high biomass density in the reactor, in contrast to the activated sludge process [29]. New technologies have combined biofilms and activated sludge in hybrid or integrated fixed-film activated sludge (IFAS) processes [17]. Such systems encourage growth of bacteria in both biofilm and suspended mode and are often implemented as an upgrade of an existing treatment plant. However, even pure biofilm reactors always contain a certain amount of suspended biomass which must be removed from the effluent before the treated water can be released to further treatment steps or into a receiving water body [15].
Mathematical models have long been used in wastewater treatment technologies to describe and analyze different processes [36]. Both activated sludge models and biofilm models have a well developed foundation and framework. Although suspended biomass is an integral part of biofilm reactor systems, it is generally neglected in biofilm models [26]. It is assumed that suspended biomass does not contribute to the biofilm system and need not be considered. We investigated the effects and role of suspended biomass in biofilm reactors in two previous papers. A one-dimensional single species single substrate model was presented in [22], where the substrate, biofilm and suspended bacteria exist in a continuous stirred tank reactor (CSTR). This model is a biofilm extension of the Freter model of microbial growth with wall attachment that has been studied in [3,4,10,18,32]. The model takes into account the exchange of biomass through attachment and detachment between the biofilm and the suspended biomass. The bacteria will either be completely washed out from the reactor or exist in both suspended and biofilm form. We studied the longterm dynamic behavior of the reactor, assuming an infinite constant flow with unlimited substrate. It was shown that suspended bacteria are relatively more efficient at substrate removal than the biofilm. The model was extended in a follow-up study in [23], where microbial complexity was introduced to represent the nitrification process for a multi-species multi-substrate case. Numerical simulations showed that suspended biomass need not be considered if the quantity of interest is the overall reactor performance with respect to ammonium removal, but it plays a significant role in the intermediate processes of nitrification.
In contrast to these earlier studies that focused on longterm reactor behavior with an infinite supply of substrates, we focus now on the question of optimizing the treatment of a finite amount of substrate in finite time. In the setup that we consider, the biofilm reactor is fed from a storage tank which is to be emptied. This leads to an optimal control problem with the reactor flow rate as control function. It should be chosen to minimize two objectives: (i) the amount of untreated substrate that is discharged from the reactor and (ii) the treatment time. These objectives are not co-operative and a suitable compromise must be established.
In the literature, several studies can be found in which optimal control and optimization are used to improve bioreactor and treatment efficiency. In [13], the substrate concentration in an activated sludge system with unlimited amount of incoming substrate and a given final time was minimized by choosing the optimal flow rate. A sequencing batch reactor (SBR), which has a cyclic filling and emptying pattern, was studied in [25]. The goal was to minimize the fill and reaction time in the reactor by varying the flow rate at which the reactor is filled, while satisfying the end condition at every cycle of achieving a minimum substrate concentration before the reactor is emptied. The unlimited substrate was supplied either in a constant manner or by a predefined function. A generic time and spatial optimal control was discussed in [31], where e.g. a maximization of the amount of biomass at a given final time in a well-mixed bioreactor by varying the flow was considered. Furthermore, in [11] the goal was to minimize the time it takes to decrease the substrate concentration in a polluted lake to a given value by manipulating the flow with which polluted water is pumped from the lake to a biological treatment reactor and clean water is pumped back into the lake. A chemostat model was used for the biological reactor.
In most of these studies it is assumed that substrate is unlimited and the goal is often to minimize process duration, i.e. a pure time-optimal problem. Our problem is not only a time-optimal problem, but also involves a second objective to increase removal efficiency in the reactor. The resulting optimal control problem is singular, leading to chattering control, and as such is difficult to treat with standard numerical methods for optimal control. Moreover, the reactor type that we consider, inspired by modern wastewater treatment technologies, contains both sessile biofilms and suspended bacteria, while the studies discussed above largely focus on reactors with suspended biomass only. Few optimization studies can be found in the literature, in which the focus is on a biofilm per se, but not the overall performance of an entire reactor. These studies are set in the context of time-dependent periodic dosing strategies for the eradication of biofilms with antibiotics [7,33].
In Section 2, we formulate the model as a dynamic system for biomass and substrate in the reactor, coupled with a two-point boundary value problem, the solution of which determines the substrate flux into the biofilm. We introduce the optimization problem for reactor performance with two objectives, for which we seek a compromise in the Edgeworth-Pareto sense. In Section 3, we remark on the problem in the traditional context of optimal control theory and on the difficulties to tackle the problem in this setting. In Section 4, based on biological considerations, we introduce a two-parametric class of control functions, and simplify the optimal control problem as a nonlinear vector optimization problem, for which we approximate the Pareto front numerically. Finally, in Section 5, we discuss briefly two alternative approaches, which, however, for computational reasons require a significant simplification of the underlying reactor model. Results and conclusion are summarized in Section 6.
2.1. Mathematical model. The single species single substrate model from [22] was constructed through a coupling between a CSTR mass balance and a traditional one-dimensional biofilm model. The bacteria are present in the reactor in suspended and biofilm form. Microbial growth is limited by a single substrate, which is continuously supplied to the reactor. Bacterial cells can detach from the biofilm and become suspended, which in turn can (re-)attach to the biofilm. It is assumed that the biofilm covers the substratum in the reactor uniformly and that the substrate diffuses into the biofilm creating concentration gradients due to diffusion and reaction. The flow, with which substrate is supplied, also removes suspended biomass and unconsumed substrate with the effluent.
Governing equations. The model is formulated with respect to the three dependent variables substrate concentration in the aqueous phase S [g/m 3 ], suspended biomass u [g] and biofilm thickness λ [m], and reads:  [3,4,10,18,32]. The biofilm detachment function d(λ) that we use is defined as where E [1/md] is the erosion parameter. This is the most frequently used detachment rate in the biofilm modeling literature. Using the same notation as in [22] the substrate net flux into the biofilm is defined through the function j(λ, S) [g/m 4 ] where C(z) [g/m 3 ] is the substrate concentration in the biofilm at depth 0 < z ≤ λ. The specific growth rates for suspended and biofilm biomass are µ u (S) [1/d] and µ λ (C(z)) [1/d] respectively, defined as where µ max u , µ max λ [1/d] are the maximum growth rates and K u , K λ [g/m 3 ] the Monod half-saturation concentrations. The local substrate concentration C(z), 0 ≤ z ≤ λ, is found as a solution to the quasi-steady state diffusion reaction equation posed as a boundary value problem Here we made use of the standard time scale argument that diffusion and consumption in the biofilm is fast compared to microbial growth.
Properties of the flux function. We observe that the function j(λ, S) in (5) depends on S through the boundary condition in (7). As a quantity that is derived from the solution of a two-point boundary value problem, j(λ, S) is difficult to evaluate. However, many of its qualitative properties are known [2,20,22]: It is well-defined and twice differentiable for λ > 0, S > 0. We have j(·, 0) = j(0, ·) = 0. Moreover, j(λ, S) is monotonically increasing in both its arguments and concave in S. Using standard comparison theorems, the flux function j(λ, S) can be estimated from below and above by The wastewater engineering literature contains several empirical formulae that approximate j(λ, S). These have been developed based on experimental data but usually contain additional parameters that are reactor specific. Using two very different approaches [2,20] one obtain the approximation that only involves parameters of the underlying two-point boundary value problem (7). Numerically it was verified in [2] that this approximation appears to work reasonably well for an extended range of realistic biofilm parameters in some realistic ranges for λ and S, however, primarily for thin biofilms. Below we will normally use the full flux calculation according to (5) in our numerical simulations and resort to the flux approximation (8) only where computational costs prevent us from using (5) in Sections 5.1 and 5.2.
With the above properties of j(λ, S), (1)-(3) is formally an ordinary differential equation, albeit with a difficult to evaluate right hand side. In [22] its longterm behavior was analyzed for a constant flow rate Q, where also easy to evaluate conditions for the stability and instability of the washout equilibrium were derived.
Solutions in the case of non-constant flow rate. In the context of optimal control problems, the flow rate Q will be the control parameter, i.e. is no longer constant. The natural setting is to search for control functions Q in the class of measurable functions. Then, (1) -(3) becomes a non-autonomous system of differential equations. For non-negative S, u, λ and measurable Q, the right hand side of (1) -(3) is continuous in the dependent variable for fixed t, and it is measurable in t for fixed S, u, λ. Thus, solutions of (1)-(3) are understood in the generalized sense of Carathéodory [6]. Many of the standard results on ODEs carry over to Carathéodory differential equations accordingly [35]. The initial value problem of (1)-(3) with positive initial data possesses a unique solution that remains nonnegative almost everywhere. Moreover, by forming a linear combination of (1), (2), (3), we obtain for non-negative initial data the upper estimate almost everywhere. Thus, biomass production is bounded by the amount of substrate supplied. Furthermore, if the initial data are positive, there is no interval (t , t ) such that either u(t) = 0, λ(t) > 0 or u(t) > 0, λ(t) = 0 for t < t < t . These results carry over directly from [22] to the setting of Carathéodory solutions, using is pumped to the treatment rea.ctor of volume V at a flow rate Q, with treated water leaving the rea.ctor at the same rate. T he t reatment reactor is aerated from the bottom and contains bacteria in the form of biofilms and suspended biomass, which ensure that the process is not oxygen limited but that substrate S is the only growth limiting factor .
2.2. Optimiza t io n p roble m . R eactor set u p . T he analysis and investigation of t he basic model (1)-(7) in [22] assumed a constant flow rate and essentially a n infinite amount of substrate to be treated; t he focus of t he investigation was on longterm behavior, i.e the asymptotics of the model for t ---+ oo. From an engineering perspective, it is useful to investigate the question of treating a finite amount of substrate in finite time. More specifically, one is interested in the problem, for example, to degrade a given amount of substrate as much as possible in as short a t ime as possible.
To this end we consider a t reatment reactor , described by (1)-(7) that is fed from a buffer reactor that contains the volume Vb(t) [m 3 ] of water with a bulk substrate concentration So [g/m 3 ], see Figure 1. T he emptying of the buffer reactor is described by the simple mass balance where Q(t) is a non-negative function such t hat J;[ Q(t)dt = Vb(O) for some T > 0 a nd Q(t) = 0 for t > T. T his restriction guarantees t hat t he buffer reactor is emptied in finite t reatment t ime T and that negative buffer volumes are avoided.
Vect o r optimizatio n p roble m . T he task at hand is now formulated as the vector optimization problem Qen T (11) where Ω is the set of admissible flow rates Q : [0, T max ] → R + 0 . This set will be made more precise below. T max is the maximum duration of the process we are prepared to tolerate.
The two objectives are not cooperative. In order to minimize the first objective, i.e. the amount of substrate in the effluent, one would expect that large process durations be advantageous, so the microbes have sufficient time to degrade almost all the substrate. This, however, obviously counters the objective to minimize process time.
The trade-off between both goals that we seek is a compromise in the sense of Edgeworth-Pareto optimality, with the positive cone as order cone [12,16]. We denote the two objectives by and the vector of objectives by Z := (z 1 , z 2 ) T . An admissible flow rate Q * (t) is said to be Edgeworth-Pareto optimal if a further improvement of one objective is only possible at the expense of increasing the other. We call then Z * = (z 1 (Q * ), z 2 (Q * )) T an EP-minimum. We must not expect that the EP-minimum is unique [12,16]. Rather, we expect that there might be infinitely many such optima. The set of these optima can be graphically represented in the z 1 -z 2 plane. They form the graph of a monotonically decreasing function, the so-called Pareto Front. Points below the graph are not attainable with an admissible flow rate, points above the graph are not optimal. To compute the Edgeworth-Pareto minima, the vector optimization problem is converted into a family of scalar optimization problems. Scalarization with linear functionals. A standard approach to compute Edgeworth Pareto optima is to scalarize the vector objective function Z with a monotonically increasing functional F : R 2 → R and to solve the scalar problem min Q∈Ω F (Z(Q)).
The solution of this scalar optimization problem is an Edgeworth-Pareto optimal solution of the vector optimization problem (11). Using linear functionals, the objective function of (12) becomes Thus, we arrive at the following optimal control problem in Bolza form: Here ) is a constant that is introduced to rescale the problem such that both objectives have the same units. This allows to replace w 1 and w 2 in (13) by a single dimensionless weight w with 0 < w < 1. The EP-optimum that is found this way as a solution of (14) depends on the choice of w. The Pareto Front is the union of these points. In the numerical realization, (14) is solved for a finite number of weights that is large enough to approximate the Pareto Front well.

Modified Polak
Algorithm. An alternative approach to determining the Pareto Front numerically, which is particularly useful and easy to implement for our problem due to the special form of z 2 , is the Modified Polak Algorithm [16].
The interval (0, T max ] is divided into a finite number of n subintervals by introducing is solved, where the sets are the sets of admissible functions such that the reactor is emptied in finite time T i . The Pareto Front can be approximated by interpolation after removing those segments of the curve that perturb monotonicity, if any. In the next sections we will make use of both approaches to discuss the solutions of (11), depending on the situation.
3. Remarks on the optimal control problem. Admissible functions in optimal control theory are often defined in the class of measurable functions. Most standard theorems and results in optimal control theory are derived for these types of functions. Furthermore, realistically there exists an upper bound Q max on the flow rate Q due to the pump capacity between the buffer reactor and the treatment reactor. We therefore proceed by reformulating our problem and studying the set of measurable, bounded flow rates.
The constraint on Q, which enforces emptying of the reactor, i.e.
T 0 Qdt = V b (0), can be restated in the form of a differential equation for the buffer volume V b . Thus, the system of equations (1)-(3) now readṡ with initial conditions and an end condition for the optimal control problem in Lagrange form where T min = V b (0)/Q max is the minimum possible time and Ω is the set of admissible control functions. Let It can be shown by using standard existence results [9] that the optimal control problem (22) with the set of admissible functions Ω, subject to (16)- (19) with initial conditions (20) and end condition (21) possesses an optimal solution Q * . The optimal solution may be characterized with Pontryagin's minimum principle [9,27]. The Hamiltonian is where λ 1 (t), . . . , λ 4 (t) are the adjoint variables given by with the transversality conditions Pontryagin's minimum principle states that for the optimal solution Q * , where S * , u * , λ * , V * b are the optimal state trajectories and λ * 1 , λ * 2 , λ * 3 , λ * 4 the optimal adjoint trajectories. Furthermore, since (22) is a free final time problem and the Hamiltonian does not explicitly depend on time, at the optimal solution the Hamiltonian must satisfy for all t ∈ [0, T * ]. The optimal control trajectory can be obtained by minimizing the Hamiltonian with respect to Q. To this end, the Hamiltonian can be rewritten as , this is now the switching function. H is linear in Q, implying a bang-bang control as the optimal solution, switching between 0 and Q max . The switching times can be determined with the signs of the switching function h.
Optimal control problems like (22) that are linear in the control variable, with bounded control sets, are known to lead to chattering [37]. Such controls are characterized by an infinite number of switches between the lower and upper bound in a finite amount of time. This type of control regime is, however, not appropriate from a practical point of view, wherefore this control set will not be considered here.
T he solutions to (22) in the class of measurable and bounded admissible functions (23) depend heavily on Qmax , the upper bound for the flow rate. Giving up this restriction, we permit Q to be unbounded, i.e. let n be the set of measurable, non-negative flow rates Q defined on [O,Tmaxl· For the optimal control problem (14), i.e., a minimum problem with a linear integrand, and with this definition of n, Fleming and Rishel [9] indicate that impulse control solutions often are found.
T hese controls are characterized by impulses of the Dirac delta function type and cause instantaneous jumps between different values in the state variables. Such a behavior is, for practical purposes, not appropriate and will, therefore, not be pursued further. 4. A finite-d ime nsion a l optimizatio n p roble m . 4.1. O ff-o n con trol function s . T he rate with which substrate can be degraded depends on both, the substrate concentration and the amount of bacterial biomass in the system. In a typical application, the bacterial populations, wall atta.ched and suspended, a re initially not fully developed, but must adjust to the environmental conditions in the reactor. Since the overall degradation rate scales with amount of biomass, in this initial period the reactor performance is limited by small population sizes. If the flow rate Q in this initial period is too high, much more substrate is supplied than can be utilized and more substrate might be washed out of the rea.ctor than degraded. Therefore, a first and simple candidate class for admissible flow rates are off-on functions. T he flow rate is initially turned off for a transient period during which the bacterial population can establish itself. Later it is switched on to a constant value until the buffer reactor is empty. With this definition the admissible functions are the functions { 0, An example of such a function is shown in Figure 2. Since for practical purposes the buffer rea.ctor cannot be emptied instantaneously, it is appropriate to introduce the minimum required treatment t ime T min . E very pair (Tswitch, T) with T min ≤ T switch < T ≤ T max defines uniquely one such function Q(t) and, conversely, every such function is defined by a unique pair (T switch , T ). Therefore, the set of admissible functions is a two-parametric family of functions and finding the optimal Q * ∈ Ω is equivalent to finding the corresponding pair (T * switch , T * ). With this choice of admissible controls, the optimal control problem (11) reduces to a vector optimization problem in two unknown parameters. We use the Modified Polak Method to scalarize the problem: we solve for selected T i , i = 0, ..., n with T min = T 0 and T max = T n , and T 0 < T 1 < T 2 < ... < T n the problem with Q according to (33). Thus, the problem reduces to finding the solutions of n+1 scalar optimization problems with one argument for every given T i . To compute the objective function requires to solve numerically (1)- (7) and to evaluate the integral (34).

4.2.
Computational setup. The model parameters for the numerical simulations are found in Table 1, where, as in [22], the parameters were taken from Benchmark Problem 1 of the International Water Association's Taskgroup on Biofilm Modelling [36], which models an aerated reactor with carbon as the only substrate, and from [22]. As in [22], we use the same parameter values for both biofilm and suspended biomass, implying that suspended biomass is treated as unattached flocs of microbial communities. All simulations, unless otherwise specified, are initiated with the following initial conditions: We consider a small laboratory scale buffer reactor with volume V b = 0.05m 3 = 50l. The minimum and maximum treatment periods were chosen as T min = T 0 = 0.1days = 2.4hours and T max = T n = 50days.
The treatment times T i were selected on a non-equidistant grid with n = 62 points. The grid points are chosen denser for small T than for large T . All simulations were executed in MATLAB R2012a using the built-in solvers ode15s (a variable order, multi-step solver based on the numerical differentiation formulas) for the solution of the system of ordinary differential equations (1)-(3), bvp4c (a finite difference code using the three-stage Lobatto IIIa formula, which is a collocation formula [24]) for the diffusion-reaction equation (7) and fminbnd, a derivative-free method based on golden section search and parabolic interpolation [24], for the minimization problem (34).

4.3.
Results and discussion. The choice of using a constant flow rate Q requires no further adjustment after initiation until the buffer reactor is emptied. It will, therefore, serve as a baseline for comparison with other flow rates to measure efficiency of a control strategy. Solution of the vector optimization problem. The Pareto Front for the optimal control problem (11), as determined by the Modified Polak Algorithm is shown in Figure 3a. This curve is monotonically decreasing from (max z 1 , min z 2 ) to (min z 1 , max z 2 ) with respect to z 2 : decreasing z 2 further will only be possible at the expense of increasing z 1 and vice versa. Increasing the treatment time T allows a   Each point on the Pareto curve is an EP minimum of (11) and represents a compromise between treatment time and reactor performance. In general it depends on user preferences which of these compromises to accept. One might argue that of all the EP minima a distinguished one is the point on the Pareto front that is closest in some suitable metric to the (unachievable) ideal target (T min , 0).
We found the solution of this goal programming problem [12,16] by normalizing the z 1 -and z 2 -axis, computing the Euclidean distance and minimizing over T . The closest point on the Pareto Front was found at T = 5.09285 days with z 1 = 0.14499. We plot the optimal performance of the off-on-controlled reactors as a function of treatment time T in Figure 3b as blue asterisks, together with the performance of the uncontrolled reactor, i.e. the reactor with constant flow rate Q = V b (0)/T . Both curves are very close to each other, with hardly any difference visible. When comparing the off-on functions with the constant flow rate Q we see a small improvement in reactor performance with off-on functions over constant Q in Figure  3c. The relative improvement ranges between approximately 0.5% and 3.5%. For short treatment times the improvement remains small and begins increasing when the treatment time is increased, reaching a maximum around T = 10 days, after which it slowly decreases.   In Figure 5 we plot the time evolutions of the three state variables S, u , A for three selected treatment times T = 2, 5, 14 days. The time-dependent functions are plotted on a normalized time axis. Initially, while the flow is turned off, the substrate in the treatment rea.ctor is degraded by suspended and wall attached biomass, which increases. The flow is switched on rather early, which is clearly visible in Figure 5a. T his brings fresh substrate into the rea.ctor. Hence S increases. Simultaneously, the flow washes out some of the suspended biomass, whence u decreases. The biofilm thickness keeps increasing throughout. Eventually all three state variables level off. T he treatment time T has an effect on the state that is attained by the reactor. T he longer the treatment period, the smaller is the final substrate concentration. Smaller bulk substrate concentrations imply that the biofilm that can be sustained will be thinner, as observed in Figure 5a. While the effect of treatment time T on the suspended biomass is clearly visible, it is not as straightforward to summarize. For one, shorter process duration implies higher substrate concentrations and therefore faster growth of suspended biomass. On the other hand, shorter process duration implies higher flow rates and therefore higher washout of suspended biomass. The latter phenomenon is also reflected by the minimum to which u dips after the flow is turned on. For T > 5, the value that u levels off at does not change notably.
Complementing the above data, we plot in Figure 6   In Figure 6 the maxima of suspended biomass and biofilm thickness in dependence on the length of the t reatment period are most likely related to substrate depletion: For small T the biomass at the end of the t reatment time increases because the more time the bacteria are given the more substrate they can remove and hence grow. However, for longer process durations, the system becomes eventually substrate depleted which explains the monotonic decrease in biomass for large T due to starvation and washout.
D ep e nde n cy of r eactor p erformance on initia l data. T he definition of off-on functions in 4.1 was motivated by an initial adaptation period for the bac-  thickness were varied. A Pareto Front was approximated for ea ch variat ion (data not shown), and for comparison, the reactor performance was also calculated for the constant flow rate. Initial suspended biomass uo was 0.005, 0.02, 0.1 and 0.5 g, compared with the standard 0.05 g; initial biofilm thickness .Xo was 10, 50, 200 a nd 500 µm, compared with the standard 100 µm. Figure 7 shows the relative improvement in objective z 1 when using optimal off-on flow rate functions compared with constant Q. T he dashed lines denote the standard values for uo and .Xo that were used throughout the previous simulations. Only one initial value was varied in each instance. We see in Figure 7a that the largest improvement is obtained for the smallest amount of initial suspended biomass. The improvement increases with increasing treatment t ime and the overall improvement ranges between 0.53 a nd 53 . The less suspended biomass is initially present, the more pronounced is the effect of the flow control compared to the uncontrolled system for process durations T > 2. This statement does not hold for smaller process durations, where the washout rate after the flow is turned on might be substantial. Increasing process duration for T > 3 leads to an increasing improvement by control, up to T ~ 10 when the rea.ctor performance levels off. T he results for varying biofilm biomass a re plotted in Figure 7b. T he overall improvement ranges approximately between 0.53 and 123 . For a very thin initial biofilm thickness .X(O) = 10µ,m the largest relative improvement can be found for a ll T. T he improvement increases monotonically with respect to T for the three thickest initial biofilms. However, for the thinner biofilm with .X(O) = 10µ,m we see first an increase to a local maximum around T = 3days after which a decrease occurs.
In both simulation experiments run here, the relative gain by control with off-on functions is biggest if the initial amount of biomass is smallest. This confirms the reasoning that led us to choose this type of control functions: T he initial off-period is required for the bacterial populations to attain a sufficiently large population size to be efficient in substrate degradation. : Sensitiv ity w it h respect to b iofilm g r owth area. In the technical realization the wastewater treatment process can be accelerated by providing a larger substratum area A for biofilm to grow on. In Moving Bed Biofilm Reactors, this is achieved by adding suspended biofilm carrier chips into the completely mixed treatment vessel. To investigate the effect of control on such an improved reactor, we repeat the numerica l control experiment for different numbers c of such carriers. T he initial data were set to the default values again.
T he results of this simulation experiment are plotted in Figure 8. The surface a rea was A min = 0.0648m 2 ~ A~ l.08m 2 = A max, corresponding to 0 ~ c ~ 150 carriers added. For T min < T < T max = 25 we plot the Pareto Front along with the substrate removed by the uncontrolled (Q = canst) reactor in subplot a), the relative efficiency of the off-on controlled reactor compared to the uncontrolled reactor in subplot b), and the optimal switching t ime T 0 *witch in subplot c). Qualitatively, the results for added carrier chips are the same as for the previously studied case. As intended, adding biofilm growth area accelerates the treatment process, i.e. quantitative differences are observed. Increasing the biofilm growth area A makes treatment more efficient: the Pareto Fronts computed for larger values of c lie under the fronts for smaller values of c. For small treatment times T , adding carriers increases the efficiency of the controlled system compared to the uncontrolled one with Q = const. The improvement, however, remains modest, below 5%. For larger treatment times T , the improvement achieved by a carrier augmented controlled system compared to the uncontrolled system becomes smaller. The optimal switching time decreases with an increased number of carrier chips added.
Discussion. Using off-on functions to control the process, the largest improvement compared to an uncontrolled reactor was found for the mid-range treatment time of 10-15 days. For short treatment times the throughput in water is high and a substantial amount of substrate might be washed out. On the other hand, for long treatment times the bacteria are exposed to the substrate for a long period of time and can deplete to a great extent. A small improvement is obtained with the off-on functions by allowing the bacteria to establish themselves by consuming the substrate that is already in the reactor before the flow rate is turned on.
The dependency of the reactor performance on the initial suspended biomass and initial biofilm thickness was investigated. A significant difference could be observed between the smallest and largest amount of biomass initially present. Nevertheless, the relative improvement in reactor performance compared with the efficacy of the constant flow rate was smaller than 12%, with most cases ranging between 1% and 5%. Overall, the efficiency of controlling the process as implemented here is not great. Moreover, the initial state that has some effect on optimal reactor performance is very difficult to assess in applications, including controlled laboratory experiments. Therefore, for practical purposes it is doubtful that reactor performance can be reliably improved by controlling the flow rate, compared to the much simpler setup of emptying the buffer vessel at a constant flow rate.

4.4.
Sensitivity of reactor performance with respect to perturbations of the constant flow rate. The previous results lead to the question how slight perturbations of the constant flow rate affect the reactor performance. To investigate this, we study for fixed T the effect of random perturbations of the constant flow rate Q on the amount of substrate in the effluent, i.e. objective z 1 . For given T we calculate the constant flow rate Q = V b /T and perform random perturbations of this value Q with up to 1%, 10% and 20%. To this end we divide the interval [0, T ] into 10 subintervals in which the flow rate is constant. We always maintain T 0 Q(t)dt = V b (0). See Figure 9 for an example. Table 2. Difference between maximum/minimum value of z 1 for random perturbations of up to 1%, 10% and 20%, and the z 1 obtained for a constant flow rate Q.

1%
10% 20% | max(z 1 ) − z 1 (Q const )|/z 1 (Q const ) (%) 0.0513 0.6837 1.7189 | min(z 1 ) − z 1 (Q const )|/z 1 (Q const ) (%) 0.0620 0.4688 0.7439 For each percentage class we randomly choose N = 1000 such perturbations of the constant flow rate. In Figure 10 we plot the histograms of the reactor performance  index. Table 2 shows the la rgest distances from z1 ( Qconst) for each of the percentage classes. La rger perturbations cause a wider range in z1, i.e. the amount of untreated substrate in the effluent. However, this range is not symmetrically centered around the value obtained with constant Q. The higher the perturbation, the more the respective histograms appear skewed with respect to the constant flow rate. A larger perturbation is more likely to worsen the value of z1 than to improve it . This effect can also be seen for the larger perturbations in Table 2, where the differences between the maximum values and constant Q are larger than the differences between the minimum values and constant Q. Nevertheless, the largest changes in z 1 are 1-2 orders of magnitude smaller than the actual perturbations in Q, i.e. reactor performance is robust with respect to small deviations from the constant flow rate. In all cases, the optimal reactor performance obtained with off-on function control is lower than the results obtained in this experiment.
5.1. Q switching between 0 and Q max in discrete time intervals. Inspired by the chattering control for measurable, bounded Q, we now define Ω as the set of non-negative piecewise constant functions, that can switch between 0 and Q max a finite number of times. For a fixed final time T the interval [0, T ] is divided into n equidistant subintervals In each subinterval, the function Q can either take the value 0 or Q max . Moreover, it has to satisfy the reactor condition T 0 Q(t)dt = V b (0). Thus, the number m of intervals in which Q = Q max is then m := nV b (0)/(T Q max ), where the parameters need to be such that m < n is an integer. We refer to this class of functions as zero-max functions Ω Z .
Numerical method. An initial guess for Q in the class of zero-max functions is chosen randomly. We employ a greedy algorithm [8] to find the optimal solution with respect to the objective function T 0 QSdt. The greedy algorithm explores the neighborhood of the current best guess, and makes a locally optimal choice to update the best guess. This process is iterated until no further improvements are possible. For this purpose, we define a neighborhood of a zero-max function as the set of all zero-max functions which differ from the current function in exactly one interval with Q = 0 and one with Q = Q max .
The Greedy Algorithm requires numerous solutions of the reactor model and is, therefore, computationally extremely expensive. The lion's share of computing time is here incurred by the computation of the substrate flux into the biofilm, requiring the solution of a two-point-boundary value problem. To ease this and make computations feasible, we use the flux approximation (8) instead of (5). The underlying model then, in fact, becomes the Freter model [10,18], albeit with different detachment and re-attachment terms.
Results. In Figure 11a we plot a randomly chosen initial guess for Q along with the optimal solution found by the Greedy Algorithm after 104 local decisions in Figure 11b, for T = 10, Q max = 0.01, n = 50, m = 25. The optimal solution starts with a short off-period with Q = 0, switches to Q = Q max until the reactor is emptied and then switches back to Q = 0. The Greedy algorithm by construction finds a local minimum. Therefore, we repeated the simulation five times with different randomly chosen initial data. In all cases the optimal solution found was the same.
Within the class of zero-max functions tested here, the optimal control function is the one that resembles the optimal off-on function found previously the most. Of course, the choice of switching time here is constrained by the discretization of the interval [0, T ]. For comparison, we computed also the optimal, grid independent switching time from 0 to Q max where it stays for a period of V b (0)/Q max . This value is found to be t = 0.1663 days, within the first subinterval of the time grid. To assess the efficiency of this control we computed for comparison the optimal offon function (33), but with the flux approximation (8) instead of (5). The optimal rea.ctor performa nce for zero-max functions with m = 25 was found to be about 1.5 times higher than the optimal reactor performance found for functions of type (33). This could be improved be refining the time grid, which however, would incur a substantia l increase in computing time. switches from 0 to Qma.x after the first subinterval and stays there until the reactor is empty. T he reactor performance is worse than for the optimal off-on functions.
For smaller T , on the other hand, the optimal zero-max control function is distinctively different, see Figure 12. It is characterized by switching several times between 0 and Qmax for small t « T. The number of such jumps in the control functions increases for smaller T , and is in part due to the smaller subintervals.
Another important aspect is the increase of Qmax for smaller T due to m = 25 being kept the same. In Figure 12 Qmax varies between 0.1and0.025 for T = 1 and T = 4, respectively. We performed additional simulations with a fixed Qmax = 0.05 and a fixed number of max-intervals m = 25 for T = 2 and T = 4. In order to obtain 25 max-intervals the number of subintervals for T = 4 was increased to n = 100. The optimal zero-max control function found by the greedy algorithm is the same for these two treatment times up to T = 2, see Figure 13. Thereafter, for T = 4 the function remains at 0 until the end. T hus, it seems that the value of Qma.x strongly influences the optimal zero-max solution rather than the length of the treatment time. T he optimal zero-max solutions for small T are characterized by many rapid switches between 0 and the maximum feasible flow rate, mimicking the chattering control with infinitely many switching points that characterizes the solutions to the original infinite-dimensional control problem (22). In this finitedimensional setup, the number of such switches possible is bounded by the number of degrees of freedom n and the minimum duration of a single on-or off-interval is bounded by T/n. T he locally optimal solution found with the greedy algorithm can depend on the initial data. However, there is not much variation in rea.ctor performance between such local optima. The Pareto Front obtained with zero-max functions with a O N OPTIMIZATION OF S UBSTRATE REMOVAL IN A BIOREACTOR 1161 fixed number of subintervals , as described here, lies above the Pareto Front of the optimal off-on functions, and above the rea ctor performance achieved with constant flow rate, see Figure 14, implying a worse control of the reactor, even compared to the uncontrolled constant flow rate. Furthermore, the increased number of switches between 0 and Qmax in Figure 12 is unrea listic for practical applications.  5.2. Optimal control softwa r e. Optima l control problems with singular arcs, state constraints and discontinuous controls are known challenges for general purpose optimal control software [5]. Such software packages, in particular those based on the direct method, can be understood loosely as searching for the optimal control within a finite dimensional class of functions, usually described by the number of basis functions used to a pproximate the optimal control, e.g. the number of subintervals into which [O,Tmax] is divided. For comparison purposes, we applied two such general purpose optimal control solvers to the problem at hand in Section 3: ACADO Toolkit (Version l.0.2613beta, June 2011, [14]) and DIDO (Version 7.3.6, [30]) . As in Section 5.1 also here the flux approximation was used. We were not able to obtain a robustly converging solution with the former software, probably because of the difficulties to mimic a chattering control. T he latter, instead of computing the solution, is designed to determine candidates for solutions which are then to be subjected to further tests. In Figure 15 we plot two such candidates that we obtained for T = 10 days, discretized with 50 and 100 nodes, wit h 0 :::; Q :S 0.05 m 3 / day. T here are oscillations in the beginning and in the end of the treatment period with a relat ively smooth arc in between. T he amount of oscillations is increased for the larger number of discretization nodes and the largest value attained by Q is~ 0.016 « 0.05 m 3 / day. Here, as for the off-on and zero-max functions, the flow rate is turned off for a few time segments in the initia l phase of the treatment. As 1162 ALMA MASIC AN D HERMANN J. E BERL in t he case of zero-max functions, the oscillatory patterns proba bly might be interpreted as an attempt to mimic chattering in the finite dimensional class of control functions and the frequency of the oscillations increases with the number of degrees of freedom. Moreover, along the trajectory found we observe for the Hamiltonian H ~ 0, which in accordance with (31) might be taken as an indication that the computed solution is optimal for a neighboring regularized problem. To investigate this, however, is beyond the scope of this study. None of the numerical candidate functions that we were able to obtain for the problem at hand were competitive, with respect to reactor performance, with the much simpler optimal off-on function of Section 4.1. On a personal computer, the t ime to compute a single point of the Pareto Front with off-on functions, using the solution of the two-point boundary value problem for the substrate concentrat ion (5) to compute the substrate flux into the biofilm, is in the order of minutes. Using instead the flux approximation (8) is in the order of seconds. On the other hand, computing the best zero-max function with the substrate flux approximation for a given T is in the order of hours; moreover, the solution depends on the initial guesses, i.e. the problem needs to be solved several t imes. Similarly, the time it takes to solve for a single weight w the optimal control problem with optimal control software is in the order of days. These numbers will multiply drastically if the full two-point boundary value problem is solved to compute substrate fluxes. Taking these considerations into account, refining the t ime discretization in the methods of Sections 5.1, 5.2 to achieve results comparable to or better than the results achieved with the much simpler off-on function of Section 4.1 becomes prohibitive.   (8) and (5), and the Pareto front for for zero-max functions with (8).
Furthermore, the flux approximation (8) was used in Sections 5.1 and 5.2 to replace (5) because it is computat ionally much cheaper. It has been shown previously that this formula is able to give a good approximation for a wide range of realistic biofilm model parameters and a wide range of biofilm thicknesses and substrate concentrations. The flux approximation of the Monod kinetics resembles the flux of a zero-order kinetics biofilrn model [2]. T his is likely a good a pproximation in the very early stages of the control period, where the substrate concentration is still relatively large and the biofilm thickness small [20]. On the other hand, the goal of the optimizat ion is to reduce the substrate concentration, i.e. to rea.ch low substrate concentrations S(t). At low substrate concentration, e.g. in the later stages of the control interval, Monod kinetics is known to be better a pproximated by first order kinetics [20]. To investigate the suitability of the approximat ion (8) for the problem at hand, we compare in Figure 16 the Pareto fronts for off-on functions computed with (8) and (5), substrate in effluent for constant flow rate Q = Vi, (O)/T computed with (8) and (5), and the Pareto front for zero-max functions with (8) . T hese results indicate, that substrate in effluent computed with the flux approximation overestimates substrate removal. This mandates that the computat ionally much more expensive (5) should be used instea.d of (8) . In the context of the a pproaches of Section 5.1, 5.2, however, this is prohibitive. Therefore, zero-max functions have no advantages compared to the off-on functions, which yield better solutions to the optimal control problem and are computationally more efficient. 6. Conclusion. In this paper, we formulated and studied an optimal control problem for a bioreactor with biofilm and suspended biomass with a limited amount of substrate and treatment time. The bioreactor model is a biofilm extension of the Freter model with wall attachment. For a one-dimensional single species single substrate model, the optimal control problem was formulated with two objectives: to minimize the treatment time and the amount of substrate in the effluent. The flow rate between the buffer reactor and the bioreactor was selected as the control for the problem. All results were compared against the constant flow rate, which characterizes the uncontrolled treatment. For computational reasons, of particular interest were the off-on functions, which allow the biomass in the bioreactor to establish itself while the flow is turned off before the flow is turned on and the buffer reactor is emptied. We reached the following conclusions: • The Pareto Front for off-on functions shows that for small treatment times the optimal removal efficiency of the reactor changes drastically, whereas it levels off for larger T . For large T almost all the substrate can be removed, whereas for small T a considerable amount leaves the reactor untreated. • The optimal switching times for the off-on functions, at which point the flow is turned on, are small relative to the treatment time. Furthermore, the improvement in reactor efficiency, compared with the constant flow rate, is modest. • Other classes of admissible functions, such as measurable and bounded or measurable and unbounded control functions, led to solutions that are inapplicable/infeasible from a practical point of view. • Calculations with an approximated expression for the substrate flux into the biofilm, used in order to reduce computational complexity, overestimated the substrate removal in the reactor. • The potential to improve reactor performance is due to an initial transient period in which the bacterial populations adapt to the environmental conditions in the reactor. Therefore, control potential depends on the initial state of the system which is typically not known, even in controlled laboratory experiments. • In summary, we come to the conclusion, that the performance increase that is gained by controlling the reactor flow rate is not a significant improvement compared to the uncontrolled, constant flow rate. Moreover, it was shown that the reactor performance is very robust with respect to small deviations from the constant flow rate. These results are based on a comprehensive numerical study for which several hundreds of scalar optimization problems were solved, each of which required numerous simulations of the underlying mathematical model. In these simulations, a small laboratory reactor was mimicked and a generic set of reaction parameters was used. It remains to be investigated whether these results carry over to other setups as well, such as full scale treatment reactors.