Optimal control of a diffusion/reaction/switching system

We consider an optimal control problem involving the use of 
bacteria for pollution removal where the model assumes the bacteria 
switch instantaneously between active and dormant states, determined 
by threshold sensitivity to the local concentration $v$ of a diffusing 
critical nutrient; compare [7], [3], [6] in which 
nutrient transport is convective. It is shown that the direct problem 
has a solution for each boundary control $ψ = ∂v/∂n$ and that 
optimal controls exist, minimizing a combination of residual pollutant 
and aggregated cost of the nutrient.


Introduction
In [7] (and [3], [6]) we considered a model of bioremediation using bacteria which are sensitive to the presence/absence of a critical nutrient in switching discontinuously between dormant and active modes.[We will consider only a single "critical nutrient" and assume any other demands are adequately supplied and do not affect the dynamics.]While we observe that the models considered here are oversimplified caricatures of realistic situations, we emphasize that our focus here is entirely on the mathematical considerations arising in this analysis.Noting that the modeling here is perhaps least realistic in the assumed instantaneity of the transitions, we comment that there is considerable variety of bacteria, so some may plausibly resemble the models here reasonably closely, and (compare, e.g., [8]) note that this instantaneity is just an idealization of the 'fast scale' of a multiscale process.
In any case, our focus here is on the mathematical interest of the problem in which switching occurs in somewhat different contexts than is typical of hybrid control.The notion of threshold induced modal switching is characteristic of hybrid systems where the modal switching is typically a control action and the thresholds relate to feedback.In the absence of a Zeno phenomenon (infinitely many transitions in a finite time interval), the system evolution can be followed separately on each of finitely many interswitching intervals.What is somewhat new here, where the switching is instead a feature of the physical modeling of the controlled process. is that each individual bacterium is subject to such modal transitions so in a continuum model there will be, globally, a continuum of modes and some form of the Zeno phenomenon becomes almost inescapable.
We focus on a distributed parameter control problem in which, assuming the initial presence of suitable bacteria in a dormant state, one provides nutrient so as to activate the biomass, thus leading to metabolism of an undesirable pollutant -or cometabolism (cf., e.g., [10], [1]), i.e., converting it to a more innocuous form.In contrast to the earlier model of [7], in which the nutrient transport was purely convective, we now consider diffusive transport with the control ψ entering as a boundary control of the diffusing nutrient concentration.Note that our concern here is with a model in which, as in [3], [7], the bacteria and the pollutant are each stationary in position.Indeed, this is essential for the treatment here.[One might alternatively wish to consider a diffusing biomass, but note that the modal transitions we describe are following the history of individual bacteria, so the mode at t would then not be expressible simply as a function of position and our present analysis would not apply.We do not attempt here to consider that situation.] The dynamics of our model are described by a coupling of four unknown functions: u = u(t, x) (for t ≥ 0, x ∈ Ω) denoting the bacterial concentration, v denoting the concentration of the critical nutrient and p denoting the pollutant, each taking values in R + = [0, ∞), together with another 'Boolean-valued' unknown function w = w(t, x) identifying the mode of the bacteria at (t, x) by setting w = 0 for dormancy and w = 1 for activity.
We will formulate the problem more precisely in Section 3, but briefly describe the dynamics here.We begin with a fairly standard diffusion/reaction equation (3.1) for the nutrient concentration v: We are assuming normal behavior (u t = γ(wu) for some growth rate γ) when the bacteria are active (w = 1), and "nothing" when dormant -although it is interesting to allow for the possibility that only some fraction r < 1 of the bacteria survive a dormant interval to be reactivated at a switching time when w switches: 0 1.The pollutant is metabolized (or cometabolized) by the (active) bacteria at a fixed rate δ so, at each x ∈ Ω, we have Our objective, of course, is to remove as much as possible of the pollutant by having enough active bacteria, accomplishing this by having a high enough nutrient concentration.
We then may consider, as a problem of optimal boundary control, introducing a controlled nutrient flux to activate the otherwise dormant bacteria to remove the undesirable pollutant by metabolism or cometabolism.Adding such a flux of new nutrient makes the boundary conditions equationvbc for the nutrient dynamics (3.1) take the form The function ψ ≥ 0 on Σ is the control.Assuming the nutrient is expensive, the natural cost functional for the control problem has the form where these L 1 -norms are on Ω, Σ, respectively, giving the amount of pollutant remaining at the terminal time and the amount of nutrient added (so we are dealing with minimization over a nonreflexive space).Scalar threshold induced modal switching by a discontinuous elementary hysteron has been well studied, but there are several versions of this; note, in particular, [2] and [11].The version used here follows [5], [8], etc., and we include in Section 2 a review of some properties we will need, especially since we are concerned here with a continuum of such hysterons corresponding to the continuum of spatial points in Ω -and this affects our treatment of the coupling between the dynamics of v and the dynamics of w.

Threshold induced switching: the hysteron
We begin with the scalar situation, here corresponding to looking at a single point x ∈ Ω, and considering the hysteron as an operator W : [y(•), ω] → w(•) (from sensor input and initial mode to modal output) which is to be determined by a set of switching rules with respect to the upper and lower thresholds 0 Thus, if the current sensor input y(t) is above the upper threshold η + or below the lower threshold η − the mode is uniquely specified, independent of past history -but it is history-dependent (the result of the most recent prior switch) when the current input lies between those thresholds.Note that (2.1-c) ensures that the output w can only switch at threshold values and (2.1-a,b) ensure that one must have modal switching: 0 1 when y crosses η + upward and 1 0 when y crosses η − downward.

Proof:
For any interswitching interval (t 1 , t 2 ) we necessarily have y(t 1 ) = η − and y(t 2 ) = η + (or vice versa) so for y(•) continuous, hence uniformly continuous on [0, T ], there must be a positive lower bound on the difference t 2 − t 1 .Hence there must be a upper bound on the number of such intervals and so on the number Var [w] of switchings.
We note that (2.1) need not uniquely determine the output w(•) for a given input y(•), even topologizing w in L 1 [0, T ] (so, e.g., its specification just at the (isolated) switching times could be ignored).The most serious nonuniqueness occurs at anomalous points when y(•) is nontransversal at a threshold.For example, suppose one might have w(t−) = 0 and y(t−) = η + (so y(t) ≤ η + on some (small) interval [τ − ε, τ ]) and suppose also that y(t) ≤ η + on some (small) interval [τ, τ + ε].Then (2.1) permits w = 0 continuing on [τ, τ + ε], but also permits a switch of w : 0 1 at t, then continuing in mode 1 until y might decrease to η − .The treatments in [2], [11], etc., show how to make selections with various nice properties of the hysteron (as isotonicity, rate independence, etc.), but we will here take W to be inclusive: given y(•) ∈ C 0 ([0, T ] and ω ∈ {0, 1} we accept all modal functions and set This defines the elementary hysteron W as an operator which will map The rationale for this definition is that the model is to be viewed as a reduced description of a multiscale process so these discontinuities represent the effects of otherwise unmodeled behavior on a faster scale than we consider.To the extent that there may be uncertainties involved in this reduced modeling, our response is to accept that any of these continuations is as consistent with our available information as any other.Further, this definition will ensure a form of well-posedness, that "the limit of solutions is a solution." In this context we point out that one way to produce such switching is to introduce and consider solutions of the singularly perturbed ODE Then W [y, ω] consists of the limits of such solutions as ε → 0+ and ŷ converges uniformly to y. It is easily seen from examples that, in general, W is neither singlevalued nor continuous.However, the key property we will need, motivating our choice of this set-valued version of the elementary hysteron, is that W is a closed operator when defined by (2.2).

Proof:
Noting the proof of Lemma 2.1, we observe that for large n the number of switching times is bounded and each w n is necessarily piecewise constant with the same values during the interswitching intervals.One then easily sees that pointwise convergence on any dense set implies that the limit w is also piecewise constant and that one has convergence of the switching times: t n,k → tk .If ȳ(t) > η + for some t, then ȳ(τ ) > η + + ε for some ε > 0 and all τ in some neighborhood of t; then uniform convergence gives y n (τ ) > η + so w n (τ ) = 1 for large n so w(t) = 1.Thus, w satisfies (2.1-a.)and, similarly, (2.1-b.).In much the same way one shows switching is only possible for w when ȳ = η ± so one also has (2.1-c.), thus ensuring that w ∈ W [ȳ, ω].
Finally, we have ] is compact when topologized, e.g., by pointwise convergence.

Proof:
Given any sequence (w n ) in W [Y, ω] we first extract a subsequence (wlog, using the same indexing) such that y n → ȳ with each w n ∈ W [y n , ω].Much as in the proof above of Lemma 2.2, each w n is piecewise constant with a bounded number of switching times t n,k and characterized by {t n,k }.By the compactness of [0, T ] one can then extract a further subsequence such that t n,k → tk for each k.Since each w n is piecewise constant with values 0 or 1, this gives pointwise convergence to some correspondingly piecewise constant w and so, by Lemma 2.2, with w We will be considering a somewhat more general hysteron W, acting on continuous functions on Q = [0, T ]×Ω to get {0, 1}-valued functions on Q by applying W independently at each x ∈ Ω.Our application assumes the initial presence of a dormant bacterial population so we have the simplification that the initial nutrient concentration is below the lower threshold whence ω = w(0, x) = 0 for each x.Thus, we will usually suppress the argument ω in considering W or W. Our definition of W is: (2.6) Note that Lemmas 2.1, etc., apply at each x ∈ Ω, but we have not yet explicitly considered properties of the operator W.
3 Formulation of the model In this section we consider the direct problem.For our model we let u = u(t, x) (for t ≥ 0, x ∈ Ω) denote the bacterial concentration at each fixed point x ∈ Ω and let v denote the diffusing concentration of the critical nutrient; these are unknown nonnegative scalar functions.Also unknown is w = w(t, x) identifying the mode of the bacteria at (t, x) by setting w = 0 for dormancy and w = 1 for activity and the pollutant concentration p = p(t, x) For v we have the fairly standard diffusion/reaction equation: where the term −αv (α = constant > 0) represents a degradation with time of the nutrient and the term −β wu represents (metabolism or) cometabolism of the nutrient by active bacteria -of course, with no such (co)metabolism for dormant bacteria.We will consider boundary control, taking as boundary conditions for (3.1): with n the unit outward normal.Thus, the control is the function ψ ≥ 0 on Σ = [0, T ] × ∂Ω, which gives the flux rate at which new nutrient is being supplied to the region.When the bacteria are active, pollutant is (co)metabolized at a fixed rate δ with no such remediative action when dormant.We then have as a family of uncoupled ODEs parametrized by x ∈ Ω; while nominally taking δ to be a positive constant, we will, somewhat artificially, set δ = 0 when p ≤ 0 since even active bacteria cannot affect a nonexistent pollutant.We are assuming normal behavior (u t = γ(u) for some growth rate γ when the bacteria are active (w = 1), and "nothing" when dormant.Thus, independently at each x ∈ Ω, the bacterial concentration u satisfies the ODE with γ(•) representing population growth when the bacteria are active.While one could let γ depend also on the nutrient level v; noting that v > η − where relevant here (i.e., where wu = 0), reasonable modeling of such dependence would not materially affect the analysis and we take the growth term to be a fairly standard nonlinearity: we need assume only that γ(•) is nonnegative and Lipschitzian so 0 ≤ γ(r) ≤ a + br.

2.]
Finally, throughout this we ask that The system (S) is then the coupling of (3.1), (3.3), (3.4) and (3.5) with the BC (3.2) and the initial conditions.Since we are here modeling pollutant and bacteria as spatially fixed, (3.3) and (3.4) are families of ODEs parametrized by x ∈ Ω with (3.5) also applied independently for each x.For our regularity hypotheses on the system (S), see the statement of Theorem 4.2.The relevant solution is discussed in the proof of that theorem, but we already note here that u, v, p will be nonnegative with u, p ∈ L ∞ (Q) and with v uniformly continuous on each The two usual ways to show existence for such a nonlinear problem would be either to seek a fixpoint of some map or to construct approximations and apply compactness.Here the discontinuity of the component W discourages use of the fixed point approach: for hybrid systems with assurance of a finite set of switching times, one usually proceeds by solving separately on each interswitching interval without concern for these discontinuities, but in our present model one expects a continuum of such times collectively although one has a finite set of switching times for each spatial point x.We will proceed by approximation.
Our strategy will be to introduce a finite partition of Ω into sets {Ω j } with characteristic functions χ j and choose x j ∈ Ω j so each x j is in the interior of Ω.We then set y(•) = y P (•) = [y 1 (•), . . ., y J (•)] = Pv with y j (•) = v(•, x j ) and replace w in the system by w P = W[Pv], i.e., w P (•, x) ∈ W [y j (•)] for x ∈ Ω j to get the approximating system (S P ) with the same data.

Existence of solutions
In this section we prove existence of solutions of the coupled system (S) for given data ψ.In general one cannot expect uniqueness of solutions, but the nature of the existence argument does show that limits of solutions are again solutions.It is almost immediate then to obtain existence of optimal controls.
We begin by presenting an existence proof for the simpler approximating system (S P ), omitting some details of the proof which will be made part of Theorem 4.2: Theorem 4.1.Each problem (S P ), as described above, has a solution.

Proof:
The key point here is that we can proceed causally on interswitching intervals.With each w j = w(•, x j ) constant on the interval, there is no difficulty solving the coupled system of DEs (the diffusion equation for v and ODEs for u, p) until (2.1-a,b) would require switching some component w j -or (2.1-c) would permit switching.We will see in proving Theorem 4.2 that each of the functions y j = v(•, x j ) is necessarily uniformly continuous on [0, T ] (no matter how we might make the switching choices if anomalous points would occur) and that there is a fixed bound, depending only on the points x j and the data estimates, on the collective number of switching times.We can therefore proceed to solve (S P ) successively on these interswitching intervals -restarting with the new (temporarily constant) value of w at each switching time -up to t = T .Note that this construction enables us to handle the possibility of reactivation failure: at each switching time where w j switches 0 1 we multiply u by ρ at each x ∈ Ω j .Since we followed the switching rules where relevant, we see that the modal indices resulting must satisfy (3.5) P , i.e., as desired, in addition to (3.1) with (3.2), (3.3) and (3.4) including reactivation failure.Thus we have a solution of (S P ).
With this in hand, we turn now to showing existence of solutions for the model system (S).Theorem 4.2.Assume Ω is bounded with ∂Ω smooth and that the initial data for u, v, p is smooth and nonnegative with v 0 ≤ η − so w t=0 ≡ 0. Let the input flux ψ be a nonnegative finite Borel measure on Σ.Then the solution set of the coupled system (S) is nonempty.This solution set is upper semicontinuous in its dependence on the data.

Proof:
We consider a sequence of problems (S ν ) = (S Pν ) using projections P = P ν with the diameters of Ω P j shrinking and the sets {x j } P becoming dense in Ω.A compactness argument will then permit extraction of convergent subsequences (so v = v ν = v P → v * , etc.) and we show the limits satisfy (S), in particular, that the limit w * is in We begin by considering the solution space of (S), i.e., the a priori regularity of solutions, Starting with v, we write v = v ψ + v u where v ψ is the (fixed) solution of (3.1), (3.2) with the term −βwu omitted; v u is the solution of (3.1) with 0 initial data and no-flux boundary conditionss (ψ = 0).Even with measure-valued ψ, the regularity for the diffusion equation localizes: function v ψ will be in C ∞ (Q) although it need not be bounded.Recalling this is bounded away from any singularity admitted by our choice of control space when ε > 0 so v ψ will certainly be uniformly continuous on any such Q ε .[For future reference, we note that a bound on ψ in [C(Σ)] * gives equi-uniform continuity on each Q ε for v ψ .]We also note that the assumed nonnegativity of initial data and of the input flux ψ ensures that v ψ ≥ 0. On the other hand, we will bound u in L ∞ (Q) and so obtain equi-uniform continuity on Q for v u with a pointwise bound whence compactness in C(Q ε ).We emphasize that this ensures a bound (depending on x ∈ Ω ε but not on the partition P ν ) on the number of switching times for the resulting w(•, x).We further can use the weak maximum principle to observe that the solutions v remain nonnegative: e.g., we take ϕ = min{0, v} as test function in the weak form of (3.1) (noting that ϕv t = ( 1 2 ϕ 2 ) t and ∇f • ∇v = |∇ϕ| 2 ae and that ϕwu = 0 [since w = 0 when 0 = ϕ = v ≤ 0 < η − ] and ϕv = 0 at t = 0) to see that ϕ = 0 on Q.
From (3.4) with γ ≥ 0 we get u > 0 if one has positive initial data, even allowing for reactivation failures which involved multiplication by ρ n with n uniformly bounded for x ∈ Ω.Using the assumed linear growth rate, we also get (pointwise) the upper bound for u noted above when u is bounded at t = 0 by comparison with u t = a + bu.From (3.3) we see that p remains bounded above by its initial value, assumed nonnegative and integrable; the diktat that δ = 0 when p ≤ 0 ensures that p also remains nonnegative.
We now let (P ν ) converge in the sense noted above -increasing the set of sensor points {x j } while further subdividing Ω, letting {x j } become dense and letting the diameters of the Ω j shrink to 0 in the limit.Since v ψ is fixed and we have suitable compactness of the set of {v u : P ν }, we can always extract a subsequence, again denoted by (P ν ), such that v Pν converges uniformly on each Q ε .
For any x ∈ Ω we have x ∈ Ω ε for some ε > 0 and may consider a sequence of sensor points x n(k) ∈ Ω ε converging to x, whence y n(k) = v P (•, x n(k) ) → v(•, x).By Lemma 2.4 we may further extract a subsequence such that w P (•, x n(k) ) converges pointwise to w * (•, x) ensuring (3.5).Given this, the ODE (3.4) gives convergence u P → u and (wu) P → wu satisfying (3.4) at x. Similarly, p P → p * as P converges with p * satisfying (3.3) at x. Since x ∈ Ω is arbitrary, this limit is a solution of (S) as desired.We note, in particular, that p * (t, •) is well-defined on Ω for each t.
Since W is, in general, set-valued, we cannot expect uniqueness.However, the argument above by compactness does show that a limit of solutions of (S) is again a solution once we note, as an implication of the localized regularity, that a bound on ψ in, e.g., [C(Σ)] * is sufficient to ensure that v ψ is in a compact subset of C(Q ε ) for each ε > 0 so the argument goes as before to show that one obtains a solution of (S) in the limit.

Existence of optimal controls
Finally, we return to the bioremediation model (S) = (S) ψ as a control problem with boundary control ψ.As noted, the cost functional J of (1.1) is a natural criterion here, balancing the running cost of supplying an expensive nutrient during [0, T ] against the economic cost of having pollutant remaining at the end of this time.We will actually enlarge the control space to the second dual, admitting measures as controls, so taking ψ ∈ [C(Σ)] * with no change in (1.1) From the arguments for Theorem 4.2 we almost immediately get existence of optimal solutions of this control problem: Theorem 5.1.There exists an optimal control ψ * for which the cost functional J of (1.1) (with fixed c 1 , c 2 > 0) is minimized, subject to ψ ≥ 0 and admitting measures as controls.

Proof:
As usual, we consider a minimizing sequence ψ n so J [ψ n ] → J * = inf J .One has coercivity there: the minimizing sequence is necessarily bounded in what is now a dual space so, applying Alaoglu's Theorem, we may then extract a subsequence which is weak- * convergent to some ψ * ∈ [C(Σ)] * .As in the last paragraph of the proof of Theorem 4.2, we have compactness of (restrictions of) v ψn in each C(Q ε ) so, as there, we may conclude that the arguments of the proof of Theorem 4.2 continue to apply: perhaps for a subsequence, we have convergence to a solution of (S) ψ * .Since the functions p are nonnegative and uniformly bounded on Ω, we have convergence of Ω p to Ω p * by the Dominated Convergence Theorem and note that Σ ψ * ≤ lim inf Σ ψ n by the weeak- * convergence.The standard argument then applies: J (ψ * ) ≤ lim inf J (ψ n ) = inf J and, since one cannot have J (ψ * ) < inf J , we attain the minimum cost with ψ * .
Remark: Essentially the same argument could be used to treat J -optimal control with a variety of other convex constraints -e.g., that the support of ψ be in some given subset of ∂Ω or that ψ be of the form K  1 ψk (t)e k (x)

Lemma 2 . 1 .
If y(•) is continuous on [0, T ], then (2.1) ensures a bound on the number of switching times, hence a bound on the variation Var [w].