*Article* **Optimal Control of a PDE Model of an Invasive Species in a River**

#### **Rebecca Pettit and Suzanne Lenhart \***

Department of Mathematics, University of Tennessee, Knoxville, TN 37996-1320, USA; beccapettit08@gmail.com **\*** Correspondence: slenhart@utk.edu

Received: 6 September 2019; Accepted: 10 October 2019; Published: 15 October 2019

**Abstract:** Managing invasive species in rivers can be assisted by appropriate adjustment of flow rates. Using a partial differential equation (PDE) model representing an invasive population in a river, we investigate controlling the water discharge rate as a management strategy. Our goal is to see how controlling the water discharge rate will affect the invasive population, and more specifically how water discharges may force the invasive population downstream. We complete the analysis of a flow control problem, which seeks to minimize the invasive population upstream while minimizing the cost of this management. Using an optimality system, consisting of our population PDE, an adjoint PDE, and corresponding optimal control characterization, we illustrate some numerical simulations in which parameters are varied to determine how far upstream the invasive population reaches. We also change the river's cross-sectional area to investigate its impact on the optimal control.

**Keywords:** optimal control; partial differential equation; invasive species in a river

#### **1. Introduction**

The need to develop methods for managing invasive species in rivers is growing [1–3], and modelling can give insight into management strategies. Jacobsen et al. [4], Jin and Lewis [5,6], and Gagnon et al. [7] used integrodifference models to represent such populations in rivers. With a more realistic model, Jacobsen et al. concluded that longer stream length and lower flow rates increase the chance of a species persisting while higher flow streams must be even longer to maintain persistence [4]. They also showed that flow velocity variability can produce persistence similar to the mean velocity [4]. Jin and Lewis [5,6] also conducted a theoretical analysis using integrodifference equations assuming periodicity with respect to time. In 2011, they first explored how the critical domain size would affect the species persistence and determined that a periodic dispersal kernel and a weighted time-averaged dispersal kernel have the same effect. However, this only holds for the estimation of the critical domain size [5]. In 2012, Jin and Lewis [6] studied spreading speeds for the same integrodifference equations in [5]. Jin and Lewis again found that the periodic dispersal kernel and the weighted time-averaged dispersal kernel both yield the same spreading speeds. Once again, this only applied for the spreading speeds, and there is no proof that the dynamics of the systems are the same [6]. In [7], Gagnon et al. took an applied approach having a specific species, *Codium fragile*. Gagnon et al. focused on wind-driven spread of floating fragments of this algae. They found that prevalent winds were vital in determining the spreading speeds as well as the fragmentation into non-buoyant and buoyant pieces was the most important feature in determining spread [7].

In a recent paper, Jin et al. built a reaction–diffusion partial differential equation (PDE) model to better represent a population in a river by using not just the population growth, but how the river itself changes over time [8]. The idea was that this new model would better represent a population trying to move up a river. Clifford et al. [9] discussed how water flow affects the ecology of the river and the shift from focusing purely on river flow dynamics to how these flow dynamics create various habitats [9]. Huang et al. [10] used a PDE system to examine how a species that usually lives in the benthic layer of a river may move in the stream using drift. They used three values to measure the persistence of a species in the river and found that the net reproductive rate determined if a population could persist in the river or not [10]. In a later paper in 2017, Jin et al. [11] explored how the cross-sectional area of a river would affect whether the population persisted. Assuming the meandering rivers are periodic, the cross-sectional area changes over time [11]. They found that increasing growth rates increases the chance of persistence by the population and allows for spreading both up and downstream while increasing flow spread decreases the chance of persistence and upstream spread but increases downstream spread.

Jin et al. [8] explored simpler forms of the models by setting the cross-sectional area and water discharge to constant values to see how a population would move in the stream. They began with traveling waves in constant environments and temporally or spatially heterogeneous environments. To see the effect on the population's movement upstream, Jin et al. then varied the water discharge and cross-sectional area between two time periods to illustrate a seasonal period of summer versus winter. The water discharge and cross-sectional area in the summer are two values, and in the winter are other values. To investigate time varying flows and corresponding management, we use optimal control theory on the water flow discharge to control the movement of an invasive species upstream.

We are interested in modeling an invasive species (like carp) which moves upstream. To manage an invasive species in a downstream position in a river, one needs an action to keep that species from moving upstream. Current methods to prevent further expansion include electric fences or nets in the river [1]. In addition, adjusting the flow rate via dams or water release mechanisms may be a reasonable containment mechanism. Poff and Zimmerman [12] mention that extreme changes in river flow can affect a populations survival success in a habitat. They discuss flows changing due to climate and controlling runoff by catchment [12]. Using the Jin et al. [8] model in our work, we choose the water discharge "flow" function to affect the invasive species movement upstream. Taking the advection coefficient as a time varying control function, we will be utilizing optimal control theory to analyze this model and this management mechanism to contain the species downstream.

In the following sections, we first formulate the problem and describe the terms of the PDE model. We define the weak solution and the solution space. We then describe the control set and the objective functional to be minimized. Using a priori bounds on the state in the weak solution space, we prove the existence of the optimal control. Then, we prove differentiability of the control-to-state solution map while giving the sensitivity system leading to the adjoint system. We differentiate the objective functional and derive an optimal control characterization using the sensitivity and adjoint systems. Finally, we present some numerical simulations to show the effects of our optimal control on the movement of the population. Specifically, we illustrate how the population moves upstream with and without control.

#### **2. Optimal Control Problem Formulation**

We use the model adapted from Jin et al. [8]. To formulate our state PDE and corresponding optimal control problem for the weak solutions, our solution spaces are *V* = *L*2(0, *T*; *H*<sup>1</sup> {0}(0, *<sup>L</sup>*)) and *V*<sup>∗</sup> = *L*2(0, *T*;(*H*<sup>1</sup> {0}(0, *<sup>L</sup>*))∗), where the subscript {0} indicates that the functions vanish at *<sup>x</sup>* <sup>=</sup> 0.

**Definition 1.** *A function N* <sup>∈</sup> *<sup>V</sup>* <sup>∩</sup> *<sup>L</sup>*∞((0, *<sup>T</sup>*) <sup>×</sup> (0, *<sup>L</sup>*)) *with Nt* <sup>∈</sup> *<sup>V</sup>*<sup>∗</sup> *is a weak solution of the state PDE:*

$$\begin{aligned} N\_t &= -A\_t(\mathbf{x}, t) \frac{N}{A(\mathbf{x}, t)} + \frac{1}{A(\mathbf{x}, t)} \left( D(\mathbf{x}, t) A(\mathbf{x}, t) N\_\mathbf{x} \right)\_\mathbf{x} - \frac{Q(t)}{A(\mathbf{x}, t)} N\_\mathbf{x} + rN \left( 1 - \frac{N}{K} \right) \\ N(0, t) &= 0 & \text{on } (0, T), \mathbf{x} = 0, \\ N\_x(L, t) &= 0 & \text{on } (0, T), \mathbf{x} = L, \\ N(\mathbf{x}, 0) &= N\_0(\mathbf{x}) & \text{on } (0, L), t = 0 \end{aligned} \tag{1}$$

*in* Λ = (0, *L*) × (0, *T*)*, if*

$$\int\_0^T \langle \mathbf{N}\_t, \boldsymbol{\phi} \rangle dt + \int\_0^T \int\_0^L \left( A\_t \frac{\mathbf{N} \boldsymbol{\phi}}{A} + D \mathbf{A} \mathbf{N}\_x \left( \frac{\boldsymbol{\phi}}{A} \right)\_x + \frac{Q}{A} \mathbf{N}\_x \boldsymbol{\phi} - r \mathbf{N} \boldsymbol{\phi} \left( 1 - \frac{N}{K} \right) \right) dxdt = 0$$

*for all <sup>φ</sup>* ∈ *V, where* , *in the integrand denotes the duality H*<sup>1</sup> {0} (0, *<sup>L</sup>*), *<sup>H</sup>*<sup>1</sup> {0} (0, *<sup>L</sup>*) ∗ *.*

The state *N*(*x*, *t*) represents the population density in the river, *A*(*x*, *t*) is the cross-sectional area of the river at location *x* at time *t*, *D*(*x*, *t*) is the diffusion coefficient, *r* is the population intrinsic growth rate, and *K* is the carrying capacity. The control function *Q*(*t*) represents the water discharge rate. An important note is that the location *x* = 0 is upstream while *x* = *L* is downstream. For our boundary conditions, at upstream where *x* = 0, we assume no population is initially there, as they are invading by moving from downstream to upstream. Likewise, at downstream, *x* = *L*, we assume no flux, meaning the population is not increased from an outside source. To simplify the calculations, the model uses a one-dimensional, spatial domain corresponding to the length of the river with the assumption of homogeneous mixing in any cross-sectional area.

Our goal is to control the water discharge rate *Q*(*t*) to prevent an invasive species from moving upstream. Our control set for 0 < *m* ≤ *M* is

$$\mathcal{U} = \{ \mathcal{Q} \in L^{\infty}(0, T) \mid m \le \mathcal{Q}(t) \le M \}.$$

Because of the damming of rivers, *Q*(*t*) is bounded above by some value *M* and bounded below by *m* positive since water discharge by definition must be positive. We assume we can fully control the water discharge rate by controlling how much water moves through a dam. Note that *Q*(*t*) is not a function of space, but just time. The idea explored in Jin et al. is that the water discharge rate being higher during certain seasons and lower in others affects the population level at various locations [8].

Our objective functional which we intend to minimize is:

$$J(Q) = \int\_0^T \int\_0^L W(\mathbf{x}) N(\mathbf{x}, t) d\mathbf{x} dt + \int\_0^T \epsilon Q^2(t) dt,$$

where  > 0 is small. The coefficient *W*(*x*) is a weight function that is large at *x* = 0 and small at *x* = *L* since our goal is to keep the population downstream (at *x* = *L*). The first integration term seeks to minimize the size of the invasive species with an emphasis on keeping the invader downstream (due to the weight function) by decreasing *N* near *x* = *L*. The *Q*2(*t*) term seeks to minimize the cost of effort needed to stop the species from moving upstream with an  coefficient to keep it small. Thus, the term with *W*(*x*)*N*(*x*, *t*) should dominate the objective functional. We seek *Q*<sup>∗</sup> ∈ *U* such that

$$J(Q^\*) = \inf\_{Q \in \mathcal{U}} J(Q).$$

Assuming *A* is a *C*<sup>1</sup> function and *D* is a continuous function on our domain, we make the following additional assumptions:

$$r > 0, K > 0,$$

$$0 < m \le Q(t) \le M\_{\prime}$$

$$0 < \epsilon\_A \le A(\mathbf{x}, t) \le M\_{A\prime}$$

$$|A\_{\mathbf{x}}(\mathbf{x}, t)| \le M\_{B\prime} \, |A\_{\mathbf{f}}(\mathbf{x}, t)| \le M\_{C\prime},$$

$$0 < \epsilon\_D \le D(\mathbf{x}, t) \le M\_{D\prime}$$

$$0 \le N\_0(\mathbf{x}, t) \le M\_1.$$

To solve the optimal control problem for a given *Q* ∈ *U*, we need the existence of a unique weak state solution to (1), and this solution must be non-negative and bounded above. The corresponding two theorems follow from results as in [13,14].

**Theorem 1** (Positivity and Existence of State Solution)**.** *Given N*<sup>0</sup> *non-negative, L*<sup>∞</sup> (0, *L*) *bounded and in H*<sup>1</sup> {0} (0, *<sup>L</sup>*)*, then, for each <sup>Q</sup>* <sup>∈</sup> *U, there is a unique non-negative weak solution <sup>N</sup>* <sup>=</sup> *<sup>N</sup>* (*Q*) *of tjhe state PDE (1) with* 0 ≤ *N*(*x*, *t*) ≤ *C, where C holds for all Q* ∈ *U.*

**Theorem 2** (A priori estimates)**.** *There exists positive constants α*<sup>1</sup> *and α*<sup>2</sup> *such that, for any Q* ∈ *U*, *and <sup>N</sup>* <sup>∈</sup> *<sup>V</sup>* <sup>∩</sup> *<sup>L</sup>*∞(Λ) *where Nt* <sup>∈</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*;(*H*<sup>1</sup> {0}(0, *<sup>L</sup>*))∗)*, the weak solution of the PDE corresponding to <sup>Q</sup> with N* ≥ 0 *a.e. satisfies*

$$\left(\iint\_{\Lambda} N\_x^2 dxdt + \iint\_{\Lambda} N^2 dxdt\right) \le a\_1 \tag{2}$$

*and NtV*<sup>∗</sup> ≤ *α*2.

We denote the weak solution of the state PDE (1) corresponding to control *Q* as *N*(*Q*).

**Theorem 3** (Existence of Optimal Control)**.** *There exists an optimal control, Q*<sup>∗</sup> ∈ *U*, *satisfying*

$$J(Q^\*) = \inf\_{Q \in \mathcal{U}} J(Q).$$

**Proof.** We choose a minimizing sequence {*Qn*} ⊂ *<sup>U</sup>* with *<sup>n</sup>* = 1, 2, 3, . . . such that

$$\lim\_{n \to \infty} J\left(Q^n\right) = \inf\_{Q \in U} J(Q).$$

and *N<sup>n</sup>* = *N* (*Qn*) the corresponding solution to tje state PDE (1). By the a priori estimates in Theorem 2, there exists *<sup>N</sup>*<sup>∗</sup> <sup>∈</sup> *<sup>V</sup> <sup>L</sup>*<sup>∞</sup> ((0, *<sup>T</sup>*) <sup>×</sup> (0, *<sup>L</sup>*)), with *<sup>N</sup>*<sup>∗</sup> *<sup>t</sup>* ∈ *V*∗, and *Q*<sup>∗</sup> ∈ *U* such that, on a subsequence, these weak convergences hold:

$$\begin{aligned} N^{\mathfrak{n}} &\rightharpoonup N^\* \text{ in } V, \\ N\_t^{\mathfrak{n}} &\rightharpoonup N\_t^\* \text{ in } V^\*, \\ Q^{\mathfrak{n}} &\rightharpoonup Q^\* \text{ in } L^2 \left( (0, T) \times (0, L) \right). \end{aligned} \tag{3}$$

We must show that *N*∗ = *N* (*Q*∗), where *N*∗ is the population solution associated with *Q*∗. Using the compactness result in [15] coupled with the convergences in (3), we obtain this strong *L*<sup>2</sup> convergence,

$$N^n \to N^\* \text{ in } L^2\left((0, T) \times (0, L)\right) \dots$$

Now, we must show that *Q*∗ is an optimal control and that *N*∗ = *N* (*Q*∗). Using the convergences in (3) and the *L*<sup>∞</sup> bounds on *Q*, the convergences above give limits for these terms:

$$\begin{aligned} &\int\_{0}^{T} \langle N\_{t}^{\mathrm{H}}, \phi \rangle dt \to \int\_{0}^{T} \langle N^{\*}, \phi \rangle dt \\ &\int\_{0}^{T} \int\_{0}^{L} \frac{A\_{t}}{A} N^{\mathrm{u}} \phi dx dt \to \int\_{0}^{T} \int\_{0}^{L} \frac{A\_{t}}{A} N^{\*} \phi dx dt \\ &\int\_{0}^{T} \int\_{0}^{L} DAN\_{x}^{\mathrm{u}} \left( \frac{\phi}{A} \right)\_{x} dx dt \to \int\_{0}^{T} \int\_{0}^{L} DAN\_{x}^{\*} \left( \frac{\phi}{A} \right)\_{x} dx dt \\ &\int\_{0}^{T} \int\_{0}^{L} -r N^{\mathrm{u}} \phi dx dt \to \int\_{0}^{T} \int\_{0}^{L} -r N^{\*} \phi dx dt \\ &\int\_{0}^{T} \int\_{0}^{L} \frac{r \Phi}{K} \left( N^{n} \right)^{2} dx dt \to \int\_{0}^{T} \int\_{0}^{L} \frac{r \Phi}{K} \left( N^{\*} \right)^{2} dx dt. \end{aligned}$$

To handle a term with *Q<sup>n</sup>* and *N<sup>n</sup> <sup>x</sup>* , integrating by parts and rearranging gives

$$\begin{split} \int\_{0}^{T} \int\_{0}^{L} \frac{1}{A} \left( Q^{n} N\_{x}^{n} - Q^{\*} N\_{x}^{\*} \right) \boldsymbol{\Phi} \mathbf{dx} dt \\ &= - \int\_{0}^{T} \int\_{0}^{L} Q^{n} \left( \frac{\boldsymbol{\Phi}}{A} \right)\_{x} \left( N^{n} - N^{\*} \right) d\mathbf{x} dt - \int\_{0}^{T} \int\_{0}^{L} N^{\*} \left( \frac{\boldsymbol{\Phi}}{A} \right)\_{x} \left( Q^{n} - Q^{\*} \right) d\mathbf{x} dt \\ &\quad + \int\_{0}^{T} \frac{Q^{n} \boldsymbol{\Phi}}{A} \left( N^{n} - N^{\*} \right) \left( L, t \right) dt + \int\_{0}^{T} \frac{N^{\*} \boldsymbol{\Phi}}{A} \left( Q^{n} - Q^{\*} \right) \left( L, t \right) dt. \end{split} \tag{4}$$

By the Cauchy–Schwartz inequality, the strong convergence of {*Nn*}, *<sup>φ</sup>* being in *<sup>L</sup>*<sup>2</sup> and the *<sup>L</sup>*<sup>∞</sup> bounds on {*Qn*}, the first and third terms of Equation (4) give

$$\begin{split} \lim\_{n \to \infty} \left| -\int\_{0}^{T} \int\_{0}^{L} Q^{n} \left( \frac{\phi}{A} \right)\_{x} \left( N^{n} - N^{\*} \right) dxdt + \int\_{0}^{T} \frac{Q^{n} \phi}{A} \left( N^{n} - N^{\*} \right) \left( L, t \right) dt \right| \\ \leq \lim\_{n \to \infty} M \left[ \left( \int\_{0}^{T} \int\_{0}^{L} |N^{n} - N^{\*}|^{2} dxdt \right)^{1/2} \left( \int\_{0}^{T} \int\_{0}^{L} \left| \frac{\phi}{A} \right|\_{x}^{2} dxdt \right)^{1/2} \\ + \left( \int\_{0}^{T} \left| \frac{\phi}{A} (L, t) \right|^{2} dt \right)^{1/2} \left( \int\_{0}^{T} |N^{n} - N^{\*}|^{2} \left( L, t \right) dt \right)^{1/2} \right] \\ = 0 \end{split}$$

using the bounds in (2) and the convergences in (3) and a result from [16] for

$$N^n \to N^\* \text{ in } L^2\left(\{L\} \times (0, T)\right) \dots$$

By the weak convergence of {*Qn*}, *<sup>φ</sup>* being fixed in *<sup>L</sup>*<sup>2</sup> and *<sup>N</sup>*<sup>∗</sup> being fixed with *<sup>L</sup>*<sup>∞</sup> bounds, the second and fourth terms of Equation (4) converge, which gives

$$\lim\_{n \to \infty} \int\_0^T \int\_0^L \frac{1}{A} \left( N\_x^n Q^n - N\_x^\* Q^\* \right) \phi dx dt = 0.$$

Hence, as *n* → ∞, *N*<sup>∗</sup> satisfies the state PDE with the control *Q*<sup>∗</sup> in (1), i.e., *N*<sup>∗</sup> = *N* (*Q*∗). Again using the weak convergences (3) in the objective functional on the minimizing sequence {*Qn*} with *<sup>L</sup>*<sup>∞</sup> bounds and with lower semicontinuity with respect to *L*<sup>2</sup> weak convergence [17],

$$\int\_0^T \left(Q^\*\right)^2 dt \le \varliminf\_{\mathfrak{u}} \int\_0^T \left(Q^{\mathfrak{u}}\right)^2 dt.$$

we obtain inf*Q*∈*<sup>U</sup> J* (*Q*) ≥ *J* (*Q*∗) and then *Q*<sup>∗</sup> is an optimal control.

We now differentiate the control-to-state map *Q* → *N*(*Q*) as a directional derivative to find the sensitivity function.

**Theorem 4** (Differentiability of solution map)**.** *The solution map*

$$Q \in \mathcal{U} \to \mathcal{N} = \mathcal{N}(\mathcal{Q}) \in V$$

*is differentiable in the following sense: There exists ψ* ∈ *V, such that*

$$\frac{N(Q+\epsilon l) - N(Q)}{\epsilon} \to \psi$$

$$\begin{array}{llll} \text{as } \boldsymbol{\epsilon} \to 0, \text{ where } \boldsymbol{Q} + \boldsymbol{\epsilon}\boldsymbol{l} \in \boldsymbol{U} \text{ and } \boldsymbol{l} \in \boldsymbol{L}^{\infty}(0,T). & \text{Moreover, } \boldsymbol{\psi} = \boldsymbol{\psi}(\boldsymbol{Q},\boldsymbol{l}) \text{ with } \boldsymbol{\psi}\_{t} \in \boldsymbol{V}^{\*} \text{ satisfies} \\ \boldsymbol{\psi}\_{t} = \frac{-A\_{t}}{A}\boldsymbol{\psi} + \frac{1}{A}\left(DA\boldsymbol{\psi}\_{\boldsymbol{x}}\right)\_{\boldsymbol{x}} - \frac{\boldsymbol{Q}}{A}\boldsymbol{\psi}\_{\boldsymbol{x}} - \frac{1}{A}\mathbf{N}\_{\boldsymbol{x}} + r\boldsymbol{\psi} - \frac{2r}{K}\boldsymbol{\psi}\boldsymbol{N} & \text{in} \quad \Lambda = (0,L) \times (0,T), \\ \boldsymbol{\psi}(0,t) = 0 & \text{on} \quad (0,T), \boldsymbol{x} = 0, \\ \boldsymbol{\psi}\_{t}(L,t) = 0 & \text{on} \quad (0,T), \boldsymbol{x} = L, \\ \boldsymbol{\psi}(\boldsymbol{x},0) = 0 & \text{on} \quad (0,T), \boldsymbol{t} = 0. \end{array}$$

**Proof.** To justify the convergence of the quotients, we use a priori estimates. Thus, we change the variables such that *w*(*x*, *t*) = *N*(*x*, *t*)*e*−*α<sup>t</sup>* and *w* (*x*, *t*) = *N* (*x*, *t*)*e*−*α<sup>t</sup>* , where *N* = *N*(*Q*) and *N*  = *N*(*Q* +  *l*). With some standard PDE estimates, we obtain the desired estimate for *α* large

$$\frac{1}{2} \int\_{0}^{L} \left(\frac{w^{\varepsilon} - w}{\varepsilon}\right)^{2} (T, x) dx + (a - \mathbb{C}\_{1}) \iint\_{\Lambda} \left(\frac{w^{\varepsilon} - w}{\varepsilon}\right)^{2} dxdt + \frac{\mathfrak{c}\_{D}}{2} \iint\_{\Lambda} \left(\frac{w^{\varepsilon} - w}{\varepsilon}\right)^{2}\_{x} dxdt \le \mathbb{C}\_{0}.$$

The bounds on *<sup>w</sup>* <sup>−</sup>*<sup>w</sup>*  give corresponding bounds on *<sup>N</sup>* <sup>−</sup>*<sup>N</sup>*  . From its PDE, we obtain a uniform estimate on *<sup>N</sup> <sup>t</sup>* −*Nt*  in *<sup>V</sup>*∗. Thus, there exists *<sup>ψ</sup>* <sup>∈</sup> *<sup>V</sup>* and *<sup>ψ</sup><sup>t</sup>* <sup>∈</sup> *<sup>V</sup>*<sup>∗</sup> such that *<sup>N</sup>*(*Q*<sup>+</sup> *<sup>l</sup>*)−*N*(*Q*)  *ψ* in *V* and *Nt*(*Q*+ *l*)−*Nt*(*Q*)  *ψ<sup>t</sup>* in *V*<sup>∗</sup> as  → 0.

Using weak convergence of the quotients, *ψ* satisfies the sensitivity PDE, by passing to the limit in

$$\begin{split} \left(\frac{N^{\varepsilon}-N}{\varepsilon}\right)\_{t} &= \frac{-A\_{t}}{A} \left(\frac{N^{\varepsilon}-N}{\varepsilon}\right) + \frac{1}{A} \left(DA\left(\frac{N^{\varepsilon}-N}{\varepsilon}\right)\_{x}\right)\_{x} \\ &- \frac{Q}{A} \left(\frac{N^{\varepsilon}-N}{\varepsilon}\right)\_{x} - \frac{\varepsilon l}{A\varepsilon}N\_{x}^{\varepsilon} + r\left(\frac{N^{\varepsilon}-N}{\varepsilon}\right) - \frac{r}{K}\left(\frac{N^{\varepsilon}-N}{\varepsilon}(N^{\varepsilon}+N)\right). \end{split}$$

Since *N*(*Q* +  *l*) → *N*(*Q*) strongly in *V*, we obtain

$$\psi\_t = \frac{-A\_t}{A}\psi + \frac{1}{A}\left(DA\psi\_x\right)\_x - \frac{Q}{A}\psi\_x - \frac{l}{A}N\_x + r\psi - \frac{2r}{K}\psi N.$$

Our initial and boundary conditions become *ψ*(0, *t*) = 0, *ψx*(*L*, *t*) = 0, and *ψ*(*x*, 0) = 0.

Next, we build our adjoint PDE using the sensitivity PDE. We arrange the sensitivity equation corresponding to an optimal state *N*∗ = *N*(*Q*∗) with direction *l* above as

$$\psi\_t + \frac{A\_t}{A}\psi - \frac{1}{A}\left(DA\psi\_x\right)\_x + \frac{Q^\*}{A}\psi\_x - r\psi + \frac{2r}{K}\psi N^\* = -\frac{l}{A}N\_x^\*$$

to obtain the linear operator on *ψ*,

$$L\psi = \psi\_t + \frac{A\_t}{A}\psi - \frac{1}{A}\left(DA\psi\_x\right)\_x + \frac{Q^\*}{A}\psi\_x - r\psi + \frac{2r}{K}\psi N^\*\psi$$

and thus the sensitivity PDE can be written as *<sup>L</sup><sup>ψ</sup>* <sup>=</sup> <sup>−</sup> *<sup>l</sup> <sup>A</sup> N*<sup>∗</sup> *x* .

For our adjoint function, we wish to find an operator *L*<sup>∗</sup> such that *λ*, *Lψ* = *L*∗*λ*, *ψ* in *V*, *V*∗ duality. Starting with the weak form of the LHS of *ψ* PDE and test function *λ* with *λ* (0, *t*) = 0 and using integration by parts, to obtain *λ*, *Lψ* = *L*∗*λ*, *ψ*, our adjoint boundary condition for *x* = *L* becomes

$$D(L,t)A(L,t)\left(\frac{\lambda(L,t)}{A(L,t)}\right)\_x + Q(t)\frac{\lambda(L,t)}{A(L,t)} = 0.5$$

Thus, we found that

$$\begin{cases} L^\*\lambda = -\lambda\_l - \left(DA\left(\frac{\lambda}{A}\right)\_x\right)\_x - Q\left(\frac{\lambda}{A}\right)\_x + \frac{A\_l\lambda}{A} - r\lambda + 2\frac{r\lambda}{X}N^\* & \text{in} \quad \Lambda\_r\\ L^\*\lambda = \mathcal{W}(\mathbf{x}) & \text{in} \quad \Lambda\_r\\ \lambda(\mathbf{x}, T) = 0 & \text{on} \quad (0, L), t = T, \\ \lambda(0, t) = 0 & \text{on} \quad (0, T), \mathbf{x} = 0, \\ D(L, t)A(L, t)\left(\frac{\lambda(L, t)}{A(L, t)}\right)\_X + Q(t)\frac{\lambda(L, t)}{A(L, t)} = 0 & \text{on} \quad (0, T), \mathbf{x} = L. \end{cases} \tag{5}$$

We choose the RHS of the adjoint PDE such that *L*∗*λ* = *W*(*x*), where *W*(*x*) is the derivative of the integrand of the objective functional with respect to *N*(*x*, *t*). We use the sensitivity and adjoint PDEs in the differentiation of map *Q* → *J*(*Q*) to obtain our optimal control characterization.

**Theorem 5.** *Given an optimal control <sup>Q</sup>*<sup>∗</sup> *in U, there exists a solution <sup>λ</sup> in <sup>V</sup>* <sup>∩</sup> *<sup>L</sup>*<sup>∞</sup> (Λ) *with <sup>λ</sup><sup>t</sup>* <sup>∈</sup> *<sup>V</sup>*<sup>∗</sup> *to the adjoint problem (5). Moreover, we obtain a characterization of our optimal control:*

$$Q^\*(t) = \min\left(M, \max\left(\frac{1}{2\epsilon} \int\_0^L \frac{\lambda}{A} N\_x^\*(x, t) dx, m\right)\right),\tag{6}$$

*where*  > 0*.*

**Proof.** Given the existence of an optimal control, the corresponding solution of the adjoint problem exists since the adjoint problem is linear in the adjoint function [18].

Suppose *<sup>Q</sup>*∗(*t*) is an optimal control. Let *<sup>l</sup>*(*t*) <sup>∈</sup> *<sup>L</sup>*∞(0, *<sup>T</sup>*) such that *<sup>Q</sup>* <sup>+</sup> *<sup>δ</sup><sup>l</sup>* <sup>∈</sup> *<sup>U</sup>* for *<sup>δ</sup>* <sup>&</sup>gt; 0. We denote *N*(*Q* + *δl*) by *Nδ*. The derivative of *J*(*Q*) with respect to *Q* at *Q*<sup>∗</sup> in the direction *l* satisfies

0 ≤ lim *δ*→0 *J*(*Q*<sup>∗</sup> + *δl*) − *J*(*Q*∗) *δ* = lim *δ*→0 1 *δ T* 0 *L* 0 *W*(*x*)*Nδ*(*x*, *t*)*dxdt* + *T* 0  (*Q*∗ + *δl*) <sup>2</sup> (*t*)*dt* − *T* 0 *L* 0 *W*(*x*)*N*∗(*x*, *t*)*dxdt* + *T* 0  (*Q*∗) <sup>2</sup> (*t*)*dt* = *T* 0 *L* 0 *W*(*x*)*ψdxdt* + *T* 0 2*l Q*∗*dt* = *T* 0 *L* 0 −*λtψ* − *DA <sup>λ</sup> A x ψx* − *Qψ λ A x* + *At <sup>A</sup> λψ* <sup>−</sup> *<sup>r</sup>λψ* <sup>+</sup> <sup>2</sup> *<sup>r</sup> <sup>K</sup> <sup>N</sup>*∗*λψ dxdt* + *T* 0 2*l Q*∗*dt* = *T* 0 *L* 0 *ψtλ* + *At <sup>A</sup> λψ* <sup>−</sup> *DA <sup>λ</sup> A x ψx* <sup>+</sup> *<sup>Q</sup> <sup>A</sup> <sup>ψ</sup>x<sup>λ</sup>* <sup>−</sup> *<sup>r</sup>ψλ* <sup>+</sup> 2*r <sup>K</sup> ψλN*<sup>∗</sup> *dxdt* + *T* 0 2*l Q*∗*dt* = *T* 0 *l* 2 *Q*<sup>∗</sup> − *L* 0 *λ <sup>A</sup> <sup>N</sup>*<sup>∗</sup> *<sup>x</sup> dx dt*.

Choosing appropriate variations for *l* gives the desired characterization of our optimal control.

The optimality system consists of the state and adjoint systems (1), (5), and the optimal control characterization (6). The uniqueness of the optimal control follows from the uniqueness of the solutions to the optimality system. Techniques that were used in [13,19] and this extra estimate *L* <sup>0</sup> *<sup>p</sup>*¯<sup>2</sup> *xdx* <sup>≤</sup> *<sup>K</sup>*5*e*2*α<sup>T</sup>* [20] gives the desired uniqueness result.

**Theorem 6** (Uniqueness of Optimal Control)**.** *For sufficiently small T, the solution to the optimality system (1), (5), and (6) is unique.*

If one has a valid generalization of this model to two of three spatial dimensions, these theoretical results for the optimal control analysis and the corresponding optimality system could be generalized, with appropriate a priori estimates . The assumptions on the diffusion coefficient and the cross-sectional area can be generalized as long as the needed a priori estimates and corresponding compactness results can be obtained. For some background on such optimal control problems, see the book by Li and Yong [21].

#### **3. Numerical Results**

In the previous sections, we have that a unique optimal control exists for our optimality system. Now, we illustrate numerical solutions for the optimality system for several scenarios. To solve for an optimal control numerically, we must solve the optimality system consisting of the state and adjoint PDEs with the optimal control characterization. For this, we use the Forward-Backward Sweep method [22,23]. We make an initial guess for the control, *Q*. Then, using the initial condition, we solve the state PDE, *N*, forward in time. Next, we solve the adjoint PDE, *λ*, backward in time. The PDEs are solved using finite differences as discussed below. We update *Q* using our calculated *N* and *λ* values in the optimal control characterization, and, finally, we check convergence of control values at successive iterations. A discussion about convergence and stability of this method can be found in [22]. In the numerical approximations, we use an upwind or downwind difference as appropriate. For an example, in our state PDE for the advection term, *Nx*, we use a downwind difference [24] since its coefficient, <sup>−</sup> *<sup>Q</sup>*(*t*) *<sup>A</sup>*(*x*,*t*), is negative. This difference goes backwards in space. Similarly, we use a downwind difference for the *At* term as well. We use an upwind difference in the adjoint advection term since *Q*(*t*) > 0, which goes forward in space [24].

#### *3.1. Cross-Sectional Area Is Constant*

We begin these numerical simulations assuming our cross-sectional area, *A*, is constant. Later, we will vary it with respect to space and time. To illustrate the spread of a population in a river with control, we use the following Base Case parameters:

$$\begin{aligned} T &= 10 \text{ (final time)}, L = 10 \text{ (length of river)}, \\ r &= 0.6, K = 200, D = 0.1, A = 20, \varepsilon = 0.05, \, 0 \le Q(t) \le 10. \end{aligned} \tag{7}$$

Note that the downstream is represented by *x* = *L*. Figure 1 shows the initial condition and the weight function used in our simulations. In Figure 1a, we place the initial population downstream to see their movement upstream. In Figure 1b, we use a weight function that has values 100 upstream (at *x* = 0) and nearly 0 upstream (at *x* = 10). Recall that we are minimizing the objective functional, so this weight function helps push the population downstream.

We first illustrate how far the population moves upstream with no control, constant control, and optimal control in Figure 2. To generate this graph, we use a lower bound of 0.5 for the population, so, for each time, we find the lowest location *x* (more upstream) with *N*(*x*, *t*) > 0.5 and record that river location. The constant control has the value of the upper bound of the optimal control. The population with no control moves further upstream than the population with constant control, which is to be expected. The population with no control moves from 8 to about 2.8, which is about two-thirds upstream. The population with constant control starts at 8 and ends about 7.1. Having a constant control at the upper bound (*Q* = 10) resulted in the least movement upstream. With the optimal control, the population moves further upstream than the constant control case. Note that the optimal control is not at its upper bound for the whole time period due to costs in *J*. The objective functional

balances keeping the population downstream with costs. The population in the optimal control case starts at 8 and ends around 5.5.

**Figure 1.** (**a**) The initial condition for the population; (**b**) the weight function in the objective functional.

**Figure 2.** Location of the population in the river for the no control case (red), for the constant control case with *Q* = 10 (blue), and the optimal control case (magenta) over the time, with parameter values fromtje Base Case list (7).

Figure 3 shows the population plots for no control (Figure 3a) and for optimal control (Figure 3b). The population in the no control case spreads farther upstream than in the optimal control case. The population also grows larger in the no control case. In Figure 3c, the populations at the final time are plotted for the length of the river. The no control population is further upstream and is approaching the carrying capacity while the optimal control population only reached a value of about 140 at any given river location. The population size was decreased, and the population was kept more downstream in the optimal control case. Our objective functional value for the no control case is 239.57, while, for the optimal control case, is 40.12, which is an 83% improvement.

**Figure 3.** Population plots for (**a**) no control case and (**b**) optimal control case; (**c**) population plots in *x* at the final time for both cases. The parameter values are from the Base Case list (7).

We varied parameters for the no control, constant control, and optimal control cases and show the changes in our objective functional *J* values. Starting with our Base Case parameter values in (7), we changed the values of *D*, *T*, *r* and *K* and investigated the resulting changes in *J*. Table 1 shows the *J* values for the Base Case along with cases with changes of *D* and *T* while Table 2 shows the variation of *K* and *r*. In all cases, the objective functional from the optimal control was smaller than the objective functional for the constant control and the no control cases, which is expected. In Figure 4, we plot the optimal controls for each of the parameters varied. Before changing any parameters for the Base Case, the optimal control is a 29% decrease in the objective functional from the constant control. This indicates that a time varying control can be significantly effective.

In Table 1, for the decrease in diffusion coefficient *D*, we see a 45% decrease in the objective functional from the constant control to the optimal control and a 74% decrease from the no control to optimal control. When we increase diffusion, there is almost no change in the *J* values for the constant control and optimal control. There is a 94% decrease from the objective functional value from the no control to the optimal control *J* value. Figure 4a shows the different optimal controls. When we increase the diffusion coefficient, the optimal control is at the upper bound for most of the time, which explains why the objective functional value was similar to the constant control one. Decreasing the diffusion creates a similar curve to the Base Case, but the optimal control values are slightly smaller. Returning to Table 1, decreasing the final time yields a 57% decrease in the objective functional from the constant control population to the optimal control population, and a 28% decrease from the no control to the optimal control. Meanwhile, increasing the final time gives a 20% decrease again from

the constant control population to the optimal control population and a 97% decrease from no control to optimal control. In Figure 4b, increasing *T* causes the optimal control to be at the upper bound for about half the time, and to then steadily decrease to 0 at the final time. If we decrease *T*, the optimal control starts midway between the upper and lower bound and decreases down to 0 at the final time. The Base Case falls directly between these two cases.

**Table 1.** Objective functional outputs for the cases tested where we changed the parameter values for *D* and *T* (given below). The base values are from the list (7). Constant control was run with *Q* = 10.


Table 2 shows a 32% decrease in the objective functional from the constant control to the optimal control when we decrease the carrying capacity, *K*. Between the no control population and the optimal control population, there is a 81% decrease in the optimal control objective functional value. Increasing *K* yields a 28% decrease of *J* from the constant control to the optimal control, while a 85% decrease in *J* is seen from the no control case to the optimal control case. Figure 4c illustrates that increasing and decreasing *K* has little change in the optimal control values.

For decreasing the growth rate, *r*, there is a 61% reduction in the objective functional from the constant control population to the optimal control population. A decrease of 63% occurs in the objective functional values from the no control case to the optimal control case. When we increased *r*, there is virtually no change in the objective functional values between the constant control and the optimal control, but there was a 95% reduction in the *J* value from the no control population to the optimal control population. Figure 4d shows the optimal control at the upper bound for most of the time when we increase *r*. Decreasing *r* gives the optimal control initially at 6 and steadily decreases to 0 at the final time.

Since our objective functional is meant to represent minimizing the population at upstream and the cost of removing the invasive species, in most cases, using an optimal control would give at least a 20% decrease in *J*(*Q*∗) from the constant control. In addition, the optimal control would give a smaller objective functional values *J*(*Q*∗) than applying no control at all.


**Table 2.** Objective functional outputs for the cases tested where we changed parameter values for *K* and *r* (given below). The base values are from the list (7). Constant control was run with *Q* = 10.

**Figure 4.** The optimal control plots for varying of the parameters *D*, *T*, *K*, and *r*. The base parameter values are from the list (7).

#### *3.2. Cross-Sectional Area Is Not Constant*

We now have the cross-sectional area dependent on space with a linear function,

$$A(\mathbf{x}) = 0.5\mathbf{x} + 25\mathbf{y}$$

and ran all the same cases (varying *D*, *T*, *K*, and *r*) as for a constant cross-sectional area. The parameters are the same as before in the Base Case list (7). We do not show the specific graphs for this case, but we will compare the various location results later.

Finally, we let the cross-sectional area vary with time and space where

$$A(\mathbf{x}, t) = (0.5\mathbf{x} + 2\mathbf{5}) + (0.2t\,(10 - t))\,\mathrm{J}$$

We chose a function so that the river was higher in the middle of the year. Again, the parameters are all the same as in the Base Case list (7). In Figure 5, the no control population starts at 8 and moves to 2.8 with an objective functional value of 238.66. The population with constant control moves from 8 to 5.5 with an objective functional value of 74.70. The population with the optimal control starts at 8 and ends at 4.9. The objective functional value for the optimal control is 66.55, which is again a 10.9% decrease from the constant control and a 72.1% decrease from the no control case. Figure 6 shows the population plots for the no control and optimal control cases, and, as expected, the population in the no control case is larger.

**Figure 5.** Location of the population in the river for the no control case (red), for the constant control case with *Q* = 10 (blue), and the optimal control case (magenta) over the time with the cross-sectional area is dependent on space and time, *A*(*x*, *t*) = (0.5*x* + 25) + (0.2*t*(10 − *t*)). The parameter values are *T* = 10, *L* = 10, *r* = 0.6, *K* = 200, *D* = 0.1, and 0 ≤ *Q*(*t*) ≤ 10.

**Figure 6.** The population plots for the (**a**) no control case and (**b**) the optimal control case with the cross-sectional area is dependent on space and time, *A*(*x*, *t*) = (0.5*x* + 25) + (0.2*t*(10 − *t*)). The parameter values are from the Base Case list (7).

We compare the optimal control results for the three different cross-sectional area functions. We have that 25 ≤ *A*(*x*) ≤ 35 and 25 ≤ *A*(*x*, *t*) ≤ 35, while the constant *A* = 20. In Figure 7, in the constant cross-sectional area case, the population moved upstream the least, at 5.5 at the furthest point. In addition, notice that *A*(*x*) and *A*(*x*, *t*) caused similar movement upstream. When the cross-sectional area is not constant, the species is able to move further upstream, possibly due to larger cross-sectional

area. The objective functional values are 40.12 for constant *A*, 60.55 for *A*(*x*), and 66.55 for *A*(*x*, *t*). The *J* values increase as our cross-sectional areas become larger and more complex. This makes sense as the effort to keep the population downstream would become more difficult with the varying larger areas.

**Figure 7.** Location of the population in the river for the optimal control where the different cross-sectional area values are: constant *A* = 20 (blue), *A*(*x*) = 0.5*x* + 25 (magenta), and *A*(*x*, *t*) = (0.5*x* + 25) + (0.2*t*(10 − *t*)) (red). The parameter values are from the Base Case list (7).

#### *3.3. Approximate Controls*

Optimal controls may change at every time step. However, changing river flow at every time step may be difficult, so we may use an optimal control to design an approximate control that is easier to implement. For example, using the optimal control as a guide, we design an approximate control where we set the control at a high value at the beginning time and a low value towards the end of our time period. In Figure 8a, we set the control, *Q*, to 10 for *t* ≤ 5 and then *Q* = 0 for *t* > 5, and plot this approximate control with the optimal control. For the approximate control, our objective functional value is 48.01, which is an increase from the 40.12 from the optimal control. In Figure 8b, another approximate control has *Q* = 8 for *t* ≤ 7 and *Q* = 0 for *t* > 7. Here, our objective value is 42.49, which is again an increase compared to the Base Case value from the optimal control. This approximate control in Figure 8b is a good approximation since its corresponding *J* value is only 6% different from the optimal *J*.

We also changed the weight *W* in the objective function to be linear function and obtained quantitatively similar results. Changing the initial condition to be somewhat larger gives a larger objective functional values but with a similar optimal control.

**Figure 8.** Control plots for both the optimal control case and two constant case. (**a**) the two constant control case where *Q* = 10 for *t* ≤ 5 and *Q* = 0 for *t* > 5, and (**b**) the two constant control case where *Q* = 8 for *t* ≤ 7 and *Q* = 0 for *t* > 7. The parameter values are from the Base Case list (7).

#### **4. Conclusions**

To investigate how an invasive population moves upstream, we began with a PDE describing the movement of a species upstream, and then completed the needed analysis results for optimal control of this PDE using the flow rate as the control. We proved the differentiability of the control-to-state map as a directional derivative, which satisfies a linear sensitivity PDE. From the sensitivity PDE, we constructed our adjoint PDE that we used to characterize our optimal control. The control was chosen to keep the population downstream while minimizing the implementation cost.

To run numerical simulations to solve our optimality system, we used a finite difference method together with the forward-backward sweep method. In our simulations, we explored how changing the parameter values would affect the population's location in the river and the *J* values for cases with no control, constant control, and optimal control while leaving our initial conditions the same for all cases and parameter changes. We first did this for a constant cross-sectional area of the river. As expected, we saw the optimal control populations were always an improvement in *J* on the no control populations and the constant control populations. The time varying control allowed for an immediate response to the changing population size. We also observed that varying the final time greatly changed the population's location and size.

Changing the cross-sectional area to be dependent on space, and later on space and time, increased the *J* values from the constant cross-sectional area case with *A*(*x*, *t*) being the largest. The non-constant *A* allowed the populations with constant and optimal controls to move further upstream than before (as seen in Figures 2, 5, and 7). The population probably moved farther upstream since there is more room for the species to move in the river for the *A*(*x*) and *A*(*x*, *t*) cases due to larger sizes of the cross-sections.

We investigated how our objective functional values changed if we implemented approximate controls that are easily implemented. We constructed an approximate control that gave a *J* value close to the optimal value.

Using optimal control theory can give strategies for controlling a flow rate to successfully manage a population in a river. Our future work will be to obtain some data for a particular invasive river species and apply these techniques. After approximating parameters and choosing appropriate cross-sectional area for a particular species in a river, we would explore how our optimal control strategy would affect the movements of the invasive population. We would also want to investigate more realistic cross-sectional area functions to better represent a river. Considering only allowing seasonal changes in controls would also be important for some situations.

**Author Contributions:** Investigation, analysis, computations and writing, S.L. and R.P.

**Funding:** This research received no external funding.

**Acknowledgments:** We acknowledge useful conservations with Mark Lewis and Yu Jin. This work was partially supported by the National Institute for Mathematical and Biological Synthesis, an Institute sponsored by the National Science Foundation through NSF Award #DBI-1300426, with additional support from the University of Tennessee, Knoxville.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
