- Email: [email protected]
A control problem in a polluted environment
E~E~ER
A Control Problem in a Polluted Environment DIANA MARIA THOMAS School of Mathematics, Georgia Tech, Atlanta, Georgia 30332 TERRY WAYNE SNELL School of Biology, Georgia Tech, Atlanta, Georgia 30332 AND
SAYED MURTAZA JAFFAR School of Mechanical Engineering, Georgia Tech, Atlanta, Georgia 30332
Received 29 August 1994; revised 16 June 1995
ABSTRACT Due to the importance of remediating polluted environments, many mathematical models have been developed to describe population-toxicant interactions. This article considers the situation where the amount of pollution in the environment is restricted to lie below a certain value M. A model is presented and viewed as a control problem in which the toxicant input into the environment is the controllable variable. An explicit formula is computed for the input function that meets certain restrictions. By analytically studying the dynamics of the population, it is possible to provide conditions on how to set M in order to guarantee survival of the population without a significant reduction of its carrying capacity. Numerical and experimental data are described that support these results.
1.
INTRODUCTION
Efforts to clean up polluted ecosystems have led to many regulations on amounts of toxicants allowed in the environment. This includes water quality criteria set by the Environmental Protection Agency (EPA) to limit certain chemicals in freshwater environments. Establishment of regulations depend on the type of behavior desired in the system. Generally, the objective is to prevent species loss and maintain ecosystem function. The aim of this paper is to study the dynamics of a single-species population in a polluted freshwater environment where the input of pollution is controlled. We focus on a specific example using rotifers, an important component of the zooplankton in many aquatic communities.
MATHEMATICAL BIOSC1ENCES 133:139-163 (1996) 0025-5564/96/$15.00 © Elsevier Science Inc., 1996 SSDI 0025-5564(95)00091-Q 655 Avenue of the Americas, New York, NY 10010
140
DIANA MARIA THOMAS ET AL.
Section 2 describes biological aspects of rotifers and presents some experimental results. In Section 3 a model is presented that describes population-toxicant interactions. The control problem is posed and solved in Section 4, and several properties of the model are consequently derived. In Section 5 we discuss the dynamical aspects of the model by conducting a stability analysis. Conditions are consequently derived that determine how to set the restrictions on the amount of toxicant in the environment in order to guarantee survival of the population without a significant reduction of its carrying capacity. In addition, some numerical results are presented. Section 6 summarizes the results of the work and poses future problems. 2. THE BIOLOGICAL EXAMPLE AND METHODS Rotifers are multicellular aquatic animals that eat single-celled algae in lakes, ponds, rivers, and wetlands [1]. Although they are small (0.1-1 mm in length), they can occur in high densities and modify the composition of phytoplankton populations by selective grazing. Most rotifers are herbivores providing a critical link in aquatic food webs between microscopic plants and their fish predators. The dynamics of laboratory and field populations of rotifers have been characterized in many ecological studies [1]. Rotifer population dynamics also have been used to examine the effects of pollution [2-4], with population growth are, r, the usual criterion. Janssen et al. illustrated how carrying capacity, K, could be a useful criterion in continuous cultures of the rotifer Brachionus calyciflorus [5]. They showed that 30-day exposures (approximately 15 generations) to 10 /zg copper/L reduced K to about half of the control, which is the unexposed culture (Figure 1). The copper concentration causing significant reduction in K is about half that of LCs0 [5]. Therefore, the effect of pollution on rotifer populations is to lower both r and K. The objective of environmental regulation is to set limits on how much pollution can be allowed into the environment without significantly reducing indigenous populations. The model that we have developed provides guidance about the maximum amount of pollution allowable before substantially lowering the carrying capacity. Using a system of differential equations based on the logistic population growth model, we derive a relationship between environmental concentrations of pollution and population dynamics. 3. THE MODEL AND PRELIMINARIES In this section, we present a modification of the model suggested in [6] that displays logistic population growth coupled with toxicant stress.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
12o--]
141
Brachionus calyciflorus _
c o p p e r L C 5 0 = 2 6 lag/L
controJ~
100t 807
I
2.5 pg Cu
~
_ _
o.ns,,, .I v
v/--I
a
I
l
I
l
l
I
I
-
i
i
l
J
l
l
l
J
l
l
l
l
l
i
J
l
l
l
l
l
I 2 3 4 5 6 7 8 9 1011 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0 2 1 2 2 2 3 2 4 2 5 2 6 2 7 2 8 2 9
Days FIG. 1. An example of population growth in continuous culture of the rotifer algal food supplied was Nannochlorisoculata at 5 × 106 cells/mL at 24°C. The population was censused daily and exposed to copper continuously. Data are from [4].
Brachionuscalyciflorus.The
The state variables are given as follows:
x(t) = CE(t) = C(t) =
biomass of the population in mg at time t, concentration of toxicant in m g / L in the environment at time t, a m o u n t of toxicant in mg in the biomass of the population at time t.
In addition, if we define
Cb(t)
=
concentration of toxicant in the biomass of the population at time t,
we have the relation
Cb={CIx
if x * O , if x = O .
(1)
142
DIANA MARIA THOMAS ET AL.
The model will be formulated by the following system of differential equations:
dX = xr( Cb ) ( 1 d---{
x) K(Cb ) ,
(2)
dCe 'dt = - aCex + b C - hC e + u,
(3)
dC dt = aCEx - IxC,
(4)
x(O) = x0 > o,
(5)
c E ( o ) : c~ > o,
(6)
c ( o ) = o.
(7)
with the initial condition
The parameters a, b, h, and /, are all positive. We shall refer to the above system as SYSI. Equation (2) is the logistic model where r(C b) represents the growth rate and K(C b) represents the carrying capacity. The functions r and K are twice continuously differentiable and have the property that
dr dCb
0
and
dK
-a~b ~-~<0.
Moreover, we shall require that K(C b) and r(C o) be bounded away from zero by some constants, ~ and 3':
0 < 6 <~K(Cb) and
O < 3" <~r( Cb ). K and r are required to be bounded away from zero for two reasons. Mathematically, we wish to have a local Lipschitz condition on the right-hand side of our system in order to guarantee local existence and uniqueness of the solutions to SYSI. Biologically, we do not wish to consider an environment that is so hostile that it has a carrying capacity near zero. The first term in (3) represents the loss of C E due to reaction with the biomass, and the second term represents the gain of C E due to metabolic losses of toxicant from the biomass. All other losses of
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
143
toxicant from the environment are described by the linear term - hC e. The exogenous rate of i n p u t / o u t p u t of toxicant into the environment is modeled by a C 1 function u = u(t, x, Ce, C). Output of toxicant from the environment would be due to a cleanup effort and would be represented by negative values of u. We shall require that u be a differentiable function. In Equation (4), aCEx is the amount of toxicant that is absorbed by the biomass from the environment, and - / x C represents all losses of C. /x can be divided into >=b+/3, where - / 3 C describes losses of the toxicant from the biomass that is not returned to the environment. Although it is not necessary for analytic purposes, we need explicit formulas for r and K for numerical simulations. For these circumstances, we define r and K by the hyperbolas
r( Cb ) = ro + rlC b Cb + 1
(8)
and
g(Cb)=
pC b -t- K o Cb+I '
(9)
where y ~< r 1 < r 0 and 6 ~< p < K 0. K 0 and r 0 are the carrying capacity and population growth rate, respectively, in the absence of the toxicant. 4.
DETERMINATION OF THE CONTROL FUNCTION
u(t, x, CE, C) Since the right-hand sides of (2)-(4) are locally Lipschitz, a unique solution to SYSI exists locally. In order to discuss further properties of SYSI, we shall need to determine which variable to control and how we wish to control it. We wish to force lim CE(t )
= M,
(10)
where M is a positive constant. M will be considered as a water quality limit determining the amount of toxicant allowed in the environment. In an ideal situation we would prefer to have the input u, to be able to drive C e to M without taking into account the state of the system. This type of system is called an open-loop system [8]. An open-loop
DIANA MARIA THOMAS ET AL.
144 Ref~ren~
Inputs+(
u
l
(SYSl)
Outout
: (x,CE,C)
!
t ~---w~"I'ev~In t Elwlrelmumt FIG. 2. Block diagram of closed-loop system governed by SYSI.
system is desirable because we will not have to measure the states and feed the information back to the controller in order to achieve (10). In a closed-loop system, however, we must feed the information of the amount of pollutant in the environment back to the controller in order for C E to reach the desired value, M. The advantage of the closed-loop system is that the process becomes self-adjusting or self-correcting [9]. We chose to model the process by the closed-loop system that is depicted in the block diagram in Figure 2. To obtain an explicit closed-form expression for u we shall employ the method of feedback linearization. Feedback linearization is a process in which a nonlinear system is algebraically transformed into a fully or partially linear system [10]. Once u is chosen to cancel the nonlinearities of the system we must shift the now linear system so that C e converges to M. Shifting the linear equation in this manner is called pole shifting [11]. Figure 3 portrays a block diagram of feedback linearization combined with pole shifting.
Reference
Inputs+.~ ~
T ~ n Sytem of [~_~ Nonnnwr
u - g(x,Cs,C)
Pale Shiffina LeaD
FIG. 3. Block diagram of feedback linearization loop combined with pole-shifting loop.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
145
We first formally make the transformation of variables,
or
11,
Then Equation (3) can be rewritten in terms of y as bC = - ax + ~
u - h + Me-----y.
(12)
Now we choose u to cancel the nonlinearities of (12) and satisfy (10): u=
(
bc
Me y ax--M-~ey + h-
hy
)
,
(13)
where A is a positive parameter. Substituting (13) into (12) yields + hy = 0
y = yo e - x t ,
where
The transformation (11) results in partially linear dynamics for SYSI; that is, Equation (3) has been transformed to a linear equation. The whole system may be transformed to a fully linear system using formulas in [11]; however, the formulas require explicit expressions for r and K. We may be able to find a simpler expression for u because we only need to cancel nonlinearities. However, the term a x C E m u s t appear in u if we are to use feedback linearization. We will show later in the paper that our dynamical results hold for any choice of u if Ce converges to M uniformly. The global existence and nonnegativity results, on the other hand, must be proved for each choice of u. Our choice of u in (13) can be written as u = - Ce - by. The first pollutant sents the that it is 4.1. This
term represents the feedback of the rate at which the level of in the environment is changing, and the second term reprepole-shifting process. If we examine the second term, we see positive if CE < M, which is the case considered in Theorem is reasonable because we wish to have this controlled input
146
DIANA MARIA THOMAS ET AL.
increase. If C e > M, the last term is negative. Thus, we need to remove toxicant to achieve our goal. For existence, uniqueness, and nonnegativity, we consider two cases. The first is the case where the initial amount of pollutant in the environment is lower than the water quality restriction, that is,
C~ < M . The second will examine the case where the initial amount of pollutant is higher than the water quality restriction and thus involves cleanup of the environment. In order to establish nonnegativity of solutions to SYSI we shall need the following theorem.
THEOREM 4.1 ( Nonnegativity of Solutions) [12] Suppose m is a positive integer, ~+ =[0,w), F = ( F / ) ~ n is a locally Lipschitz function from ~'~ into R '~ that is quasi-positive: m
~=(~i)1 e~m
and
~=O~Fk(~)~>O.
Under these circumstances, the ordinary differential equation y = F(y),
y(0) = ~" e ~m,
t >/0,
has a unique, noncontinuable, nonnegative solution, y(-, ~") on an interval [O,b¢), where O
THEOREM 4.2 Assume that C~ < M. Let u be defined as in (13). Then, (i) Unique nonnegative solutions to ( SYSI ) exist globally. (ii) The solution to (3) is given by C E =
Me y,
and hence lim CE(t )
= M.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
Proof.
147
First notice that u(O, x ( O ) , C E ( O ) , C ( O ) ) = MeYo[aXo + h - Zy0] > 0.
Thus, there is an e 1 > 0 such that u >1 0 on [0, el). So on [0, 81) , the right-hand side of SYSI is quasi-positive, and so by T h e o r e m 4.1 there is an e 2 such that the solution to SYSI is unique and nonnegative on [0, e2).
W e claim that e z = + oo. T o see this, note that by (13) the solution to (3) on [0, e2) is given by C E = M e y. Also on [0, e 2) we have the differential inequalities
(x)
x7 1-~
(
<<.Jc
and a C ~ x - tzC <<.(~ <~aMx - IzC. Therefore,
o < t(xO) <. x <~L ( x o ) and 0 < a C ~ l ( x ° ) ( 1 - e - ~ ' ) <~C <<.a M L ( x ° ) tx tz where l(xo)=min[y,6
]
and
L(xo)=max[xo,K(O)].
Thus
Ixl+ ICEI+ Ifl <~
as t -o ~2"
By T h e o r e m 4.1, e 2 = +o0. Since nonnegative solutions exist on [0, oo), we have that C e = Me y on [0, ~).
148
DIANA MARIA THOMAS ET AL.
Since C E = Me y, we are able to eliminate C e from the system and consider the reduced problem,
dX d--[ = xr( Cb ) ( 1
x )) K(Q
'
(14)
dC d--t = aMeYx - IzC,
(15)
x ( 0 ) = x 0 > 0,
(16)
C ( 0 ) = 0.
(17)
We shall refer to the reduced system as SYSII. N o w we shall look at the case where C~ > M.
THEOREM 4.3 Suppose that C~ > M. Let u be defined as in (13) and choose h so that axo + h
-
-
Y0
> A.
(lS)
Then, (i) Unique nonnegative solutions to (SYSI) exist globally. (ii) The solution to (3) is given by
C E = Me y, and hence lim CE( t )
Proof.
= M.
By condition (19) we have that
u(O) > o. Therefore, we are able to make the same a r g u m e n t as in the p r o o f of T h e o r e m 4.2 to obtain a unique, nonnegative solution on an interval [0, e). To show that e = + ~ we use the fact that
M<<.CE=MeY<~CE
for t >~ 0.
(19)
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
149
Using (19) we are able to obtain the differential inequalities
xy(l-~)~<~
K(O))
and
aMx - ~C < C < aC~ x - txC, which gives us that o < t(Xo)
x
L(xo)
and
aMl( xo ) 0 < - - ( 1/x- e - " t )
<~C <~
aC~ L ( xo ) /x
From here on the proof is similar to he proof of Theorem 4.2. 5.
T H E DYNAMICAL B E H A V I O R OF T H E BIOMASS
5.1. THE ASYMPTOTICALLY AUTONOMOUS SYSTEM Since SYSI is a reduction of a system of three differential equations to a system of two differential equations, we are able to apply theorems available for planar systems. However, we gain a complication in the fact that our reduced system is nonautonomous. Since SYSII is an asymptotically autonomous system, we are able to get around this difficulty.
DEFINITION A nonautonomous system of differential equations in R n,
x'=f(t,x),
(20)
is said to be asymptotically autonomous with limit equation
y'=g(y)
(21)
if f(t, x) ~ g(x) as t ~ 0%where the convergence is uniform on compact sets.
150
DIANA MARIA THOMAS ET AL.
The limiting equations for SYSII are given by
dr=xr(Cb)( 1 dt
x) K ( Cb) '
(22)
dC - aMx - tzC, dt x(O) = x 0 > 0,
(23) (24)
c ( o ) = o.
(25)
We shall refer to this system of as SYSIII. Naturally, we would hope that the solutions to SYSIII and SYSII have the same asymptotic behavior. However, it is not always true that the solution to an asymptotically autonomous system will have a structure similar to that of the solution to its limiting system as seen in examples from [13]. Thus, are require more conditions on the two systems to guarantee that their long time behavior will be the same. These conditions are provided in the Dulac criterion (see [13], e.g.).
DULAC CRITERION Let X, Y, and D be subsets of ~z, D open and simply connected, with the following properties: (i) Every bounded forward orbit of (21) in X has its ~-limit set in Y. (ii) All possible periodic orbits of (22) in Y and the closures of all possible orbits of (22) that chain equilibria of (22) are contained in D. (iii) g is continuously differentiable on D, and there is a real-valued continuously differentiable function p on D such that the divergence of Pg,
VX(pg)(xl,x2) =-~a(Pg)(xl,x2)+
(Pg)(xl,x2),
is either strictly positive almost everywhere on D or strictly negative almost everywhere on D. Under these conditions we have the following corollary from [13].
COROLLARY Let the g equilibria in X be isolated in ~2, and assume that the Dulac criterion holds. Then every bounded forward g orbit in X and every bounded forward f orbit in X converges toward a g equilibrium. We shall see in Section 5.2 that the g equilibria of our system are indeed isolated. Thus to apply the above corollary, we need only show that the Dulac criterion holds.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
151
LEMMA 5.1 Let X = Y = (0, K(0)] ×(0, aMK(O)/tz], and let D = int(X). Assume that ar'( Cb) + ¢]K'( Cb) > - / z ,
(26)
where = S(xo),
r(O)K(O)
~2
and
S(xo)
aMK(O) ml(xo) "
Then the Dulac criterion holds.
Proof.
By the differential inequalities
xy(l-~)<~Yc<~xr(O)(1
x
K(O))
and
aC~ x - tzC < C < aMx - tzC, which we obtained in the proof of T h e o r e m 4.2, we have that
0
aC~ 1( x o) /x
( 1 - - e -'~t) ~
aML( x o)
But if x o < K(O) we have that L(x o) = K(O), so 0 < x < K(O) and O
aMK(O) /z
/.~
DIANA MARIA THOMAS ET AL.
152
for t > 0 and 0 ~
x C Br ( C b)
K(C )
Now we make the change of variables, x = e w,
which produces the system
=r(Cb)
( eW) 1
K(Co )
Cb = a M - ~ C b - Cbr( Cb ) +
=g,(w,Co),
eWCbr(Cb) K(Cb)
g2(w,Q).
Setting p = 1 we have that
V ' g ( w , C D = Ogl t w , C bJ~ + -U~b(w,Cb) Og2 Ow ~ -
eWr(C b) K(Cb )
tz - Cbrb(Cb) -- r ( C b )
+ eW[ C a r ' ( C b ) K ( C 6 )
[ <~ - tz - S ( x o ) r ' ( C b ) -
_ r(Ca)K'(Cb)
K (Q) K(O) r ( O ) K ' ( C b ) 62
+
r(Ca) ]
(27)
By condition (26), expression (27) is less than zero. Thus g satisfies (iii), and this completes the proof. • Condition (26) is a restriction on the rate at which r and K may change. Since we wish to preserve the dynamics of the biomass in the absence of toxicant, (26) must be met in order to guarantee avoidance of periodic behavior. We discussed the nonuniqueness of the control function at the beginning of this article. We now note that the asymptotic results will hold for any control function u for which C e converges to M uniformly.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
153
This is due to the fact that the only condition we need to satisfy to obtain our results is that the system be asymptotically autonomous. 5.2. STABILITY ANALYSIS OF THE LIMITING SYSTEM SYSIII has two equilibria, (0,0) and ( K ( a M / I X ) , ( a M / I X ) K ( a M / I x ) ) . The steady state x = 0 and C = 0 represents the extinction of the population. The following theorem provides us with the conditions under which (0, 0) would be an asymptotically stable equilibrium point of SYSIII. THEOREM 5.4 The equilibrium point (0, O) o f S Y S I I is asymptotically stable if (i) Ix - r ( 0 ) > 0
and (ii) - Ixr(0)- aMr'(O) > O. Proof
The Jacobian matrix of SYSIII at (0,0) is given by
J(O,O)
=Jr(O) r'(O)] [aM
- IX
(28)
'
with characteristic equation }~2 _+_
[ /.Z - - r(0)] A - Ixr(0) - aMr'(O).
(29)
The Hurwitz matrix of the polynomial (29) is H(0,0) = [IX-;(0)
0 - Ixr(O) - aMr'(O)
],
(30)
and applying the Hurwitz criterion (see [8]) we have that (29) has roots with negative real part if and only if Ix - r(0) > 0
(31)
[ - Ixr(O) - aMr'(O)] [ Ix - r(O)] > o.
(32)
and
This completes the proof.
•
Now, for (i) and (ii) to hold, r(0) must be small a n d / o r the rate of loss of toxicant from the biomass must be large and in addition a, M
154
DIANA MARIA THOMAS ET AL.
a n d / o r Ir'(0)l must be large. Biologically, this says that if a slowly growing population is exposed to toxicant stress, and the rate of toxicant uptake by the biomass is fast, and the initial decrease of the growth rate is fast, then the water quality limit should not be set too high because the result will be the extinction of the population. Assuming that (i) holds, condition (ii) provides us with a way to choose M in order to avoid extinction of the population. However, if M has to be very small, toxicant release might have to be prohibited to preserve the population. We have no control over condition (i). If the population is close to extinction prior to the introduction of the toxicant, extinction is likely to be the result of toxicant stress. Now, we wish to examine the case where the population does not go extinct. In this situation, the biomass tends to a new, but lower, carrying capacity. The conditions under which this steady state will be asymptotically stable are provided in the following theorem.
THEOREM5.5 The steady state
is an asymptotically stable equilibrium point of SYSIII if # + 0 > 0,
(33)
where 0 Proof
r(aM/ tx) K(aM/tz)
a___~M ) + aM K,( aM ] ]" [K(
/z
k /x l]
For notational ease let
~.=K(AM]
(a__~_M),
The Jacobian matrix of SYSIII at
---
K---7 K,(a__~_M)
aM
(K, BK) is
l -° aM
K
- I~
,
(34)
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
155
with characteristic equation A2 + ( / ~ + 0 ) A + / z 0 -aMr-r-r-~ K ' . K
(35)
The Hurwitz matrix for (35) is
H(K., -B--K) =
1
txO - a M ~ ~ K
'
(36)
and once again applying the Hurwitz criterion gives us that the roots of (35) have negative real part if and only if /z + 0 > 0
(37)
and ( tx + O)[ ~ O - a M K - K r ]
>O.
(38)
We claim that condition (37) implies condition (38). To see this, assume that (37) holds. Then we need only show that txO - a M r K ' > O. K
Now, txO - a M L K ' = O - B r-r-r-r-~K ' K K
(39)
=/x~ > O. Thus, txO - a M r K is a positive quantity, and this completes the proof.
•
156
DIANA MARIA THOMAS ET AL.
The condition ~ + O > 0 can e written as
+ ~ > - B-- K'. K
(40)
Inequality (40) will tend to hold if the magnitude of the rate of decrease of the carrying capacity at C b = B is not too large, the carrying capacity at Cb = B is not too small, the water quality limit M is not set too high, a n d / o r the rate of uptake of toxicant from the environment, a, is not too fast. Thus, we may choose M so that x tends to R as t ~ ~. In fact, we can even control the value of the new carrying capacity, ~', by choosing M appropriately. The next section displays some numerical results of SYSI, SYSII, and SYSIII. 6.
N U M E R I C A L RESULTS
Formulas (8) and (9) are used in all of the numerical simulations. Figures 4 - 7 were generated by the package Differential Systems in which a fifth-order Runge-Kutta method was used. The parameter values corresponding to each figure are presented in Table 1. Figure 4 shows several trajectories of the biomass and the amount of toxicant in the biomass for SYSIII. We are able to see that the trajectories of the system are converging to the asymptotically stable equilibrium point (K, BK). Figure 5 is a phase portrait of SYSIII that also shows orbits converging to the equilibrium value. Figure 6 displays trajectories of the biomass and the total amount of toxicant in the biomass for the nonautonomous system, SYSII. As our analysis from Section 5 indicated, the trajectories of SYSII are converging to the stable equilibrium point of the limiting system SYSIII. Further trajectories for SYSII of the biomass are shown in Figure 7 for varying parameter values. Here again we numerically see the asymptotic results of Section 5. Figure 8 is a numerical simulation of several trajectories of SYSI. The solid line represents biomass, the dashed line represents the concentration of toxicant in the environment, and the dash-dot line represents the amount of toxicant in the biomass. Figure 9 is a graph of the control function versus time. Since x, Ce, and C all converge, we have that u becomes relatively constant for large time.
A C O N T R O L P R O B L E M IN A P O L L U T E D E N V I R O N M E N T
t.r
200
490
157
590
800
1000
500
80g
100|
]
(a) 200
460
-i
i"
~T
.~ 21
(b) FIG. 4. Several trajectories of the biomass and the amount of toxicant in the biomass for SYSIII.
158
DIANA MARIA THOMAS ET AL. 9
I
Z
3
4
5
5
7
8
g
Ill
'¢1
!
"
'
,
BLomasa
FIG. 5. Phase portrait of SYSIII.
Figure 10 is a plot of the biomass versus time in the absence of toxicant modeled by logistic growth (dashed line) and a plot of the biomass versus time in the presence of toxicant modeled by SYSI. We are able to see that the effect of the toxicant on the biomass is a lowered carrying capacity. This inhibited carrying capacity effect also occurred in experimental data as seen in Figure 1. 7.
CONCLUSIONS
From the theoretical analysis of the model it was found that if the toxicant input into the environment is appropriately controlled, the biomass will tend to a new lowered carrying capacity. The value of this new carrying capacity can be determined from the water quality criteria. It was also shown that in cases where the population was stressed before the toxicant was introduced, prohibition of the toxicant may be the only solution to preserve the population. The analytical results are supported by experimental data. Lowered carrying capacity was seen in rotifer populations that were exposed to low concentrations of copper. There is much to be done for future work. A more realistic control problem should require that C E ~< M. There are extensive results on
A C O N T R O L PROBLEM IN A POLLUTED E N V I R O N M E N T 200
400
159
6OO
800
I BOll
600
800
lOON
qr
mu,)
i
Time
(a) 200
400
t
Time
(b) FIG. 6. Several trajectories of the biomass and the amount of toxicant in the biomass for SYSII.
160
D I A N A M A R I A T H O M A S ET AL 500
800
ION
eee
iooe
Ttme
(a)
(b) FIG. 7. Further trajectories of the biomass with different parameter values.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
161
TABLE 1 Parameter Values Corresponding to Figures 4-10 Parameter Figure 4 5 6 7a 7b 8 9 10
ro
rI
p Ko
0.1 0.1 0.1 0.01 0.01 0.1 0.1 0.1
0.01 0.01 0.01 0.1 0.3 0.01 0.01 0.01
2 2 2 3 3 2 2 2
5 5 5 6 4 5 5 5
a
M
/x
}t
Yo
b
h
g
BK
0.1 0.1 0.1 0.01 0.001 0.01 0.1 0.1
0.5 0.5 0.5 0.4 0.5 0.5 0.5 0.5
0.3 0.3 0.3 0.2 0.2 0.2 0.3 0.3
0.01 0.01 0.01 0.1 0.1 0.01 0.01 0.01
-3.1 -3.1 -3.1 -0.29 -0.51 NA -3.1 -3.1
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
4.57 4.57 4.57 5.94 4.008 4.93 4.57 4.57
0.762 0.762 0.762 0.12 0.01 0.123 0.762 0.762
these types of c o n s t r a i n t s in control theory literature. Also, p e r h a p s a n a g e - d e p e n d e n t m o d e l or delay m o d e l should be i n c o r p o r a t e d to a c c o u n t for the oscillation in the rotifer b i o m a s s growth that is seen in Figure 1. Again, there are m a n y resources available for this type of p r o b l e m . A n adaptive p a r a m e t e r c o n t r o l p r o b l e m should be investigated b e c a u s e the p a r a m e t e r s are n e v e r m e a s u r e d exactly. T h e r e are m a n y t e c h n i q u e s available to adjust for errors in m e a s u r e m e n t s in the r e a l m of control
r
\ .
.
.
.
.
.
200
.
.
.
.
.
.
400
.
.
.
.
.
.
600
.
.
.
.
m "
800
| I000
FIG. 8. Tr~ectories of SYSI. Solid curves are x; dashed curves are dashed-dot curves are C.
time
CE; and
162
DIANA MARIA THOMAS ET AL.
0,12
0.1
0.08
0,06
0.04
O. 02
100
200
300
400
500
600
700
FiG. 9. Plot of the control function u versus time.
x 5
3
2
1'
lOO
200
300
400
soo
~oo
700
FIG. 10. Plots of biomass behavior in the absence of toxicant (dashed curve) and in the presence of toxicant (solid curve) versus time.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
163
theory. Finally, o n e should c o n s i d e r the o p t i m a l c o n t r o l p r o b l e m , w h e r e we m i n i m i z e the cost associated with controlling the p o l l u t i o n in the environment.
We thank Almut Burchard and Hugo Leiva for the time spent discussing the material and for their suggestions. We also thank Dr. H. I. Freedman for his insight on modifications for our model. REFERENCES 1 R.L. Wallace and T. W. Snell, Rotifera, in Ecology and Systematics of North American FreshwaterInvertebrates, J. H. Thorp and A. P. Covich Eds., Academic, New York, 1991, pp. 187-248. 2 T. W. Snell and B. D. Moffat, A two day life cycle test with the rotifer Brachionus calyciflorus, Environ. Toxicol. Chem., 11:1249-1257 (1992). 3 C.R. Janssen, G. Personne, and T. W. Snell, Ecotoxicological studies with the freshwater rotifer Brachionus calyciflorus. VI. Short chronic toxicity tests with rotifers, Aquat. Toxicol., in press (1994). 4 T.W. Snell and C. Janssen, Rotifers in ecotoxicology: a review, Hydrobiologia, in press (1994). 5 C. Janssen, The use of sublethal criteria for toxicity test with the freshwater rotifer Brachionus ca~ciflorus (Pallas), Ph.D. Thesis, Univ. Ghent, Belgium, 1992. 6 T.G. Hallam, C. E. Clark, and G. S. Jordan, Effects of toxicants on populations: a qualitative approach. II. First order kinetics, J. Math. Biol., 18:25-37 (1983). 7 T.G. Hallam and C. E. Clark, Effects of toxicants on populations: a qualitative approach. I. Equilibrium environmental exposure, Ecol. Modeling, 18:291-304 (1983). 8 S. Barnett and R. G. Cameron, Introduction to Mathematical Control Theory, Oxford Univ. Press, New York, 1985, pp. 351-379. 9 E.B. Lee and L. Markus, Foundations of Optimal Control Theory, Wiley, New Yori, 1967, p. 24. 10 J.-J. E. Slotine and W. Li, Applied Nonlinear Control, Prentice-Hall, Englewood Cliffs, NJ, 1967, p. 207. 11 H.K. Khalil, Nonlinear Systems, Macmillan, New York, 1992. 12 R.H. Martin and M. Pierre, Nonlinear Reaction Diffusion Systems--Nonlinear Equations in the Applied Sciences, Academic, New York, 1992. 13 H.R. Thieme, Asymptotically autonomous differential equations in the plane, Rocky Mt. J. Math., 21:351-379 (1994).
A Control Problem in a Polluted Environment DIANA MARIA THOMAS School of Mathematics, Georgia Tech, Atlanta, Georgia 30332 TERRY WAYNE SNELL School of Biology, Georgia Tech, Atlanta, Georgia 30332 AND
SAYED MURTAZA JAFFAR School of Mechanical Engineering, Georgia Tech, Atlanta, Georgia 30332
Received 29 August 1994; revised 16 June 1995
ABSTRACT Due to the importance of remediating polluted environments, many mathematical models have been developed to describe population-toxicant interactions. This article considers the situation where the amount of pollution in the environment is restricted to lie below a certain value M. A model is presented and viewed as a control problem in which the toxicant input into the environment is the controllable variable. An explicit formula is computed for the input function that meets certain restrictions. By analytically studying the dynamics of the population, it is possible to provide conditions on how to set M in order to guarantee survival of the population without a significant reduction of its carrying capacity. Numerical and experimental data are described that support these results.
1.
INTRODUCTION
Efforts to clean up polluted ecosystems have led to many regulations on amounts of toxicants allowed in the environment. This includes water quality criteria set by the Environmental Protection Agency (EPA) to limit certain chemicals in freshwater environments. Establishment of regulations depend on the type of behavior desired in the system. Generally, the objective is to prevent species loss and maintain ecosystem function. The aim of this paper is to study the dynamics of a single-species population in a polluted freshwater environment where the input of pollution is controlled. We focus on a specific example using rotifers, an important component of the zooplankton in many aquatic communities.
MATHEMATICAL BIOSC1ENCES 133:139-163 (1996) 0025-5564/96/$15.00 © Elsevier Science Inc., 1996 SSDI 0025-5564(95)00091-Q 655 Avenue of the Americas, New York, NY 10010
140
DIANA MARIA THOMAS ET AL.
Section 2 describes biological aspects of rotifers and presents some experimental results. In Section 3 a model is presented that describes population-toxicant interactions. The control problem is posed and solved in Section 4, and several properties of the model are consequently derived. In Section 5 we discuss the dynamical aspects of the model by conducting a stability analysis. Conditions are consequently derived that determine how to set the restrictions on the amount of toxicant in the environment in order to guarantee survival of the population without a significant reduction of its carrying capacity. In addition, some numerical results are presented. Section 6 summarizes the results of the work and poses future problems. 2. THE BIOLOGICAL EXAMPLE AND METHODS Rotifers are multicellular aquatic animals that eat single-celled algae in lakes, ponds, rivers, and wetlands [1]. Although they are small (0.1-1 mm in length), they can occur in high densities and modify the composition of phytoplankton populations by selective grazing. Most rotifers are herbivores providing a critical link in aquatic food webs between microscopic plants and their fish predators. The dynamics of laboratory and field populations of rotifers have been characterized in many ecological studies [1]. Rotifer population dynamics also have been used to examine the effects of pollution [2-4], with population growth are, r, the usual criterion. Janssen et al. illustrated how carrying capacity, K, could be a useful criterion in continuous cultures of the rotifer Brachionus calyciflorus [5]. They showed that 30-day exposures (approximately 15 generations) to 10 /zg copper/L reduced K to about half of the control, which is the unexposed culture (Figure 1). The copper concentration causing significant reduction in K is about half that of LCs0 [5]. Therefore, the effect of pollution on rotifer populations is to lower both r and K. The objective of environmental regulation is to set limits on how much pollution can be allowed into the environment without significantly reducing indigenous populations. The model that we have developed provides guidance about the maximum amount of pollution allowable before substantially lowering the carrying capacity. Using a system of differential equations based on the logistic population growth model, we derive a relationship between environmental concentrations of pollution and population dynamics. 3. THE MODEL AND PRELIMINARIES In this section, we present a modification of the model suggested in [6] that displays logistic population growth coupled with toxicant stress.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
12o--]
141
Brachionus calyciflorus _
c o p p e r L C 5 0 = 2 6 lag/L
controJ~
100t 807
I
2.5 pg Cu
~
_ _
o.ns,,, .I v
v/--I
a
I
l
I
l
l
I
I
-
i
i
l
J
l
l
l
J
l
l
l
l
l
i
J
l
l
l
l
l
I 2 3 4 5 6 7 8 9 1011 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0 2 1 2 2 2 3 2 4 2 5 2 6 2 7 2 8 2 9
Days FIG. 1. An example of population growth in continuous culture of the rotifer algal food supplied was Nannochlorisoculata at 5 × 106 cells/mL at 24°C. The population was censused daily and exposed to copper continuously. Data are from [4].
Brachionuscalyciflorus.The
The state variables are given as follows:
x(t) = CE(t) = C(t) =
biomass of the population in mg at time t, concentration of toxicant in m g / L in the environment at time t, a m o u n t of toxicant in mg in the biomass of the population at time t.
In addition, if we define
Cb(t)
=
concentration of toxicant in the biomass of the population at time t,
we have the relation
Cb={CIx
if x * O , if x = O .
(1)
142
DIANA MARIA THOMAS ET AL.
The model will be formulated by the following system of differential equations:
dX = xr( Cb ) ( 1 d---{
x) K(Cb ) ,
(2)
dCe 'dt = - aCex + b C - hC e + u,
(3)
dC dt = aCEx - IxC,
(4)
x(O) = x0 > o,
(5)
c E ( o ) : c~ > o,
(6)
c ( o ) = o.
(7)
with the initial condition
The parameters a, b, h, and /, are all positive. We shall refer to the above system as SYSI. Equation (2) is the logistic model where r(C b) represents the growth rate and K(C b) represents the carrying capacity. The functions r and K are twice continuously differentiable and have the property that
dr dCb
0
and
dK
-a~b ~-~<0.
Moreover, we shall require that K(C b) and r(C o) be bounded away from zero by some constants, ~ and 3':
0 < 6 <~K(Cb) and
O < 3" <~r( Cb ). K and r are required to be bounded away from zero for two reasons. Mathematically, we wish to have a local Lipschitz condition on the right-hand side of our system in order to guarantee local existence and uniqueness of the solutions to SYSI. Biologically, we do not wish to consider an environment that is so hostile that it has a carrying capacity near zero. The first term in (3) represents the loss of C E due to reaction with the biomass, and the second term represents the gain of C E due to metabolic losses of toxicant from the biomass. All other losses of
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
143
toxicant from the environment are described by the linear term - hC e. The exogenous rate of i n p u t / o u t p u t of toxicant into the environment is modeled by a C 1 function u = u(t, x, Ce, C). Output of toxicant from the environment would be due to a cleanup effort and would be represented by negative values of u. We shall require that u be a differentiable function. In Equation (4), aCEx is the amount of toxicant that is absorbed by the biomass from the environment, and - / x C represents all losses of C. /x can be divided into >=b+/3, where - / 3 C describes losses of the toxicant from the biomass that is not returned to the environment. Although it is not necessary for analytic purposes, we need explicit formulas for r and K for numerical simulations. For these circumstances, we define r and K by the hyperbolas
r( Cb ) = ro + rlC b Cb + 1
(8)
and
g(Cb)=
pC b -t- K o Cb+I '
(9)
where y ~< r 1 < r 0 and 6 ~< p < K 0. K 0 and r 0 are the carrying capacity and population growth rate, respectively, in the absence of the toxicant. 4.
DETERMINATION OF THE CONTROL FUNCTION
u(t, x, CE, C) Since the right-hand sides of (2)-(4) are locally Lipschitz, a unique solution to SYSI exists locally. In order to discuss further properties of SYSI, we shall need to determine which variable to control and how we wish to control it. We wish to force lim CE(t )
= M,
(10)
where M is a positive constant. M will be considered as a water quality limit determining the amount of toxicant allowed in the environment. In an ideal situation we would prefer to have the input u, to be able to drive C e to M without taking into account the state of the system. This type of system is called an open-loop system [8]. An open-loop
DIANA MARIA THOMAS ET AL.
144 Ref~ren~
Inputs+(
u
l
(SYSl)
Outout
: (x,CE,C)
!
t ~---w~"I'ev~In t Elwlrelmumt FIG. 2. Block diagram of closed-loop system governed by SYSI.
system is desirable because we will not have to measure the states and feed the information back to the controller in order to achieve (10). In a closed-loop system, however, we must feed the information of the amount of pollutant in the environment back to the controller in order for C E to reach the desired value, M. The advantage of the closed-loop system is that the process becomes self-adjusting or self-correcting [9]. We chose to model the process by the closed-loop system that is depicted in the block diagram in Figure 2. To obtain an explicit closed-form expression for u we shall employ the method of feedback linearization. Feedback linearization is a process in which a nonlinear system is algebraically transformed into a fully or partially linear system [10]. Once u is chosen to cancel the nonlinearities of the system we must shift the now linear system so that C e converges to M. Shifting the linear equation in this manner is called pole shifting [11]. Figure 3 portrays a block diagram of feedback linearization combined with pole shifting.
Reference
Inputs+.~ ~
T ~ n Sytem of [~_~ Nonnnwr
u - g(x,Cs,C)
Pale Shiffina LeaD
FIG. 3. Block diagram of feedback linearization loop combined with pole-shifting loop.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
145
We first formally make the transformation of variables,
or
11,
Then Equation (3) can be rewritten in terms of y as bC = - ax + ~
u - h + Me-----y.
(12)
Now we choose u to cancel the nonlinearities of (12) and satisfy (10): u=
(
bc
Me y ax--M-~ey + h-
hy
)
,
(13)
where A is a positive parameter. Substituting (13) into (12) yields + hy = 0
y = yo e - x t ,
where
The transformation (11) results in partially linear dynamics for SYSI; that is, Equation (3) has been transformed to a linear equation. The whole system may be transformed to a fully linear system using formulas in [11]; however, the formulas require explicit expressions for r and K. We may be able to find a simpler expression for u because we only need to cancel nonlinearities. However, the term a x C E m u s t appear in u if we are to use feedback linearization. We will show later in the paper that our dynamical results hold for any choice of u if Ce converges to M uniformly. The global existence and nonnegativity results, on the other hand, must be proved for each choice of u. Our choice of u in (13) can be written as u = - Ce - by. The first pollutant sents the that it is 4.1. This
term represents the feedback of the rate at which the level of in the environment is changing, and the second term reprepole-shifting process. If we examine the second term, we see positive if CE < M, which is the case considered in Theorem is reasonable because we wish to have this controlled input
146
DIANA MARIA THOMAS ET AL.
increase. If C e > M, the last term is negative. Thus, we need to remove toxicant to achieve our goal. For existence, uniqueness, and nonnegativity, we consider two cases. The first is the case where the initial amount of pollutant in the environment is lower than the water quality restriction, that is,
C~ < M . The second will examine the case where the initial amount of pollutant is higher than the water quality restriction and thus involves cleanup of the environment. In order to establish nonnegativity of solutions to SYSI we shall need the following theorem.
THEOREM 4.1 ( Nonnegativity of Solutions) [12] Suppose m is a positive integer, ~+ =[0,w), F = ( F / ) ~ n is a locally Lipschitz function from ~'~ into R '~ that is quasi-positive: m
~=(~i)1 e~m
and
~=O~Fk(~)~>O.
Under these circumstances, the ordinary differential equation y = F(y),
y(0) = ~" e ~m,
t >/0,
has a unique, noncontinuable, nonnegative solution, y(-, ~") on an interval [O,b¢), where O
THEOREM 4.2 Assume that C~ < M. Let u be defined as in (13). Then, (i) Unique nonnegative solutions to ( SYSI ) exist globally. (ii) The solution to (3) is given by C E =
Me y,
and hence lim CE(t )
= M.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
Proof.
147
First notice that u(O, x ( O ) , C E ( O ) , C ( O ) ) = MeYo[aXo + h - Zy0] > 0.
Thus, there is an e 1 > 0 such that u >1 0 on [0, el). So on [0, 81) , the right-hand side of SYSI is quasi-positive, and so by T h e o r e m 4.1 there is an e 2 such that the solution to SYSI is unique and nonnegative on [0, e2).
W e claim that e z = + oo. T o see this, note that by (13) the solution to (3) on [0, e2) is given by C E = M e y. Also on [0, e 2) we have the differential inequalities
(x)
x7 1-~
(
<<.Jc
and a C ~ x - tzC <<.(~ <~aMx - IzC. Therefore,
o < t(xO) <. x <~L ( x o ) and 0 < a C ~ l ( x ° ) ( 1 - e - ~ ' ) <~C <<.a M L ( x ° ) tx tz where l(xo)=min[y,6
]
and
L(xo)=max[xo,K(O)].
Thus
Ixl+ ICEI+ Ifl <~
as t -o ~2"
By T h e o r e m 4.1, e 2 = +o0. Since nonnegative solutions exist on [0, oo), we have that C e = Me y on [0, ~).
148
DIANA MARIA THOMAS ET AL.
Since C E = Me y, we are able to eliminate C e from the system and consider the reduced problem,
dX d--[ = xr( Cb ) ( 1
x )) K(Q
'
(14)
dC d--t = aMeYx - IzC,
(15)
x ( 0 ) = x 0 > 0,
(16)
C ( 0 ) = 0.
(17)
We shall refer to the reduced system as SYSII. N o w we shall look at the case where C~ > M.
THEOREM 4.3 Suppose that C~ > M. Let u be defined as in (13) and choose h so that axo + h
-
-
Y0
> A.
(lS)
Then, (i) Unique nonnegative solutions to (SYSI) exist globally. (ii) The solution to (3) is given by
C E = Me y, and hence lim CE( t )
Proof.
= M.
By condition (19) we have that
u(O) > o. Therefore, we are able to make the same a r g u m e n t as in the p r o o f of T h e o r e m 4.2 to obtain a unique, nonnegative solution on an interval [0, e). To show that e = + ~ we use the fact that
M<<.CE=MeY<~CE
for t >~ 0.
(19)
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
149
Using (19) we are able to obtain the differential inequalities
xy(l-~)~<~
K(O))
and
aMx - ~C < C < aC~ x - txC, which gives us that o < t(Xo)
x
L(xo)
and
aMl( xo ) 0 < - - ( 1/x- e - " t )
<~C <~
aC~ L ( xo ) /x
From here on the proof is similar to he proof of Theorem 4.2. 5.
T H E DYNAMICAL B E H A V I O R OF T H E BIOMASS
5.1. THE ASYMPTOTICALLY AUTONOMOUS SYSTEM Since SYSI is a reduction of a system of three differential equations to a system of two differential equations, we are able to apply theorems available for planar systems. However, we gain a complication in the fact that our reduced system is nonautonomous. Since SYSII is an asymptotically autonomous system, we are able to get around this difficulty.
DEFINITION A nonautonomous system of differential equations in R n,
x'=f(t,x),
(20)
is said to be asymptotically autonomous with limit equation
y'=g(y)
(21)
if f(t, x) ~ g(x) as t ~ 0%where the convergence is uniform on compact sets.
150
DIANA MARIA THOMAS ET AL.
The limiting equations for SYSII are given by
dr=xr(Cb)( 1 dt
x) K ( Cb) '
(22)
dC - aMx - tzC, dt x(O) = x 0 > 0,
(23) (24)
c ( o ) = o.
(25)
We shall refer to this system of as SYSIII. Naturally, we would hope that the solutions to SYSIII and SYSII have the same asymptotic behavior. However, it is not always true that the solution to an asymptotically autonomous system will have a structure similar to that of the solution to its limiting system as seen in examples from [13]. Thus, are require more conditions on the two systems to guarantee that their long time behavior will be the same. These conditions are provided in the Dulac criterion (see [13], e.g.).
DULAC CRITERION Let X, Y, and D be subsets of ~z, D open and simply connected, with the following properties: (i) Every bounded forward orbit of (21) in X has its ~-limit set in Y. (ii) All possible periodic orbits of (22) in Y and the closures of all possible orbits of (22) that chain equilibria of (22) are contained in D. (iii) g is continuously differentiable on D, and there is a real-valued continuously differentiable function p on D such that the divergence of Pg,
VX(pg)(xl,x2) =-~a(Pg)(xl,x2)+
(Pg)(xl,x2),
is either strictly positive almost everywhere on D or strictly negative almost everywhere on D. Under these conditions we have the following corollary from [13].
COROLLARY Let the g equilibria in X be isolated in ~2, and assume that the Dulac criterion holds. Then every bounded forward g orbit in X and every bounded forward f orbit in X converges toward a g equilibrium. We shall see in Section 5.2 that the g equilibria of our system are indeed isolated. Thus to apply the above corollary, we need only show that the Dulac criterion holds.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
151
LEMMA 5.1 Let X = Y = (0, K(0)] ×(0, aMK(O)/tz], and let D = int(X). Assume that ar'( Cb) + ¢]K'( Cb) > - / z ,
(26)
where = S(xo),
r(O)K(O)
~2
and
S(xo)
aMK(O) ml(xo) "
Then the Dulac criterion holds.
Proof.
By the differential inequalities
xy(l-~)<~Yc<~xr(O)(1
x
K(O))
and
aC~ x - tzC < C < aMx - tzC, which we obtained in the proof of T h e o r e m 4.2, we have that
0
aC~ 1( x o) /x
( 1 - - e -'~t) ~
aML( x o)
But if x o < K(O) we have that L(x o) = K(O), so 0 < x < K(O) and O
aMK(O) /z
/.~
DIANA MARIA THOMAS ET AL.
152
for t > 0 and 0 ~
x C Br ( C b)
K(C )
Now we make the change of variables, x = e w,
which produces the system
=r(Cb)
( eW) 1
K(Co )
Cb = a M - ~ C b - Cbr( Cb ) +
=g,(w,Co),
eWCbr(Cb) K(Cb)
g2(w,Q).
Setting p = 1 we have that
V ' g ( w , C D = Ogl t w , C bJ~ + -U~b(w,Cb) Og2 Ow ~ -
eWr(C b) K(Cb )
tz - Cbrb(Cb) -- r ( C b )
+ eW[ C a r ' ( C b ) K ( C 6 )
[ <~ - tz - S ( x o ) r ' ( C b ) -
_ r(Ca)K'(Cb)
K (Q) K(O) r ( O ) K ' ( C b ) 62
+
r(Ca) ]
(27)
By condition (26), expression (27) is less than zero. Thus g satisfies (iii), and this completes the proof. • Condition (26) is a restriction on the rate at which r and K may change. Since we wish to preserve the dynamics of the biomass in the absence of toxicant, (26) must be met in order to guarantee avoidance of periodic behavior. We discussed the nonuniqueness of the control function at the beginning of this article. We now note that the asymptotic results will hold for any control function u for which C e converges to M uniformly.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
153
This is due to the fact that the only condition we need to satisfy to obtain our results is that the system be asymptotically autonomous. 5.2. STABILITY ANALYSIS OF THE LIMITING SYSTEM SYSIII has two equilibria, (0,0) and ( K ( a M / I X ) , ( a M / I X ) K ( a M / I x ) ) . The steady state x = 0 and C = 0 represents the extinction of the population. The following theorem provides us with the conditions under which (0, 0) would be an asymptotically stable equilibrium point of SYSIII. THEOREM 5.4 The equilibrium point (0, O) o f S Y S I I is asymptotically stable if (i) Ix - r ( 0 ) > 0
and (ii) - Ixr(0)- aMr'(O) > O. Proof
The Jacobian matrix of SYSIII at (0,0) is given by
J(O,O)
=Jr(O) r'(O)] [aM
- IX
(28)
'
with characteristic equation }~2 _+_
[ /.Z - - r(0)] A - Ixr(0) - aMr'(O).
(29)
The Hurwitz matrix of the polynomial (29) is H(0,0) = [IX-;(0)
0 - Ixr(O) - aMr'(O)
],
(30)
and applying the Hurwitz criterion (see [8]) we have that (29) has roots with negative real part if and only if Ix - r(0) > 0
(31)
[ - Ixr(O) - aMr'(O)] [ Ix - r(O)] > o.
(32)
and
This completes the proof.
•
Now, for (i) and (ii) to hold, r(0) must be small a n d / o r the rate of loss of toxicant from the biomass must be large and in addition a, M
154
DIANA MARIA THOMAS ET AL.
a n d / o r Ir'(0)l must be large. Biologically, this says that if a slowly growing population is exposed to toxicant stress, and the rate of toxicant uptake by the biomass is fast, and the initial decrease of the growth rate is fast, then the water quality limit should not be set too high because the result will be the extinction of the population. Assuming that (i) holds, condition (ii) provides us with a way to choose M in order to avoid extinction of the population. However, if M has to be very small, toxicant release might have to be prohibited to preserve the population. We have no control over condition (i). If the population is close to extinction prior to the introduction of the toxicant, extinction is likely to be the result of toxicant stress. Now, we wish to examine the case where the population does not go extinct. In this situation, the biomass tends to a new, but lower, carrying capacity. The conditions under which this steady state will be asymptotically stable are provided in the following theorem.
THEOREM5.5 The steady state
is an asymptotically stable equilibrium point of SYSIII if # + 0 > 0,
(33)
where 0 Proof
r(aM/ tx) K(aM/tz)
a___~M ) + aM K,( aM ] ]" [K(
/z
k /x l]
For notational ease let
~.=K(AM]
(a__~_M),
The Jacobian matrix of SYSIII at
---
K---7 K,(a__~_M)
aM
(K, BK) is
l -° aM
K
- I~
,
(34)
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
155
with characteristic equation A2 + ( / ~ + 0 ) A + / z 0 -aMr-r-r-~ K ' . K
(35)
The Hurwitz matrix for (35) is
H(K., -B--K) =
1
txO - a M ~ ~ K
'
(36)
and once again applying the Hurwitz criterion gives us that the roots of (35) have negative real part if and only if /z + 0 > 0
(37)
and ( tx + O)[ ~ O - a M K - K r ]
>O.
(38)
We claim that condition (37) implies condition (38). To see this, assume that (37) holds. Then we need only show that txO - a M r K ' > O. K
Now, txO - a M L K ' = O - B r-r-r-r-~K ' K K
(39)
=/x~ > O. Thus, txO - a M r K is a positive quantity, and this completes the proof.
•
156
DIANA MARIA THOMAS ET AL.
The condition ~ + O > 0 can e written as
+ ~ > - B-- K'. K
(40)
Inequality (40) will tend to hold if the magnitude of the rate of decrease of the carrying capacity at C b = B is not too large, the carrying capacity at Cb = B is not too small, the water quality limit M is not set too high, a n d / o r the rate of uptake of toxicant from the environment, a, is not too fast. Thus, we may choose M so that x tends to R as t ~ ~. In fact, we can even control the value of the new carrying capacity, ~', by choosing M appropriately. The next section displays some numerical results of SYSI, SYSII, and SYSIII. 6.
N U M E R I C A L RESULTS
Formulas (8) and (9) are used in all of the numerical simulations. Figures 4 - 7 were generated by the package Differential Systems in which a fifth-order Runge-Kutta method was used. The parameter values corresponding to each figure are presented in Table 1. Figure 4 shows several trajectories of the biomass and the amount of toxicant in the biomass for SYSIII. We are able to see that the trajectories of the system are converging to the asymptotically stable equilibrium point (K, BK). Figure 5 is a phase portrait of SYSIII that also shows orbits converging to the equilibrium value. Figure 6 displays trajectories of the biomass and the total amount of toxicant in the biomass for the nonautonomous system, SYSII. As our analysis from Section 5 indicated, the trajectories of SYSII are converging to the stable equilibrium point of the limiting system SYSIII. Further trajectories for SYSII of the biomass are shown in Figure 7 for varying parameter values. Here again we numerically see the asymptotic results of Section 5. Figure 8 is a numerical simulation of several trajectories of SYSI. The solid line represents biomass, the dashed line represents the concentration of toxicant in the environment, and the dash-dot line represents the amount of toxicant in the biomass. Figure 9 is a graph of the control function versus time. Since x, Ce, and C all converge, we have that u becomes relatively constant for large time.
A C O N T R O L P R O B L E M IN A P O L L U T E D E N V I R O N M E N T
t.r
200
490
157
590
800
1000
500
80g
100|
]
(a) 200
460
-i
i"
~T
.~ 21
(b) FIG. 4. Several trajectories of the biomass and the amount of toxicant in the biomass for SYSIII.
158
DIANA MARIA THOMAS ET AL. 9
I
Z
3
4
5
5
7
8
g
Ill
'¢1
!
"
'
,
BLomasa
FIG. 5. Phase portrait of SYSIII.
Figure 10 is a plot of the biomass versus time in the absence of toxicant modeled by logistic growth (dashed line) and a plot of the biomass versus time in the presence of toxicant modeled by SYSI. We are able to see that the effect of the toxicant on the biomass is a lowered carrying capacity. This inhibited carrying capacity effect also occurred in experimental data as seen in Figure 1. 7.
CONCLUSIONS
From the theoretical analysis of the model it was found that if the toxicant input into the environment is appropriately controlled, the biomass will tend to a new lowered carrying capacity. The value of this new carrying capacity can be determined from the water quality criteria. It was also shown that in cases where the population was stressed before the toxicant was introduced, prohibition of the toxicant may be the only solution to preserve the population. The analytical results are supported by experimental data. Lowered carrying capacity was seen in rotifer populations that were exposed to low concentrations of copper. There is much to be done for future work. A more realistic control problem should require that C E ~< M. There are extensive results on
A C O N T R O L PROBLEM IN A POLLUTED E N V I R O N M E N T 200
400
159
6OO
800
I BOll
600
800
lOON
qr
mu,)
i
Time
(a) 200
400
t
Time
(b) FIG. 6. Several trajectories of the biomass and the amount of toxicant in the biomass for SYSII.
160
D I A N A M A R I A T H O M A S ET AL 500
800
ION
eee
iooe
Ttme
(a)
(b) FIG. 7. Further trajectories of the biomass with different parameter values.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
161
TABLE 1 Parameter Values Corresponding to Figures 4-10 Parameter Figure 4 5 6 7a 7b 8 9 10
ro
rI
p Ko
0.1 0.1 0.1 0.01 0.01 0.1 0.1 0.1
0.01 0.01 0.01 0.1 0.3 0.01 0.01 0.01
2 2 2 3 3 2 2 2
5 5 5 6 4 5 5 5
a
M
/x
}t
Yo
b
h
g
BK
0.1 0.1 0.1 0.01 0.001 0.01 0.1 0.1
0.5 0.5 0.5 0.4 0.5 0.5 0.5 0.5
0.3 0.3 0.3 0.2 0.2 0.2 0.3 0.3
0.01 0.01 0.01 0.1 0.1 0.01 0.01 0.01
-3.1 -3.1 -3.1 -0.29 -0.51 NA -3.1 -3.1
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
4.57 4.57 4.57 5.94 4.008 4.93 4.57 4.57
0.762 0.762 0.762 0.12 0.01 0.123 0.762 0.762
these types of c o n s t r a i n t s in control theory literature. Also, p e r h a p s a n a g e - d e p e n d e n t m o d e l or delay m o d e l should be i n c o r p o r a t e d to a c c o u n t for the oscillation in the rotifer b i o m a s s growth that is seen in Figure 1. Again, there are m a n y resources available for this type of p r o b l e m . A n adaptive p a r a m e t e r c o n t r o l p r o b l e m should be investigated b e c a u s e the p a r a m e t e r s are n e v e r m e a s u r e d exactly. T h e r e are m a n y t e c h n i q u e s available to adjust for errors in m e a s u r e m e n t s in the r e a l m of control
r
\ .
.
.
.
.
.
200
.
.
.
.
.
.
400
.
.
.
.
.
.
600
.
.
.
.
m "
800
| I000
FIG. 8. Tr~ectories of SYSI. Solid curves are x; dashed curves are dashed-dot curves are C.
time
CE; and
162
DIANA MARIA THOMAS ET AL.
0,12
0.1
0.08
0,06
0.04
O. 02
100
200
300
400
500
600
700
FiG. 9. Plot of the control function u versus time.
x 5
3
2
1'
lOO
200
300
400
soo
~oo
700
FIG. 10. Plots of biomass behavior in the absence of toxicant (dashed curve) and in the presence of toxicant (solid curve) versus time.
A CONTROL PROBLEM IN A POLLUTED ENVIRONMENT
163
theory. Finally, o n e should c o n s i d e r the o p t i m a l c o n t r o l p r o b l e m , w h e r e we m i n i m i z e the cost associated with controlling the p o l l u t i o n in the environment.
We thank Almut Burchard and Hugo Leiva for the time spent discussing the material and for their suggestions. We also thank Dr. H. I. Freedman for his insight on modifications for our model. REFERENCES 1 R.L. Wallace and T. W. Snell, Rotifera, in Ecology and Systematics of North American FreshwaterInvertebrates, J. H. Thorp and A. P. Covich Eds., Academic, New York, 1991, pp. 187-248. 2 T. W. Snell and B. D. Moffat, A two day life cycle test with the rotifer Brachionus calyciflorus, Environ. Toxicol. Chem., 11:1249-1257 (1992). 3 C.R. Janssen, G. Personne, and T. W. Snell, Ecotoxicological studies with the freshwater rotifer Brachionus calyciflorus. VI. Short chronic toxicity tests with rotifers, Aquat. Toxicol., in press (1994). 4 T.W. Snell and C. Janssen, Rotifers in ecotoxicology: a review, Hydrobiologia, in press (1994). 5 C. Janssen, The use of sublethal criteria for toxicity test with the freshwater rotifer Brachionus ca~ciflorus (Pallas), Ph.D. Thesis, Univ. Ghent, Belgium, 1992. 6 T.G. Hallam, C. E. Clark, and G. S. Jordan, Effects of toxicants on populations: a qualitative approach. II. First order kinetics, J. Math. Biol., 18:25-37 (1983). 7 T.G. Hallam and C. E. Clark, Effects of toxicants on populations: a qualitative approach. I. Equilibrium environmental exposure, Ecol. Modeling, 18:291-304 (1983). 8 S. Barnett and R. G. Cameron, Introduction to Mathematical Control Theory, Oxford Univ. Press, New York, 1985, pp. 351-379. 9 E.B. Lee and L. Markus, Foundations of Optimal Control Theory, Wiley, New Yori, 1967, p. 24. 10 J.-J. E. Slotine and W. Li, Applied Nonlinear Control, Prentice-Hall, Englewood Cliffs, NJ, 1967, p. 207. 11 H.K. Khalil, Nonlinear Systems, Macmillan, New York, 1992. 12 R.H. Martin and M. Pierre, Nonlinear Reaction Diffusion Systems--Nonlinear Equations in the Applied Sciences, Academic, New York, 1992. 13 H.R. Thieme, Asymptotically autonomous differential equations in the plane, Rocky Mt. J. Math., 21:351-379 (1994).