A Geometrical Approach for the Optimal Control of Sequencing Batch Bio-Reactors

In this work, we consider an optimal control problem of a biological sequencing batch reactor (SBR) for the treatment of pollutants in wastewater. This model includes two biological reactions, one being aerobic while the other is anoxic. The objective is to find an optimal oxygen-injecting strategy to reach, in minimal time and in a minimal time/energy compromise, a target where the pollutants concentrations must fulfill normative constraints. Using a geometrical approach, we solve a more general optimal control problem and thanks to Pontryagin’s Maximum Principle, we explicitly give the complete optimal strategy.


Introduction
The biological treatment of organic and/or chemical pollutants contained in wastewaters is the transformation of the biodegradable material (also called biomass) in sludge. The principle of biological treatment is to put together microorganisms and pollutants in reactors in which it is possible to control environmental conditions. Due to the simplicity of its implementation, the biological pathway for the treatment of wastewaters is widely used. In this process, bacteria agglomerate into flocs in the reactive part of the system while in a second part, in the absence of agitation, they settle under the effect of gravity. The sludge is thus separated from the treated water and recycled or withdrawn, while the clean water is rejected in the environment.
Several technologies using these principles for the biological treatment of wastewaters have been developed. One such technology relies on the biological Sequencing Batch Reactors (SBR). By separating in time rather than in space the above-mentioned treatment and separation phases, the SBRs, compared to their continuous counterparts, allow a better control of the process and therefore a better quality of the treated effluent.
The price to pay is that the considered system is operating in batch rather than in continuous mode. Therefore, it requires the use of an upstream storage basin for wastewaters which arrive continuously to the treatment plant.
As a consequence, a major issue linked to this process is to minimize the time during which the process is unavailable, i.e. the total reaction time required to treat a batch of wastewater, under the assumption that the settling time is constant. The biological treatment of waters requires different aeration conditions in order to remove both organic carbon and nitrogen compounds. Several optimal control algorithms for the optimization of SBRs have 370 A GEOMETRICAL APPROACH FOR THE OPTIMAL CONTROL where the state s = (s 1 , s 2 ) belongs to (0, +∞) 2 and the control variable u belongs to an interval [0, u max ] with u max > 0. The admissible controls are the L ∞ functions u defined on an interval [0, t(u)] with values in [0, u max ], we write for short u(·) ∈ L ∞ ([0, t(u)], [0, u max ]). The functions ρ i and f i , i = 1, 2, are assumed to satisfy the following hypothesis.

Hypothesis H
• ρ 1  Clearly, this hypothesis guarantees existence and unicity of the solutions of (1) for fixed initial conditions and control function, as well as the invariance of the positive quadrant (0, +∞) 2 . Note also that the variables s i are always non increasing.
We consider a square target C 0 defined by: where a 1 and a 2 are positive numbers, and a continuous nonnegative function J defined on [0, u max ]. Given an initial condition s 0 ∈ (0, +∞) 2 , we introduce the following optimal control problem in free time, When J ≡ 0, this is a minimum time problem. However even when J ̸ ≡ 0, it is possible to consider (P 0 ) as a minimum time problem up to a modification of the control system. And actually we can simplify drastically the problem thanks to change of variables and of time that we describe now. Consider the change of variables s → S defined by: In coordinates S the target writes as (see Figure 1) Now, for a given u(·), set: Using τ as time-parameter, the function S(τ ) is solution of the control systeṁ Noticing that τ (t(u)) = C(u), the optimal control problem (P 0 ) can be written as: S(·) solution of (3) s.t. S(0) = 0 and S(τ (u)) ∈ C S .
Note that the problems (P 0 ) and (P) are equivalent in the sense that the minimum value of their objective functions are equal, and that the optimal solutions u(·) (if they exists) are the same for both problems up to timereparameterization. It is then sufficient to solve the simplified problem (P), which will be the object of the next section. Before doing so, let us make two remarks.
In this case the problem is trivial (one-dimensional problem). Moreover, in the applications the target is considered to be small w.r.t. the initial values of the state. Thus in the sequel we always assume S i > 0 for i = 1, 2, that is, S ∈ (0, +∞) 2 . • The problem (P) is not convex in general. As a consequence, the existence of optimal solutions is not guaranteed by the usual theorems.
The second remark suggests to solve first the problem for the convexified system, where conv(A) denotes the convex hull of the set A.

Resolution by convexification
In this section we first solve a class of convex minimum time problems (subsections 3.1 to 3.3), and we use these results in 3.4 to solve Problems (P) and (P 0 ).

The convex problem (P co )
Let V be a compact and convex subset of R 2 + having a non-empty intersection with the positive quadrant, i.e. V ∩ (0, +∞) 2 ̸ = ∅. We introduce the control system defined by: and the associated minimum time problem, Note that for the control system (4), the set of states that are reachable from 0 is the cone R + V . Since we have assumed S ∈ (0, +∞) 2 and V ∩ (0, +∞) 2 ̸ = ∅, this reachable set intersects the target C S . It then results from standard results (for instance [10, Th. 6.2.1]) that there exists optimal solutions to the problem (P co ).
Let us now apply Pontryagin's Maximum Principle (see [11] for instance). We introduce first the Hamiltonian H := H(S, p, p 0 , v) where p = (p 1 , p 2 ) ∈ R 2 and p 0 ∈ R, as Consider now an optimal control v(·) and the associated trajectory S(·). There exist t(v) Moreover, the following maximization condition holds: The two conditions above imply Besides, the following transversality condition is satisfied: which implies that p ∈ (R + ) 2 .
To summarize, with every optimal trajectory S(·) is associated a four-tuple which satisfies all preceding conditions. Such a four-tuple is called an extremal, and its first component S(·) an extremal trajectory. An extremal is said to be normal if p 0 ̸ = 0, abnormal otherwise.

Lemma 1
There are no abnormal extremals for problem (P co ).

Proof
Assume that p 0 = 0. Then p ̸ = 0, and from (7), we have p T v(·) ≡ 0. It follows that the adjoint vector p ∈ (R + ) 2 is orthogonal to v, which implies that v(t) is either horizontal or vertical for all t ∈ [0, t f ]. In both cases, the extremal trajectory starting at the point S(0) = 0 stays on one of the axis {S 1 = 0} or {S 2 = 0}, and so can not reach the target.
Hence p 0 > 0 and, by (7), p ̸ = 0. We can multiply the Hamiltonian H by 1/∥p∥ and assume that p belongs to the quarter of the unit circle S 1+ = S 1 ∩ (R + ) 2 . The condition (6) on p can be written as Geometrically, ∂V p is the subset of points in V that have the largest coordinate p in the base (p, p ⊥ ). In other words, ∂V p is the part of the boundary of V where the outward normal vector is p, see Figure 2. Since V is convex, ∂V p is either a singleton or a line segment. Note also that, when p is one of the vectors e 1 , e 2 of the canonical basis,

Remark 1
We could also define ∂V p in terms of supporting plane. Recall that the support function of V is the function and that, for p ∈ S 1 , the set H(p) = {v ∈ R n : v T p = h V (p)} is called a supporting hyperplane with exterior unit normal vector p. The intersection H(p) ∩ V is always non empty. Thus the maximized Hamiltonian can be written Let us now analyse in details the possible form of the extremal trajectories. From (8), the target point S(t(v)) of an extremal trajectory is always on the edge of the target. Two cases may occur: either S(t(v)) is on the corner of the target or it is interior to the edge. We will begin by the second case and introduce first a definition.

Definition 1
We say that the point S is under (respectively above) the cone R + ∂V p if either S ∈ R + ∂V p or S lies between the cone R + ∂V p and the horizontal axis (OS 1 ) (respectively the vertical axis (OS 2 )).

Lemma 2
Let S(·) be a minimizing trajectory of (P co ) with control v(·) such that S(t(v)) ̸ = S.
In that case, S is above the cone R + ∂V e2 .

Proof
We only give the proof for the case (i), the same argument can be used to prove (ii). Let us consider an extremal trajectory S(·) that reaches the target at a point on the open half-line {S 1 = S 1 } × {S 2 > S 2 } and p an adjoint vector associated to S(·). Then, from the transversality condition and since ∥p∥ = 1, we derive that p = e 1 = (1, 0). This implies that v(t) ∈ ∂V e1 a.e. and that v 1 (t) = v max 1 . As a consequence, the trajectory S(·) satisfies: Integrating equation (10)

Lemma 3
Let S(·) be a minimizing trajectory of (P co ) with control v(·) such that S(t(v)) = S. Then, there exists p ∈ S 1+ such that S ∈ R + ∂V p . Moreover, t(v) = ∥S∥ ∥v S ∥ and we have the following alternative: Two cases must be distinguished.
The uniqueness of p 0 allows to conclude that t(v) does not depend on the control v and that all optimal trajectories reach the point S at the same time. In particular, ). This implies t(v) = ∥S∥ ∥v S ∥ . Two possibility occurs: either ∂V p is a segment containing in its interior the point v S , and then v(t) ∈ ∂V p for a.e. t ∈ [0, t(v)]; or ∂V p is the singleton {v S } and in this last case v(·) ≡ v S . • If ∂V is not smooth in v S , then there exist θ 0 and θ 1 in [0, π/2] such that the normal cone to V at v S is written as: If p ∈ {e iθ : θ ∈ (θ 0 , θ 1 )}, then ∂V p = {v S } and so v(·) ≡ v S . The trajectory associated to such a constant control reaches the target in time t(v) = ∥S∥ ∥v S ∥ . In the case where p = e iθ0 , two situations may occur: either ∂V e iθ 0 is equal to {v S } and in this case v(·) ≡ v S ; or ∂V e iθ 0 is a segment. In the latter case, v S is an extremity of the segment and necessarily v(·) = v S (it is the only control with values in ∂V e iθ 0 that allows to reach S). As a consequence t(v) = ∥S∥ ∥v S ∥ . The case where p = e iθ1 is treated in the same way.

Optimal synthesis of (P co )
We are now in a position to construct the optimal synthesis of the convex problem (P co ). We distinguish three different subsets of the positive quadrant (0, +∞) 2 .

Definition 2
We define the partition Z 0 ∪ Z 1 ∪ Z 2 of (0, +∞) 2 by setting: • Z 0 is the union of all cones R + ∂V p for p different from e 1 and e 2 , that is, • Z 1 is the set of points of (0, +∞) 2 lying under the cone R + ∂V e1 ; • Z 2 is the set of points of (0, +∞) 2 lying above the cone R + ∂V e2 . Figure 3. Partition of (R + ) 2 into three regions Z 2 (in yellow), Z 0 (in green) et Z 1 (in blue).

Theorem 1
Let S(·) be a minimizing trajectory of (P co ), and v(·) be the corresponding control defined on [0, t(v)].
where ∂V p is the largest segment included in ∂V and containing v S in his interior; the trajectory S(·) reaches the target C S at S, i.e. S(t(v)) = S.

Remark 2
If ∂V contains no segment then the optimal control is constant, and its value only depends on the position of the point S with respect to the set V . It can be determined geometrically, see Figure 4.

Proof
If S ∈ Z 0 , then S is neither under R + ∂V e1 nor above R + ∂V e2 , which implies that S(·) satisfy the hypothesis of Lemma 2. As a consequence, Lemma 3 applies and we get the conclusion.
Let us consider the case S ∈ Z 1 (the case S ∈ Z 2 can be treated in the same manner). If S ̸ ∈ R + ∂V e1 , then only the conclusion of Lemma 2-(i) applies and the result follows. Assume now that S ∈ R + ∂V e1 . Two situations 376 A GEOMETRICAL APPROACH FOR THE OPTIMAL CONTROL may occur: either S(t(v)) = S, and then from Lemma 3 we deduce that v(t) ∈ ∂V e1 a.e. and that t(v) = ∥S∥ ∥v S ∥ ; or S(t(v)) ∈ {S 1 = S 1 } × {S 2 > S 2 }, and it results from Lemma 2-(i) that v(t) ∈ ∂V e1 a.e. and that t(v) = S1 v max 1 .
Since ∥S∥ ∥v S ∥ = S1 v max 1 , both cases give the same conclusion., and the proof is completed.

Back to problem (P)
We come back now to problem (P). We choose V to be the convex hull of the set {F (u) : u ∈ [0, u max ]} and we denote by (P co ) the associated minimum time problem. Thus we are considering two minimum time problems with the same initial condition S(0) = 0 and the same target C S but the control system in (P) iṡ

Proposition 1
Let t min (P co ) and t min (P) be the minimal times corresponding to problems (P co ) and (P), respectively. Then, t min (P co ) = t min (P).

Proof
Since the set of admissible velocities for (P co ) contains the one of (P), we have t min (P co ) ≤ t min (P), so we have to prove the converse inequality. From Theorem 1, (P co ) admits at least one optimal solution S * (·) associated with a constant control The solution S(·) ofṠ(t) = F (u(t)), S(0) = 0, is an admissible trajectory for the problem P and satisfies S(t * ) = S * (t * ) ∈ C S . We deduce that t min (P) ≤ t * = t min (P co ), which concludes the proof.
This result allows to solve the original problem (P 0 ). Indeed, let us define the homeomorphism φ by: ) .
Then the optimal value of the problem (P 0 (s 0 )) obtained by starting from an initial condition s 0 is equal to the minimum time of the problem (P co (φ(s 0 ))) obtained for the target C S with S = φ(s 0 ). As for the minimizing controls, note that u(·) is minimizing for (P(φ(s 0 ))) (and then for (P 0 (s 0 ))) if and only if v(t) = F (u(t)) is minimizing for (P co (φ(s 0 ))). We then obtained the corresponding optimal control of (P 0 (s 0 )) through the time-reparameterization (2). In particular, the construction (13) implies the following property.

Corollary 1
For all s 0 ∈ R 2 + , there exists an optimal piecewise constant control for problem (P 0 ) with at most one discontinuity.

Application
With respect to the model proposed by Mazouni,[7], we assume here that nitrification and denitrification are well controlled by oxygen and that the constraint on the organic matter needed for denitrification is ignored. This simplification allows us to handle a differential system in two dimensions and to consider an optimal control problem in minimal time and energy.
Our objective is as follows: given normative constraints s i > 0, from any initial condition x i (0) > 0, s i (0) > s i , one seeks to synthesize the control u(t) such that s i be smaller or equal to s i > 0 in minimizing a function of both time and energy.
The functions µ 1 and µ 2 are the degradation rates of the substrates which concentrations must be controlled. The choice of growth functions essentially depends on the nature of the treated substrates. In the context of wastewater treatment, the limitation of the nitrogen oxidation by oxygen (or nitrification) is generally modeled by a Monod-type function µ max O2 Ks+O2 multiplied by the growth function. Similarly, the substrate inhibition by oxygen is described by the µ max Ks Ks+O2 . Thus, the limitation of the substrate degradation is modeled by the product of the last function in oxygen and a Monod-type function in nitrogen.

Remark 3
If a substrate s is inhibiting the process at high concentrations, the inhibition can be modeled by the Haldane-type growth function µ max s Ks+s+KI /s 2 .

The optimal control problem
We consider cost functions of the form: where α is a given nonnegative parameter. We take J(u) = J α (u) = αu 2 , in problem P 0 of Section 2. This cost function expresses a compromise between the time to reach the target and the oxygen energy consumed during this duration. For α ≥ 0, we define then the following optimal control problem: s(·) solution of (17) s.t. s(0) = s 0 and s(t(u)) ∈ T .
T is the target defined by where a 1 and a 2 represent the threshold allowed to substrate concentrations in the rejected water. With the previous choices of f 1 , f 2 and J, all the curves in the set are strictly convex curves. Moreover, from Remark 2, the optimal control corresponding to problem (P α ) is constant. The region Z α 2 is empty, and Z α 0 and Z α 1 form a partition of R 2 + . Then, the optimal control value depends on the region S belongs to. The initial concentrations s 0 1 = s 1 (0) et s 0 2 = s 2 (0) are both above the threshold which defines the target T , see Section 2. This means that s 0 1 ≥ a 1 and s 0 2 ≥ a 2 . That's why, we will restrict our study only to the subspace D defined by: D = {s 0 ∈ R 2 + , s 0 1 ≥ a 1 and s 0 2 ≥ a 2 }. Therefore, we consider the restriction of the homeomorphism φ on D defined by: ) . We These regions define a partition of D. For all s 0 ∈ Z α 1 , the minimizing trajectory s u (·) starting from s 0 and corresponding to the control: reaches the target T at a point in {s 1 = a 1 } × {s 2 < a 2 }. If s 0 ̸ ∈ Z α 1 , the minimizing trajectories starting from s 0 reach the target at the corner (a 1 , a 2 ).

Numerical results
In the following, the function f 1 is of Monod-type, the function f 2 expresses the effect of oxygen inhibition at high concentrations in oxygen. They are defined on [0,1] by: The target T is given by T = [0, 0.1] × [0, 0.1], the initial condition is (s 0 1 , s 0 2 ) = (0.65, 0.5). In the following, we consider the cases α = 0 and α = 20, first when µ 1 et µ 2 are of Monod-type and then when µ 1 is of Haldane-type and µ 2 is of Monod-type.

4.3.1.
Case where µ 1 and µ 2 are of Monod-type In this part, we assume that the growth functions µ 1 and µ 2 are defined on [0,1] and are of Monod-type: We present in Figure 5 and Figure 6 the regions Z α 0 and Z α 1 with the corresponding regions Z α 0 and Z α 1 , for α = 0 and α = 20, respectively.
In the case α = 0, the initial condition s 0 is in Z 0 0 , the minimizing trajectory reaches the target at the point (a 1 , a 2 ). But, when α = 20, s 0 ∈ Z 20 1 the minimizing trajectory reaches the target on the edge {s 1 = 0.1}×]0, 0.1[. If α = 20, the minimizing trajectory takes more time to reach the target than in the case of minimal time trajectories (α = 0), as shown in Figure 7. In this case, since we take into account the energy, the control value is lower than in the case α = 0. The growth of x 2 is thus less inhibited than in the case α = 0. The rate at which s 2 decreases is then faster than the rate at which s 1 decreases.
Interestingly, we note the existence of a critical value α * such that t α is constant up to α * and increases for α ≥ α * . Indeed, for α ≤ α * , s 0 is in the region D\Z α 1 . Then, the minimal time to reach the target is given by ∥v∥. This value does not depend on α. Now, if α > α * then s 0 ∈ Z α 1 and the time taken for the minimizing trajectory to reach the target is given by

4.3.2.
Case where µ 1 is of Haldane-type and µ 2 is of Monod-type Figures 8 and 9 are obtained with a growth function µ 1 corresponding to the inhibition of the substrate (function Haldane-type) and a growth function µ 2 of Monod-type. In this case, the region Z 0 1 (respectively Z 0 1 ) is more important than the region Z 0 0 (respectively Z 0 0 ). Moreover, In this case, the shapes of the trajectory discussed in the previous section are amplified in the sense the decrease of s 2 is still faster. Indeed, s 1 inhibits the growth of x 1 which diminish its decrease rate. Putting the inhibition on x 2 would lead to the opposite tendency. From the different cases we considered, we noticed that the regions Z α 1 increase when α increases. When working with two species of different nature, a first one with an inhibitory growth function in substrate and a limited one in oxygen, and the second species with limited growth functions in substrate and oxygen, the reduction of oxygen consumption leads to an increase in the duration of the reaction, for a large set of initial conditions. In fact, for s 0 ∈ Z α 1 , the necessary time to reach the target depends on α and is no longer minimal. More general comments are difficult to establish since the optimal trajectories depend not only on the objective function (and thus on α) but also on model structure parameters. Finally, the most important result is that in any case, the optimal control is constant.

Conclusion
In this paper, we propose a model that involves two biological reactions, one being aerobic and the other anoxic, required for processing two types of substrates. The control is the dissolved oxygen considered as an input taking continuous values. The objective is to minimize a compromise between the total reaction time and the energy consumption. The dynamics are described by a differential system coupled in the control variable and decoupled in the state variables. Our approach consists in studying first a more general optimal control problem with a cost function being the integral of a positive function J, using a geometric technique. Due to the maximum principle, we have characterized the extremal trajectories and we determined the minimal time to reach the target and the corresponding control, according to the position of the target point on the edge of the target. This technique allows us to solve the problem in a very general framework (only the assumption of continuity of the growth functions in oxygen is required). Using these results, we determined the optimal strategy to control the nitrificationdenitrification reactions by oxygen, under the assumption that the organic matter is not needed for denitrification. The mathematical study is followed by numerical simulations which illustrate the obtained results. The results of this first study will help us tackle the problem in three dimensions where the variables are coupled by an additional substrate (carbon).