Applied Numerical Mathematics 59 (2009) 845858

www.elsevier.com/locate/apnum
An application of optimal control theory to river pollution
remediation
L.J. Alvarez-Vzquez
a,
, A. Martnez
a
, M.E. Vzquez-Mndez
b
, M.A. Vilar
b
a
Departamento de Matemtica Aplicada II, ETSI Telecomunicacin, Universidad de Vigo, 36310 Vigo, Spain
b
Departamento de Matemtica Aplicada, Escola Politcnica Superior, Universidad de Santiago de Compostela, 27002 Lugo, Spain
Available online 29 March 2008
Abstract
The main goal of this work is to use mathematical modelling and numerical optimization to obtain the optimal purication of
a polluted section of a river. The most common strategy consists of the injection of clear water from a reservoir in a nearby point.
In this process, the main problem consists, once the injection point is selected by geophysical reasons, of nding the minimum
quantity of water which is needed to be injected into the river in order to purify it up to a xed level: this will be the aim of
this paper. We formulate this problem as a hyperbolic optimal control problem with control constraints, where the state system
is given by the 1D shallow water equations coupled with the pollutant concentration equation, the control is the ux of injected
water, and the objective function is related to the total quantity of injected water and the pollution thresholds. We present the
mathematical formulation of the environmental problem, deriving a rst order optimality condition in order to characterize the
optimal solutions (via the introduction of the adjoint system for the computation of the gradient of the cost function). Finally,
we deal with the numerical resolution of a realistic problem, where a nite element/nite difference discretization is used, an
optimization algorithm is proposed, and computational results are provided.
2008 IMACS. Published by Elsevier B.V. All rights reserved.
MSC: 49K20; 76D55; 93C20
Keywords: Numerical optimization; Optimal control; Remediation; River pollution
1. Introduction
During last decades rivers have been changed into the main receivers for the disposal of wastewater generated
by the growing industrialization and the population explosions in urban areas. Wastewater discharges in our rivers
have increased up to a so alarming level, that they are nowadays one of the more serious water pollution problems.
Depending on its source, wastewater discharges in rivers are classied under three main types: domestic wastewater
(those coming from a purifying plant which collect the water from an urban sewer system), industrial wastewater
(those coming from an industrial plant), and agricultural wastewater (those coming from irrigation and fertilizing). In
the last times, very strict legislative requirements concerning to the wastewater disposal in rivers have been established
*
Corresponding author.
E-mail addresses: lino@dma.uvigo.es (L.J. Alvarez-Vzquez), aurea@dma.uvigo.es (A. Martnez), ernesto@usc.es (M.E. Vzquez-Mndez),
miguel@usc.es (M.A. Vilar).
0168-9274/$30.00 2008 IMACS. Published by Elsevier B.V. All rights reserved.
doi:10.1016/j.apnum.2008.03.027
846 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
by the governments of developed countries, so that all wastewater discharged into a river must rst be treated in a
purifying plant, in order to reduce its pollutants levels.
Despite all these legal regulations, many of our rivers (or, specically, several sections of these rivers) keep levels of
pollutants higher than the allowed thresholds, since the river basins are not capable of assimilating all the wastewater
disposed there. The most common practice for the remediation of these polluted river sections consists of the injection
of a given amount of clear water from a reservoir in a nearby point. This strategy of ow regulation by means of
reservoirs presents a high efciency in order to purify polluted rivers in a short period of time. In this process of
increasing the river ow by controlled releases from reservoirs, the main problem consists (once the injection point is
chosen by geophysical reasons) of nding the minimum quantity of water which is needed to be injected into the river
section in order to purify it up to a xed level.
In this paper we make use of the powerful tools of the mathematical modelling and the optimal control of partial
differential equations in order to determine this minimal quantity of injected water such that we can be able to ensure
that the contaminant concentration in the river section is lower than a xed threshold. As far as we know, this is the
rst time that this environmental problem is studied from a rigorous mathematical point of view, but the authors have
previously analyzed and solved a related problem (the water conveyance problem [2]) for the case of a boundary
control problem posed on a two-dimensional area of shallow water. The only (and early) attempt addressed to deal
with a very simplied related problem can be found in the early works of Bogobowicz and Sokolowski [10,9] where
the water quality in a river is analyzed in terms of dissolved oxygen (DO) and biochemical oxygen demand (BOD),
but always from a completely linear point of view, that is, the state system is a linearization of the shallow water
equations, which avoids the main difculties of the problem. Among other works devoted to similar problems we
have also to cite, for instance, those of Lee and Leitmann [17] or Zheng et al. [25] where DO/BOD are studied in the
framework of the optimal control of delay systems of ordinary differential equations. In the eld of optimal location
and identication of pollution sources in a river we can mention, among others, the interesting works of Smith [23],
Piasecki and Katopodes [20,21], Kernevez [15] or Bhal et al. [8], but all of them formulated in terms of water height
instead of wet section, as is our case.
Our approach, here applied to a fully one-dimensional problem, can be easily extended to the case of 2D shallow
water equations, allowing us to analyze the optimal remediation of a lake or a bay, for instance. In this case, purifying
water could be injected trough a small part of the boundary, and the problem should be formulated as a boundary
control problem with the ow rate through that part of the boundary as the design variable. The rst steps in this
direction were given by the authors in their works studying the optimal purication of polluted regions of an estuary
by means of clear water injection [2], or the optimal management of a river shway [4].
So, in Section 2 we formulate our environmental problem as a hyperbolic optimal control problem with control
constraints, where the state system is given by the one-dimensional inviscid shallow water equations (formulated in
terms of averaged velocity and wet section) coupled with the pollutant concentration equation (we must remark that
in this work we deal with a unique generic pollutant for instance, in an urban scenario the most commonly used
are pathogenic microorganisms or chlorides but it is direct its generalization to several pollutants). Moreover, the
control is given by the ux of injected water, and the objective function is related to the total quantity of injected water
and the satisfaction of the pollution constraints.
Section 3 is devoted to the derivation of a rst order optimality condition, necessary in order to characterize the
optimal solutions, which will be obtained via the introduction of the adjoint system for the computation of the approx-
imate gradient of the cost function.
In Section 4 we deal with the numerical resolution of the previous problem, where a nite element/nite difference
discretization is used, and an optimization algorithm is proposed. Finally, the computational results for a realistic
problem, and the conclusions are provided in the latter sections.
2. Mathematical formulation of the environmental problem
Let us take a river L meters in length, as given in Fig. 1, and consider E tributaries owing into the river (located
at points e
1
, . . . , e
E
), V wastewater discharges (located at points v
1
, . . . , v
V
) coming from purifying plants, and a
point p where clear water is discharged from a nearby reservoir.
Having in mind that we want to control pollution in the river section corresponding to [p, L] for a time interval
of T seconds, and since we are going to consider only one-dimensional changes along the direction of ow in the river,
L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858 847
Fig. 1. General diagram of a river.
for each (x, t ) [0, L] [0, T ] let us denote by A(x, t ) the area of the section occupied by water (usually known as
wet section) that is assumed to remain positive for any time t [0, T ]; let us denote by u(x, t ) the average velocity
in the wet section x meters from the source and t seconds from the moment the control is initiated; let us denote
by q(x, t ) the ow rate across the section (that is, q(x, t ) = A(x, t )u(x, t )); and let us denote by c(x, t ) the quantity
of a generic pollutant in the wet section measured in units per meter (that is, c(x, t ) = A(x, t )(x, t ), with (x, t )
the averaged pollutant concentration). The evolution of the wet area A(x, t ), the ow rate q(x, t ) and the quantity of
pollutant c(x, t ) is given by the following hyperbolic initial-boundary value problem (see, for instance, [5] or [7]):
A
t
+
q
x
=Q(x p) +
E

j=1
q
j
(x e
j
) +
V

k=1
p
k
(x v
k
)
. ,, .
=g
1
(x,t )
in (0, L) (0, T ),
q
t
+

x
_
q
2
A
_
+gA

x
=QW cos( )(x p)
+
E

j=1
q
j
U
j
cos(
j
)(x e
j
) +
V

k=1
p
k
V
k
cos(
k
)(x v
k
) +S
f
. ,, .
=g
2
(x,t )
in (0, L) (0, T ),
c
t
+

x
_
qc
A
_
+kc =
E

j=1
n
j
(x e
j
) +
V

k=1
m
k
(x v
k
)
. ,, .
=g
3
(x,t )
in (0, L) (0, T ),
A(L, t ) =A
L
(t ) in [0, T ],
q(0, t ) =q
0
(t ) in [0, T ],
c(0, t ) =c
0
(t ) in [0, T ],
A(x, 0) =A
0
(x) in [0, L],
q(x, 0) =q
0
(x) in [0, L],
c(x, 0) =c
0
(x) in [0, L],

(1)
where:
(x b) denotes de Dirac measure at a generic point b [0, L].
for j = 1, 2, . . . , E, e
j
(0, L) is the point where the mouth of the j-th tributary is located, q
j
(t ) is the corre-
sponding ow rate, U
j
(t ) is its velocity,
j
is the angle between the j-th tributary and the main river, and n
j
(t )
is its mass pollutant ow rate.
for k = 1, 2, . . . , V , v
k
(0, L) is the point where the k-th wastewater discharge is located, p
k
(t ) is the corre-
sponding ow rate, V
k
(t ) is its velocity,
k
is the angle between the k-th discharge and the river, and m
k
(t ) is its
mass pollutant ow rate.
p (0, L) is the point where clear water is discharged, Q(t ) is the corresponding ow rate (which will be our
control), W(t ) is its velocity, and is the angle between the discharge and the river. (We must also remark that,
since we are injecting clear water, this term does not appear in the second member of the pollutant equation.)
848 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
Fig. 2. Longitudinal cut of a river at time t .
g is gravity.
S
f
denotes the bottom friction stress, which is given, for instance, by the Chzy law.
(x, t ) =H(x, t ) +b(x) is the height of water with respect to a xed reference level (see Fig. 2), where H(x, t )
represents the height of the water column and b(x) geometrically describes the river bottom.
k is the loss rate for pollutant.
Remark 1. At a rst glance at system (1), four unknowns can be detected: A(x, t ), q(x, t ), (x, t ), and c(x, t ).
However, it is obvious that, if the geometry of the river is known, A(x, t ) can be calculated from (x, t ). In effect, for
each x [0, L], the geometry of the river gives us a smooth, strictly increasing and positive function S(, x) verifying
S(0, x) = 0 and
S
_
H(x, t ), x
_
=A(x, t ) in [0, L] [0, T ]. (2)
In this way, since we are dealing with a system of conservation laws with source terms (also known as balance laws)
whose conservative variables are A, q and c, if we write in terms of A, the non-conservative unknown can be
suppressed in (1). Explicitly, we have

x
(x, t ) =

x
_
B
_
A(x, t ), x
__
+b

(x), (3)
where, for each x [0, L], B(, x) denotes the inverse of the function S(, x), that is, B(A(x, t ), x) =H(x, t ), which
will also be a smooth, strictly increasing and positive function.
Moreover, it can be shown (cf. [7]) that
A(x, t )

x
(x, t ) =

x
_
G
_
A(x, t ), x
__
F
_
A(x, t ), x
_
+A(x, t )b

(x), (4)
where
G(A, x) =
B(A,x)
_
0
S(r, x) dr,
F(A, x) =
B(A,x)
_
0
S
x
(r, x) dr =
A
_
0
B
x
(s, x) ds.
Remark 2. Bearing in mind (4), the state system (1) can be written in the standard notation for conservation laws:
U
t
+

x
K(U, x) =C(U, x, Q), (5)
where
U =
_
A
q
c
_
, K(U, x) =

q
q
2
A
+gG(A, x)
qc
A

,
C(U, x, Q) =
_
Q(x p) +g
1
gF(A, x) gAb

+QW cos( )(x p) +g

2
kc +g
3
_
.
L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858 849
Since the Jacobian matrix of K
DK(U, x) =

0 1 0

q
2
A
2
+g
G
A
(A, x) 2
q
A
0

qc
A
2
c
A
q
A

has the three eigenvalues

1
(U, x) =
q
A
(A, x),
2
(U, x) =
q
A
+(A, x),
3
(U, x) =
q
A
,
where (A, x) =
_
g
G
A
(A, x) is the celerity of the system, and since (from the denition of G in previous remark)
G
A
(A, x) =S
_
B(A, x), x
_
B
A
(A, x) =A
B
A
(A, x) >0 (6)
due to the strictly increasing character of B, we have that the system of conservation laws (5) is strictly hyperbolic,
with the corresponding right eigenvectors given, respectively, by
R
1
(U, x) =
_
1
q
A
(A, x)
c
A
_
, R
2
(U, x) =
_
1
q
A
+(A, x)
c
A
_
, R
3
(U, x) =
_
0
0
1
_
.
One can easily check that the third characteristic eld is linearly degenerate (since R
3
(U, x) D
3
(U, x) = 0)
and, consequently, associated with contact waves. However, the rst and second characteristic elds are genuinely
non-linear (in the sense that R
k
(U, x) D
k
(U, x) = 0, k = 1, 2) and related to shock or rarefaction waves (for more
details see, for instance, the classic books [11,12,22] or [24], and the references therein).
Since the state system presents Dirac measures in the source term, and its solution may develop discontinuities
after a nite time, weak solution should be considered instead of classical solutions (for the study of the effect of the
Dirac measures in the source term, the interested reader is referred to the detailed analysis of Diehl [13]). Thus, we
will say that a bounded function U L
p
(0, T ; L
p
(0, L))
3
, for 1 p , is a weak solution of (5) if it satises (5)
in the sense of distributions, that is, if
T
_
0
L
_
0
_
U

t
+K(U, x)

x
+C(U, x, Q)
_
= 0
holds for all D((0, L) (0, T ))
3
. In this case, initial and boundary conditions will be naturally interpreted in the
sense of traces. We must note that any weak solution which is regular enough is necessarily a regular solution (see,
for instance, Dafermos [12] or Lu [18]).
Recalling the mathematical formulation of the environmental problem, by technological reasons we are led to
consider only the positive uxes in the set of admissible controls
U
ad
=
_
Q L
2
(0, T ): 0 QQ
max
_
, (7)
since we are just injecting (not extracting) clear water, and the quantity of injected water must be bounded.
Then, we formulate the control problem considering as the cost functional the total amount of clear water injected
through the point p together with a measure in the region of the river starting from point p of the contaminant
concentration which remains higher than the xed threshold c
max
. Thus, we dene the following cost function (similar
to the objective function used by the authors in the related work [2]):
J(Q) =

2
T
_
0
Q(t )
2
dt +

2
T
_
0
L
_
p
_
c(x, t ) c
max
_
2
+
dx dt (8)
where and are two weight parameters, and (c c
max
)
+
denotes the positive part of c c
max
, that is,
(c c
max
)
+
= max{c c
max
, 0}.
850 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
So the problem, denoted by (P), of the optimal water injection for the purication of a polluted section in a river
consists of nding the control ux Q U
ad
of injected clear water in such a way that, verifying the state system (1),
minimizes the cost function J given by (8). Thus, the problem can be written as:
(P) min
QU
ad
J(Q).
3. Analysis of the optimal control problem
The question regarding the existence of solution for problem (P) is still an open problem and will not be discussed
here. Moreover, the problem will be non-convex because of the non-linearity of the state system, so uniqueness of
solution is not expected.
In this section we will center our attention in obtaining a formal rst-order optimality condition satised by the
solutions of problem (P). In order to express this necessary optimality condition in the simplest way we introduce the
functions (w, P, s) solutions of the adjoint system (see, for instance, [20]):

w
t
+
q
2
A
2
P
x
+
qc
A
2
s
x
g
B
A
(A, x)
(AP)
x
= g

x
P in (0, L) (0, T ),

P
t
2
q
A
P
x

w
x
=
c
A
s
x
in (0, L) (0, T ),

s
t

q
A
s
x
+ks =
[p,L]
(c c
max
)
+
in (0, L) (0, T ),
_
q
2
A
2
P +
qc
A
2
s g
B
A
(A, x)AP
_
(0, t ) = 0 in [0, T ],
_
w +2
q
A
P +
c
A
s
_
(L, t ) = 0 in [0, T ],
_
q
A
s
_
(L, t ) = 0 in [0, T ],
w(x, T ) = 0 in [0, L],
P(x, T ) = 0 in [0, L],
s(x, T ) = 0 in [0, L],

(9)
where
A
denotes the indicator function of a generic set A, that is,

A
(x) =
_
1, if x A,
0, otherwise.
Remark 3. The adjoint system (9) can be rewritten in the standard notation for linear transport systems:
V
t
+

A(U, x)
V
x
=

B(U, x)V +

C(U, x), (10)
where
V =
_
w
P
s
_
,

A(U, x) =

0
q
2
A
2
+g
B
A
(A, x)A
qc
A
2
1 2
q
A
c
A
0 0
q
A

B(U, x) =
_
0 g

x
g
B
A
(A, x)
A
x
0
0 0 0
0 0 k
_
,

C(U, x) =
_
0
0

[p,L]
(c c
max
)
+
_
.
Obviously, from (6), matrix

A(U, x) is the transpose of DK(U, x) and, consequently, it has the same three eigen-
values that for the state system (that is,
1
(U, x),
2
(U, x) and
3
(U, x)), which are independent on V .
L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858 851
Now, assuming the solvability of the state system (1) and of the adjoint system (9), we can prove the following
result:
Theorem 1. Let Q U
ad
be a solution of the control problem (P). Let (A, q, c) and (w, P, s) be, respectively, the
corresponding solutions of the state system (1) and of the adjoint system (9). Then, the following relation is veried:
T
_
0
_
Q(t ) w(p, t ) P(p, t )W(t ) cos( )
__
R(t ) Q(t )
_
dt 0, R U
ad
. (11)
Proof. Since Q is a solution of the minimization problem (P), the following inequality holds:
DJ(Q) (R Q) 0, R U
ad
. (12)
Let (A, q, c) be the state corresponding to the optimal control, then we have:
DJ(Q) (R Q) =
T
_
0
Q(R Q) +
T
_
0
L
_
p
(c c
max
)
+
c,
where (

A, q, c) =
D
DQ
(A, q, c)(Q) (R Q) is given by the linearized system:

A
t
+
q
x
=(R Q)(x p) in (0, L) (0, T ),
q
t
+

x
_
2q q
A
_

x
_
q
2
A
2

A
_
+g

A

x
= gA

x
_
B
A
(A, x)

A
_
+(R Q)W cos( )(x p) in (0, L) (0, T ),
c
t
+

x
_
qc
A
_
+

x
_
q c
A
_

x
_
qc
A
2

A
_
+k c = 0 in (0, L) (0, T ),

A(L, t ) = 0 in [0, T ],
q(0, t ) = 0 in [0, T ],
c(0, t ) = 0 in [0, T ],

A(x, 0) = 0 in [0, L],

q(x, 0) = 0 in [0, L],
c(x, 0) = 0 in [0, L].

(13)
Thus, we have:
DJ(Q) (R Q) =
T
_
0
Q(R Q) +
T
_
0
L
_
0

[p,L]
(c c
max
)
+
c
=
T
_
0
Q(R Q) +
T
_
0
L
_
0
_
s
t
+
q
A
s
x
ks
_
c
=
T
_
0
Q(R Q)
T
_
0
L
_
0
s
_
c
t
++

x
_
q c
A
_
+k c
_
=
T
_
0
Q(R Q) +
T
_
0
L
_
0
s
_

x
_
qc
A
_

x
_
qc
A
2

A
__
852 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
=
T
_
0
Q(R Q) +
T
_
0
L
_
0
_
s
x
qc
A
2

A
s
x
qc
A
_
+
T
_
0
_
s
qc
A
_
(L, t ) +
T
_
0
_
s
qc
A
2

A
_
(0, t )
=
T
_
0
Q(R Q) +
T
_
0
L
_
0

A
_
w
t

q
2
A
2
P
x
+g
B
A
(A, x)
(AP)
x
g

x
P
_
+
T
_
0
L
_
0
q
_
P
t
+2
q
A
P
x
+
w
x
_
+
T
_
0
_
s
qc
A
_
(L, t ) +
T
_
0
_
s
qc
A
2

A
_
(0, t )
=
T
_
0
Q(R Q) +
T
_
0
L
_
0
_

q
t
P

x
_
2
q
A
q
_

q
x
w
_
+
T
_
0
L
_
0
_

A
t
w +
P
x
_
q
2
A
2

A
_
P g

x
_
B
A
(A, x)

A
_
AP g

AP
_
=
T
_
0
Q(R Q)
T
_
0
L
_
0
(R Q)(x p)w
T
_
0
L
_
0
(R Q)W cos( )(x p)P.
So, we obtain:
DJ(Q) (R Q) =
T
_
0
Q(t )
_
R(t ) Q(t )
_
dt
T
_
0
_
R(t ) Q(t )
_
w(p, t ) dt

T
_
0
_
R(t ) Q(t )
_
W(t ) cos( )P(p, t ) dt
=
T
_
0
_
Q(t ) w(p, t ) P(p, t )W(t ) cos( )
__
R(t ) Q(t )
_
dt. (14)
Now, taking expression (14) to (12) we derive the optimality condition (11).
Remark 4. The optimality condition (11) requires, in order to make sense, that P(x, t ) and w(x, t ) be continuous in
the space variable (so that P(p, t ) and w(p, t ) can be dened), but since we are dealing with the adjoint states P and w
in the space L
2
(0, T ; H
1
(0, L)) and, from the RellichKondrashov Theorem (see, for instance, the classical mono-
graph [1]) we know that H
1
(0, L) is continuously embedded in C([0, L]), this requirement is fully accomplished.
4. Numerical resolution of the discretized problem
In order to compute the value of the cost function J at any control Q, we rstly need to obtain the numerical solution
of the state system (1), that is, we have to compute the values of the state variables A, q and c. Since variable c is
not coupled with the two rst equations of system (1), we will proceed to solve it sequentially: the rst step in the
numerical resolving of this problem will be to obtain, by solving the two rst equations of system (1), the wet section
A(x, t ) and the ow rate q(x, t ) in every point x and at every time t . As a second step, we will proceed to solve the
third equation of system (1), which will give us the quantity of pollutant c(x, t ), necessary in order to compute the
value of the cost function. Finally, the third step consists of computing and minimizing the value of J(Q).
Due to technological reasons (ow control mechanisms cannot act upon water ow in a continuous way, but discon-
tinuously at short time periods) we choose to seek the admissible controls Q among the piecewise-constant L
2
(0, T )
L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858 853
functions. So, for the time interval [0, T ] we choose a number K N, we consider the time step = T/K and we
dene the discrete times
m
= m for m = 0, 1, . . . , K. Thus, a function Q U
ad
which is constant at each subin-
terval determined by the grid {
0
,
1
, . . . ,
K
} is completely xed by the set of values Q

=(Q
0
, Q
1
, . . . , Q
K1
)
[0, Q
max
]
K
R
K
, where Q
m
=Q(
m
), m= 0, . . . , K 1.
Thus, to solve the state system (1) corresponding to a given control Q, we consider the following discretization:
for N N given, we dene t = T/N and take t
n
= nt , for n = 0, 1, . . . , N. Moreover, we choose
h
= {x
0
=
0, x
1
, x
2
, . . . , x
M
=L} a partition of interval [0, L] in M subintervals I
k
= [x
k1
, x
k
], k = 1, 2, . . . , M, such that there
exists P {1, 2, . . . , M 1} verifying x
P
=p.
4.1. Computation of the discretized A and q
We begin by making an implicit time discretization of the rst equation of system (1). It leads us to write, for
n = 0, 1, . . . , N 1
A(x, t
n+1
) A(x, t
n
) +t
_
Q(t
n+1
)(x p) +g
1
(x, t
n+1
)
q
x
(x, t
n+1
)
_
in (0, L). (15)
For the time discretization of the second equation of system (1) we are going to use the method of characteristics
(see [6]), which stems from considering the following equality
D(Vq)
Dt
(x, t ) =
q
t
(x, t ) +
(uq)
x
(x, t ), (16)
where
D(Vq)
Dt
(x, t ) =

_
V(x, t ; )q
_
X(x, t ; ),
__
|=t
,
with
X(x, t ; ) the characteristic line, providing the position at time of the particle that occupied the position x at
time t . Thus, the characteristic line is the unique solution of the following ordinary differential equation:
dX
d
(x, t ; ) =u
_
X(x, t ; ),
_
X(x, t ; t ) =x

(17)
V(x, t ; ) the evolution of the element of volume, which is given by the solution of the following ordinary differ-
ential equation:
dV
d
(x, t ; ) =
u
x
_
X(x, t ; ),
_
V(x, t ; )
V(x, t ; t ) = 1

(18)
Taking into account (3), the expression (16) for uq =q
2
/A changes the second equation of system (1) into
D(Vq)
Dt
(x, t ) +gA(x, t )
_
B
x
_
A(x, t ), x
_
+b

(x)
_
=Q(t )W(t ) cos( )(x p) +g
2
(x, t ) in (0, L) (0, T ).
If we denote X
n
(x) = X(x, t
n+1
; t
n
) and V
n
(x) = V(x, t
n+1
; t
n
), the characteristics method puts forward the ap-
proximation
D(Vq)
Dt
(x, t
n+1
)
q(x, t
n+1
) q(X
n
(x), t
n
)V
n
(x)
t
, (19)
which leads us to the following discretization for the second equation of system (1): for each n = 0, 1, . . . , N 1
q(x, t
n+1
)
t
+gA(x, t
n+1
)
_
B
x
_
A(x, t
n+1
), x
_
+b

(x)
_
=
q(X
n
(x), t
n
)V
n
(x)
t
+Q(t
n+1
)W(t
n+1
) cos( )(x p) +g
2
(x, t
n+1
) in (0, L). (20)
854 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
Remark 5. Now, if we use the approximation of A(x, t
n+1
) given by (15), the unique unknown of equation (20) is
q(x, t
n+1
). So, the problem given by (20) and the boundary conditions q(0, t
n+1
) =q
0
(t
n+1
), A(L, t
n+1
) =A
L
(t
n+1
)
can be solved, for example, by using Lagrange P
1
nite elements, in order to obtain an approximate q
n+1
h
(x) Y
h
=
{y C
0
([0, L]) such that y
|I
k
P
1
, k = 1, 2, . . . , M} (see [7] for further details).
With the method above described, we obtain, for each n = 0, 1, . . . , N, the functions q
n
h
, A
n
h
verifying q
n
h
(x)
q(x, t
n
), A
n
h
(x) A(x, t
n
) in [0, L]. Thus, we can summarize the method in the following diagram:
Step 0: Take k = 0, q
0
h
(x) =q
0
(x), A
0
h
(x) =A
0
(x).
Step 1: Taking
q
k
h
(x)
A
k
h
(x)
u(x, t
k
) as data, compute X
k
h
(x) X
k
(x) and V
k
h
(x) V
k
(x) by using the backward
Euler scheme to solve (17) and (18).
Step 2: Taking q
k
h
(x), A
k
h
(x), X
k
h
(x) and V
k
h
(x) as data, solve (20) and obtain q
k+1
h
(x) q(x, t
k+1
) (see Re-
mark 5).
Step 3: Taking q
k+1
h
(x) and A
k
h
(x) as data, compute A
k+1
h
(x) A(x, t
k+1
) from (15).
Step 4: If k +1 =N stop, else k =k +1 and go to Step 1.
4.2. Computation of the discretized c
The third equation of system (1) can be solved by using an implicit upwind nite difference scheme, which does
not need to solve any linear equations system. In order to do it, because of the Dirac measures characterizing the
sources, we consider the following approximations: for each k = 0, 1, . . . , M, we dene
hk
: [0, L] [0, +) by

hk
(b) =

b x
k1
(x
k
x
k1
)
2
if b [x
k1
, x
k
],
x
k+1
b
(x
k+1
x
k
)
2
if b [x
k
, x
k+1
],
0 otherwise,
where
hk
(b) (x
k
b) for a generic point b [0, L].
So, taking {c
i
0
=c
0
(t
i
), i = 0, 1, . . . , N} and {c
0
j
=c
0
(x
j
), j = 0, 1, . . . , M} as data, for each n = 0, 1, . . . , N 1
and for each k = 1, . . . , M we compute c
n+1
k
from the following expression:
c
n+1
k
c
n
k
t
+
(q
n+1
h
(x
k
)/A
n+1
h
(x
k
))c
n+1
k
(q
n+1
h
(x
k1
)/A
n+1
h
(x
k1
))c
n+1
k1
x
k
x
k1
+k(t
n+1
, x
k
)c
n+1
k
=
E

j=1
n
j
(t
n+1
)
hk
(e
j
) +
V

j=1
m
j
(t
n+1
)
hk
(v
j
).
Finally, for each n = 0, 1, . . . , N 1, we approach c(x, t
n+1
) by the unique continuous function c
n+1
h
(x) Y
h
verifying c
n+1
h
(x
k
) =c
n+1
k
, for each k = 0, 1, . . . , M.
4.3. Computation of the discretized J
We denote by U

ad
the set of all admissible controls in U
ad
which are piecewise-constant in the partition of the
time interval [0, T ] given by the time step . Then, for any given control Q

ad
, we will need to compute
the value of the cost function J(Q

). By previous steps we have obtained an approximated eld for the quantity of

pollutant c
n
h
, which will be used to approach the objective function value J(Q

) by the following trapezoidal rule

approximation:
L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858 855
J
t
h
_
Q

_
=

2

K1

m=0
(Q
m
)
2
+(Q
m+1
)
2
2
+

2
t
N1

n=0
M1

k=P
x
k+1
_
x
k
(c
n
h
(x) c
max
)
2
+
+(c
n+1
h
(x) c
max
)
2
+
2
dx. (21)
Thus, instead of solving the original minimization problem (P), we are led to solve the following discretized
problem:
(

P) min
Q

ad
J
t
h
_
Q

_
.
4.4. Minimization of the discretized problem
In order to minimize the objective function J
t
h
in the set of admissible controls U

ad
, the standard proposal
consists of using a gradient-type algorithm, where the gradient J
t
h
(Q

) can be directly obtained from a suitable

approximation of expression (14) via the numerical resolution of the adjoint system (9). However, due to the high
computational cost arisen from this proposal, in this paper we alternatively present a gradient-free algorithm for
solving the discretized optimization problem (

P), where the adjoint system will not be involved in. (This algorithm
has already given very good results in several environmental control problems previously studied by the authors [3,4].)
To do this, we will change our discretized problem (

P) into an unconstrained optimization problem by introducing
a penalty function involving the constraints appearing in the denition of the set of admissible controls (7), that is,
Q0 and QQ
max
0.
Thus, we dene the penalty function

J in the following way:

J
_
Q

_
=J
t
h
_
Q

_
+
K1

m=0
max
_
Q
m
, Q
m
Q
max
, 0
_
(22)
where the parameter > 0 determines the relative contribution of the objective function and the penalty terms.
Function

J is an exact penalty function in the sense that, for sufciently large , the solutions of our constrained
problem (

P) are equivalent to the minimizers of function

J in R
K
.
For computing a minimum of this penalty function

J we propose the use of a direct search algorithm: the Nelder
Mead simplex method [19]. This is a gradient-free method, which merely compares function values; the values of
the objective function being taken from a set of sample points (simplex) are used to continue the sampling. A short
description of the above algorithm can be found, for instance, in the nice paper of Kelley [14]. Although the Nelder
Mead algorithm is not guaranteed to converge in the general case, it presents good convergence properties in low
dimensions (cf. [16] for a detailed analysis of its convergence under convexity requirements), which is our case.
Moreover, to prevent stagnation at a non-optimal point, we use a modication proposed by Kelley: when stagnation is
detected, we modify the simplex by an oriented restart, replacing it by a new smaller simplex.
5. Numerical example
In this section we present numerical results obtained by using above method to determine the optimal inow ux
in a river which is L = 1000 m in length, and where we consider E = 3 tributaries, V = 2 domestic wastewater
discharges, and one clear water discharge from a reservoir (the diagram and data can be seen in Fig. 3). Moreover, we
consider a parabolic river bed with a non-constant bottom (see Fig. 4) in such a way that:
A=S(H, x) =
4

H
3
3
,
b(x) =
_
500 x
200
if 0 x 500,
0 if 500 x 1000.
856 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
Fig. 3. Diagram of the river and data for the numerical example.
Fig. 4. River bed for the numerical example.
Both initial and boundary conditions were taken as constant, particularly, A
L
(t ) =
4

125
3
m
2
, q
0
(t ) = 1 m
3
s
1
,
c
0
(t ) = 0 um
1
, A
0
(x) =
4

125
3
m
2
, q
0
(x) = 1 m
3
s
1
and c
0
(x) = 0 um
1
. The time interval to control the pollution
was T = 3600 s. Moreover, the loss rate for pollutant was considered constant (k(x, t ) = 10
4
s
1
), and the bottom
friction stress was neglected (S
f
= 0).
Out of the several numerical experiences developed by the authors, we present here only one example correspond-
ing to the case of K = 4 time subintervals. For the objective function we have chosen the threshold c
max
= 5.5 um
1
,
the bound Q
max
= 25 m
3
s
1
, and the weight parameters = 10
3
, = 10 and = 10
5
. For the time discretization
we have taken N = 6000 (that is, a time step of t = 0.6 s), and for the space discretization we have tried a regular
partition of [0, L] in M = 1000 subintervals (consequently, the clear water inow point p = 700 m corresponds to the
node x
P
=x
700
).
Then, applying the NelderMead algorithm, we have passed, after 962 function evaluations, from an initial random
cost

J = 2006.24 to the minimum cost

J = 861.02, corresponding to the optimal ow rate Q
0
= 9.98 m
3
s
1
, Q
1
=
6.91 m
3
s
1
, Q
2
= 5.45 m
3
s
1
, Q
3
= 3.99 m
3
s
1
. Fig. 5 shows the quantity of pollutant c in the nodes x
750
= 750
and x
950
= 950 along the whole time interval [0, 3600]. These nodes are subsequent to the injection point p = 700 m,
so we can see how, after an initial transitory time, the quantity of pollutant ceases into its natural trend for increasing
(due to the accumulation of waste), and remains almost every time lower than the desired threshold c
max
= 5.5 um
1
.
L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858 857
Fig. 5. Controlled quantity of pollutant in the points x = 750 m (solid line) and x = 950 m (dashed line) along the time interval.
6. Conclusions
In this paper the authors have formulated, analyzed and solved an optimal control problem related to river pollution
control, mainly, the optimal remediation of a polluted river by injection of clear water from a nearby reservoir. Once
the environmental problem is mathematically well posed in terms of wet section, averaged water velocity and quantity
of pollutant, a numerical discretization method is proposed for solving the system of partial differential equations
involved in the modelling. A direct search method (NelderMead algorithm) is proposed for solving the discrete
optimization problem, because of its good performance in other related environmental problem previously analyzed
by the authors. The efciency of the algorithm is conrmed by the numerical experiments developed by the authors,
where good results have been obtained with lower computational efforts than with other gradient-free methods also
implemented by the authors.
Acknowledgements
The research contained in this work was supported by Research Project MTM 2006-01177 of Ministerio de Edu-
cacin y Ciencia (Spain). The authors also thank the anonymous referees for their interesting remarks.
References
[1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975.
[2] L.J. Alvarez-Vzquez, A. Martnez, R. Muoz-Sola, C. Rodrguez, M.E. Vzquez-Mndez, The water conveyance problem: Optimal puri-
cation of polluted waters, Math. Models Meth. Appl. Sci. 15 (2005) 13931416.
[3] L.J. Alvarez-Vzquez, A. Martnez, C. Rodrguez, M.E. Vzquez-Mndez, Numerical optimization for the location of wastewater outfalls,
Comput. Optim. Appl. 22 (2002) 399417.
[4] L.J. Alvarez-Vzquez, A. Martnez, M.E. Vzquez-Mndez, M.A. Vilar, Optimal shape design for shways in rivers, Math. Comput. Simul. 76
(2007) 218222.
[5] E. Audusse, M.-O. Bristeau, Transport of pollutant in shallow water: a two time steps kinetic method, M2AN Math. Model. Numer. Anal. 37
(2003) 389416.
[6] M. Bercovier, O. Pironneau, V. Sastri, Finite elements and characteristics for some parabolic-hyperbolic problems, Appl. Math. Modelling 7
(1983) 8996.
[7] A. Bermdez, R. Muoz-Sola, C. Rodrguez, M.A. Vilar, Theoretical and numerical study of an implicit discretization of a 1D inviscid model
for river ows, Math. Mod. Meth. Appl. Sci. 16 (2006) 375395.
[8] M.G. Bhat, K.R. Fister, S. Lenhart, An optimal control model for the surface runoff contamination of a large river basin, Natur. Resource
Modeling 12 (1999) 175195.
[9] A. Bogobowicz, Theoretical aspects of modeling and control of water quality in a river section, Appl. Math. Comp. 41 (1991) 3560.
858 L.J. Alvarez-Vzquez et al. / Applied Numerical Mathematics 59 (2009) 845858
[10] A. Bogobowicz, J. Sokolowski, Modelling and control of water quality in a river section, in: System Modelling and Optimization, Lecture
Notes in Control and Inform. Sci., vol. 59, Springer, Berlin, 1984, pp. 403414.
[11] A. Bressan, Hyperbolic Systems of Conservation Laws: The One-Dimensional Cauchy Problem, Oxford University Press, Oxford, 2000.
[12] C.M. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer-Verlag, Berlin, 2000.
[13] S. Diehl, A conservation law with point source and discontinuous ux function modelling continuous sedimentation, SIAM J. Appl. Math. 56
(1996) 388419.
[14] C.T. Kelley, Detection and remediation of stagnation in the NelderMead algorithm using a sufcient decrease condition, SIAM J. Optim. 10
(1999) 4355.
[15] J.P. Kernevez, The Sentinel Method and its Application to Environmental Pollution Problems, CRC Press, Boca Raton, FL, 1997.
[16] J.C. Lagarias, J.A. Reeds, M.H. Wright, P.E. Wright, Convergence properties of the NelderMead simplex algorithm in low dimensions, SIAM
J. Optim. 9 (1998) 112147.
[17] C.S. Lee, G. Leitmann, Continuous feedback guaranteeing uniform ultimate boundedness for uncertain linear delay systems: an application to
river pollution control, Comput. Math. Appl. 41 (1988) 929938.
[18] Y.-G. Lu, Hyperbolic Conservation Laws and the Compensated Compactness Method, Chapman & Hall/CRC, Boca Raton, FL, 2003.
[19] J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308313.
[20] M. Piasecki, N.D. Katopodes, Control of contaminant releases in rivers and estuaries. Part I: Adjoint sensitivity analysis, ASCE J. Hydraulic
Eng. 123 (1997) 486492.
[21] M. Piasecki, N.D. Katopodes, Control of contaminant releases in rivers and estuaries. Part II: Optimal design, ASCE J. Hydraulic Eng. 123
(1997) 493503.
[22] D. Serre, Systems of Conservation Laws, Vol. 1, Cambridge University Press, Cambridge, 1999.
[23] R. Smith, Where to put a steady discharge in a river, J. Fluid Mech. 115 (1982) 111.
[24] J. Smoller, Shock Waves and ReactionDiffusion Equations, Springer-Verlag, New York, 1994.
[25] F. Zheng, Q.G. Wang, T.H. Lee, Adaptive robust control of uncertain time delay systems, Automatica J. IFAC 41 (2005) 13751383.