WATER POLLUTION CONTROL RESEARCH SERIES • 16110 EGQ 04/72
EXTENSIONS OF MATHEMATICAL
PROGRAMMING FOR REGIONAL
WATER QUALITY MANAGEMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes the
results and progress in the control and abatement of pollution
in our Nation's waters. They provide a central source of
information on the research, development, and demonstration
activities in the water research program of the Environmental
Protection Agency, through inhouse research and grants and
contracts with Federal, State, and local agencies, research
institutions, and industrial organizations.
Inquiries pertaining to Water Pollution Control Research
Reports should be directed to the Chief, Publications Branch
(Water), Research Information Division, R&M, Environmental
Protection Agency, Washington, B.C. 20460.
-------
EXTENSIONS OF MATHEMATICAL PROGRAMMING
FOR REGIONAL WATER QUALITY MANAGEMENT
by
University of California
Graduate School, of Management
Los Angeles, California 90024
for the
OFFICE OF RESEARCH AND MONITORING
ENVIRONMENTAL PROTECTION AGENCY
Project #16110 EGQ
April 1972
For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402 - Price $1.00
-------
EPA Review Notice
This report has been reviewed by the Environmental
Protection Agency and approved for publication.
Approval does not signify that the contents necessarily
reflect the views and policies of the Environmental
Protection Agency, nor does mention of trade names or
commercial products constitute endorsement or recommenda-
tion for use.
11
-------
ABSTRACT
This report extends the work contained in the earlier work,
"Mathematical Programming for Regional Water Quality Manage-
ment." A mixed integer, continuous variable non-linear pro-
gramming model is developed which promises to be much more
realistic and effective in selecting an optimal configuration
of regional treatment plants and pipe juncture nodes. A new
Bender's type decomposition for the non-linear model is also
presented as an appropriate and feasible solution technique.
In the area of the physical response of the environment, the
earlier work linked the dissolved oxygen profile and the con-
trol options through a constant transfer matrix obtained as a
steady state solution to a system of differential equations.
This work shows how the transfer coefficients can be recalcu-
lated periodically in the course of the system optimization.
Full scale computations are carried out demonstrating this
capability. The theoretical framework for including two addi-
tional treatment methods, mechanical re-aeration and flow aug-
mentation via reservoirs is also developed.' In the social and
economic area, a new proposal for paying for a regional system
is presented. Illustrative computations using the new proposal
are included. This report was submitted in fulfillment of
project 16110EGQ under the partial sponsorship of the Environ-
mental Protection Agency.
111
-------
TABLE OF CONTENTS
Section I.
Section II.
Section III.
Section IV.
Section V.
Section VI.
Section VII.
List of Figures
List of Tables
Conclusions and Recommendations
Introduction
Selection of a Waste Treatment
Plant Configuration
Transfer Functions
Estimation of the Depth and Velocity
of a River for a Given Rate of Flow
Transient Effects in an Estuary
Pollution Taxes and Charges in
Perspective
Section VIII. References
Page
vi
vii
1
3
5
15
37
61
67
81
V
-------
LIST OP FIGURES
FIGURE TITLE PAGE
1 Schematic of General System 6
2 Optimization Process,, Assuming Fixed Transfer 20
Matrix
3 Iterative Procedure for Finding Stable Transfer 21
Matrix
4 Solution of By-Pass Piping Problem Using New 33
Fixed Transfer Matrix
5 Solution of By-Pass Piping Problem Using New 35
Variable Transfer Matrix
6 Schematic and Notations for Spot Inflows 44
7 Schematic and Notations for Derivatives of Spot 45
Inflows
8 Schematic and Notations for Spot Outflows 47
9 Critical, Super-Critical and Sub-Critical Flows as 50
Depicted by the Depth-Specific Energy Relation
10 Schematic and Notations for Continuous Inflow 57
vz
-------
LIST OF TABLES
TABLE TITLE PAGE
1 Net Flows into Section i 23
2 Eddy Exchange Coefficients Between Sections 23
i-1 and i
3 Volumes of Section i 25
4 Reaeration Rate in Section i 25
5 Transfer Matrix 26
6 Thomann's Transfer Matrix 29
7 Transient Period of Section i, Corresponding 65
Starting Conditions and Computational Error
-------
I. CONCLUSIONS AND RECOMMENDATIONS
The benefits to be derived from the use of regional optimiza-
tion models for water quality management would seem to be enor-
mous. Regional optimization models are applicable to bays,
estuaries and rivers. They are capable of dealing with any
quality goal for which an appropriate system of diffusion, diff-
erential equations can be obtained to represent the physical
response of the environment. In particular, they are capable
of dealing with dissolved oxygen, salinity and thermal contamin-
ation. The control options include treatment at the source,
by-pass piping, regional treatment plants, mechanical re-aera-
tion and flow maintenance augmentation. The'output from the
model gives the least cost solutions to achieve the pre-selected
goals using the full range of control options. The optimal
solutions would specify the treatment levels at the major sour-
ces of contamination, the complete regional pipe network, the
location and levels of treatment of an optimal number of region-
al plants, the location and level of mechanical re-aeration and
the location and level of flow augmentation. The model could
be used in conjunction with regional growth of demand projec-
tions to plan the long term treatment capacity expansion. It
could ultimately be used in a real time regional control pro-
cess utilizing physical and chemical treatment as well as the
traditional biological treatment to capitalize on seasonal and
daily flucuation to greatly reduce yearly and daily treatment
costs. This report develops the theoretical possibility of in-
cluding two additional treatment methods, mechanical re-aeration
and flow maintenance or augmentation via reservoirs. In order
to deal with the question of the effect on the re-aeration co-
efficient of radical shifts in flow, Section VI investigates
the effect of flow on the depth and velocity of a river. Fin-
ally, Section VII takes a preliminary look at the transient
effects in an estuary as opposed to the long term steady state
conditions.
In the area of system optimization, Section IV examines the
question of an optimal configuration of regional treatment
plants and pipe junction points. The original model deals with
this question by forcing a solution using only the regional
plants. This forces the plants to operate at sufficiently high
levels as to utilize their economies of scale. When the other
treatment methods are then introduced, they are compelled to
compete against the efficiently operating plants. Where econo-
mically justified they then force out and close the non-compe-
titive plants. This is computationally wasteful as the plants
must first be opened then some of them closed. It is also
mathematically doubtful as to the global optimality of this pro-
cedure. The mathematical foundations of a potentially much
more efficient procedure using an integer non-linear programming
model are presented. Computational effectiveness hinges upon
the power of a new decomposition technique and the capability
of rapidly solving integer programming problems.
- 1-
-------
the realism of the models would seem to be adequate at this
time although they may ultimately be considereably improved.
The computational feasibility for full scale problems has been
demonstrated. It would seem entirely practical to proceed to
actual pilot studies at this time with the tools already at
hand. Further development of the theoretical possibilities
should proceed concomitant with the pilot studies. In parti-
cular computer implementation of the new mixed integer, con-
tinuous variable non-linear model should be undertaken to sub-
stantiate its potential value. The possibility of including
the additional treatment methods of mechanical re-aeration and
flow augmentation should also be exploited. On a longer hori-
zon the dynamic or transient behavior of the physical environ-
ment should be studied more closely. It would be entirely possi-
ble to realize substantial gains from an optimal holding and
discharge policy based on the transient behavior of the environ-
ment. Finally the mathematical programming model for regional
water quality management should be coupled with a regional de-
mand model. Rational planning over a realistic amortization
period requires fairly good estimates of projected demand over
the period. The capital budgeting problem associated with the
optimal timing and amount of treatment capacity expansion should
also be carefully studied.
- 2 -
-------
II. INTRODUCTION
This report continues and extends the work of the earlier re-
port, "Mathematical Programming for Regional Water Quality
Management." The goal of this work is to develop a general
practical planning tool to provide optimal solutions for the
complex choices involved in balancing alternative methods for
attaining water quality goals9 The work encompasses three
broad related areas. The first area is the social and econo-
mic question of what are desirable and feasible water quality
goals and what is the best way to achieve and pay for them.
The second area is the complex technical engineering problem
of the physical response of the environment to various control
strategies. The third area is the mathematical problem of ob-
taining least cost solutions in the total regional water system
that will satisfy the desired goals.
In the area of physical response our earlier work presented a
general mathematical model relating the control options of by-
pass piping, treatment at regional treatment plants and treat-
ment at the source to the-quality goal specified as a dissolved
oxygen profile on an estuary. The original model was limited
in several respects. The coupling of the dissolved oxygen pro-
file goal and the control options was through a constant trans-
fer matrix obtained as a steady state solution to a system of
differential equations first given by Thomann [10]. This is
inadequate when by-pass piping and regional treatment plants
substantially alter the flow. Section IV of this report re-
examines the question of the transfer functions. This work
shows how the transfer coefficients can be re-calculated period-
ically in the course of the system optimization. It further
develops the theoretical possibility of including two additional
treatment methods, mechanical re-aeration and flow augmentation
via reservoirs. In order to deal with the question of the ef-
fect on the re-aeration coefficient of radical shifts in flow,
Section V investigates the effect of flow on the depth and velo-
city of a river. Finally, Section VI takes a preliminary look
at the transient effects in an estuary as opposed to the long
term steady state conditions.
In the area of system optimization. Section III examines the
question of an optimal configuration of regional treatment
plants and pipe junction points. The original model deals
with this question by forcing a solution using only the re-
gional plants. This forces the plants to operate at suffici-
ently high levels as to utilize their economies of scale.
When the other treatment methods are then introduced, they
are compelled to compete against the efficiently operating
plants. Where economically justified, they then force out
and close the non-competitive plants. This is computationally
wasteful as the plants must first be opened then some of them
closed. It is also mathematically doubtful as to the global
optimality of this procedure. The mathematical foundations
- 3 -
-------
of a potentially much more efficient procedure using an integer
non-linear programming model are presented. Computational ef-
fectiveness hinges upon the power of a new decomposition tech-
nique and the capability of rapidly solving integer programming
problems.
In the social and economic area. Section VII discusses economic
policy for regional water quality control. A new proposal for
paying for a regional system is presented. Illustrative compu-
tations using the new proposal are included.
- 4 -
-------
Til. SELECTION OF A WASTE TREATMENT PLANT CONFIGURATION
In our earlier report, "Mathematical Programming for Regional
Water Quality Management/1 the model presented dealt with the
difficult question of the best treatment plant configuration
in a very limited manner. Flow of effluent was restricted flow
direct to a discharge point or flow to a single plant and then
from the plant to a discharge point. There was no flow between
plants and hence no way to achieve say primary treatment for
a large percentage of the effluent at huge, cheap primary plants
and subsequent treatment for a much smaller percentage at smaller
more expensive secondary or tertiary plants. Since pipe juncture
nodes can be considered plants without treatment, the flow pattern
was restricted to single juncture points. This resulted in prac-
tice in manual pipe consolidation and weakened results.
The solution procedure for the original model required forcing
a solution using only the regional plants. This forces the plants
to operate at sufficiently high levels as to utilize their econ-
omies of scale. When the other treatment methods are then intro-
duced, they are compelled to compete against the efficiently
operating plants. Where economically justified they then force
out and close the non-competitive plants. This is computationally
wasteful as the plants must first be opened and driven down and
closed. There are also difficult problems with infinite deriva-
tives in the cost function with the initial very low or zero
flows. These problems can be handled much more effectively using
discrete integer variables in the formulation of the model. This
possibility was not exploited originally becuase of the limita-
tions of the solution techniques for large mixed intefer and con-
tinuous variable non-linear programming problems. Recent research
has greatly expanded our capability to solve integer programming
problems. This extended capability coupled with new and more
powerful decomposition techniques change the picture radically-
It now seems attractive to switch to the mixed integer, continu-
ous variable formulation.
Mixed Integer Regional Treatment Model
In this model we can use a much more extensive and realistic
supporting flow network. This is illustrated by the schematic
shown in Figure 1 on the following page.
The nodes in the network are of three types. A node is a pollu-
ter or source, a sink or a discharge point, or an intermediate
node in the network consisting of a pipe junction or treatment
plant. The sources are fixed at the present sites of the iden-
tified major polluters. The discharge nodes are the presently
used discharge sections of the estuary plus any additional po-
tentially attractive discharge sections. The intermediate nodes
consist of any existing treatment plants or pipe junctures plus
any additional potentially attractive sites. Combining existing
- 5 -
-------
(alternate
s i te s)
(alternate
capacities)
FIGURE 1
SCHEMATIC OF GENERAL SYSTEM
and now sites into an optimal total network is particularly
nice in this approach because the fixed or capital costs can be
separated from the variable or operating costs. This permits
a true comparison of the marginal cost of a new plant to an
existing plant. The arcs in the network are pipes. They can
consist of existing pipes and potentially attractive additional
pipes. Again the fixed and variable costs can be separated.
Also it is possible to have discrete alternate capacities and
enforce a choice of pipe size. Similarly, we can enforce a
- 6 -
-------
choice between alternate plant capacities at a given site or
between alternate plant sites.
In the mathematical statement of the model we use the following:
Notation
Physical Quantities
A - matrix of transfer coefficients
(mg/i per Ib./day)
C - vector of D. 0. goals (mg//)
r. -%removalofB.O.D. at node i
q. • - flow from node i to node j (MGD)
J. - concentration leaving node i (Ib/MG)
m. - mass of B.O.D. leaving node i (Ib/D)
Costs
f. - capital cost of plant i of a given capacity ($)
v. - operating costs of plant i ($/MG)/% removal)
f.. - capital cost of pipe between node i and node j
of given capacity ($)
v.. - operating cost of pipe between node i and node j
of given capacity ($/MG)
p. - cost of additional treatment at source ($/% removal)
Index Sets
R = {j I j is a section or reach of the estuary]
S = {i I node i is a source]
D = {i | node i is a discharge point]
I = {i | node i not a source or sink]
N = S U D u I set of nodes; a c N x N set of arcs
a (x) = {y e N | (x,y) e a] nodes "after" x
b (x) = {y e N (y^x) e a] nodes "before" x
— 7 —
-------
Integer Variables
r0 plant i is closed
Zi ~ *• 1 plant i is open
_ rO pipe (ijj) is closed
ij *-1 pipe (ijj) is open
The mathematical model is:
Subject to
Conservation of Flow
Z q . - Z qk. = 0 i e I
jea(i) 1J keb(i) kl
Z q--=q. i e S (q current discharge MG/D)
jea(i) 13 L
Conservation of Concentration
J. -
keb(i)
m. _
J. = ^- (1 - r.) i e S (m . current discharge Ib)
qi
Capacity at Plant (&. - lower , (j, . - upper capacity)
2i • J&i ^ Z qki ^ zi^i
1 1 keb(i) kl 1 1
(Note that when the plant is closed or not chosen for
the network all flow is shut off and when it is opened
it functions in a restricted range. By breaking up
the range non-convexity can be eliminated and more truly
global optimum achieved. This would require merely add-
ing some:
Alternative Constraints (Plantsj
Z z = 1
These constraints would force a choice of only one of
the plants in the set E..)
_ Q _
-------
Capacity of Pipes (i, . - lower, ^. . - upper capacity)
(Again the pipes can be eliminated or forced to
function in a restricted range. By breaking up
the range non-convexity can be eliminated. Again
we would employ:
Alternative Constraints (Pipes)
2 z. = 1
These constraints would force a choice of one of the
alternatives in the set E...)
The quality goals are enforced by introducing the auxili-
ary variables
m = J • ( 2 q ) t e D
^ *• keb(t) Kt
m = present discharge in Ib's
and the following:
Quality Constraints
2 a., (m -in.) <^ -c. i e R
teD lt: t C
The objective is to minimize the total cost function
minimize
T(zfc,zts,q,r) = PL(zfc,q) + PI(zts,q) + PO(r)
which is broken down into:
Plant Costs
zfc - (ft+vt(r)- kj(t)qkt)
Zts ' (fts
Pipe Casts
Pl(zts^) = ,. Z, a
V t , 5 ) 6A
Polluter Costs
P0(r) =2 r - p
teS
- 9 -
-------
Solution Technique
The foregoing complex mixed integer continuous variable non-linear
programming model presents a truly formidable computational chal-
lenge. Although a direct implicit enumeration of the integer
variables coupled with the bounding off effect of solving the
remaining non-linear programming problems might be attempted,
the outcome of such a venture would be extremely doubtful. Re-
cent advances in Bender's type (see [2]) decomposition offer a
much more promising avenue of attack. A Bender's type decomposi-
tion applicable to a general non- linear programming problem
Subject to
gk(y)
min
is possible. In this type of decomposition a subset of the
variables say the y~ vector can be used to parametrize a re-
laxed Sub-Problem
Subject to
g1(Y1-Y2) £ o
m / — \
mm g (Y1,y2)
and the parameters y_ need not be continuous variables.
Whenever a fixed choice of parametes y9 satisfying the remain-
ing constraints
gk(72) *o
permits a feasible solution to the sub-problem, we obtain a
feasible solution to the complete problem. _In order to avoid
a random search through plausible choices of y~ , we employ a
class of lower bound functions D. In a stepwise procedure we
generate members of this class. By ensuring that at least each
lower bound function already generated indicates a better possi-
ble solution than the best known incumbent feasible solution we
guide the choice of the parameter variables. In point of fact
since it is a necessary condition that each lower bound func-
tion be less than the incumbent when there exists a better solu-
tion, the lower bound functions provide us with a convergence
criteria.
- 10 -
-------
In order to generate the class of lower bound functions of
interest here we require the following expansions of the
functions.
Expansion 1 (With respect to y, )
(i) gi(y,y) =gi(YJy) + vgi(yJy) - AY
Now in order to eliminate the gradients and utilize informa-
tion furnished by the dual variables in the sub-problems, we
must make:
Expansion 2 (With respect to y )
(2) Vg1(y°,y2) - AY]_ = vg1(y°,y°) • AYX
where 2 i
5 g
is a matrix of mixed partial derivatives. Combining the error
terms
(3)
AY2
1
yields the desired representation of each function
(4) g1(y1,y2) = g1(Yj,y2) + vg1(y°Jy°) - AyL + Rx(Ayr *y2)
Using this representation of our functions and the following
relation,
(5) Z x± vgy^y^) = vgm(y°,y°)
which is derived as equation (14) on page 17 of our previous
report, "Mathematical Programming for Regional Water Quality
- 11 -
-------
Management",;yieIds our class of lower bound functions.
The argument employed in establishing this class of lower
bound functions is similar to the argument used to establish
the "Weak Duality Theorem" of linear programming. Any e-
feasible solution of the sub-problems satisfies
(6) gi(y1,y2) = gi(y°;Y2) + vgNy^y^) ' &i + ^(AY^AY;,) ^ el
or
(7)
Using (5) and (7) we demonstrate that as a function of
For any value of y,
(8) gm(y1,y2) = gm(Y^Y2) + vgm(yj,y°) • ^ + Rm(Ay-L
AY
, AY2)
+ 2 x. -e.
For convenience we restate the conditions for a stationary
point of the non-linear algorithm presented in our previous
report [4] on page 19.
Local Linear Stationary Point
a a set of dual variables x 3
P
1) 2 (-x ) > B-. (insufficient resolution)
P P
12 -
-------
a a w and H 3
2) 2 (-x ) • gp(y°) > -gw(y°) - e (Inconsistency)
peH p
a a y e F ^
3) m~^
.2 (-Xi) - g1(y°) > - e (Optimality)
Now choose £-,_ = e/B^ and assume that condition 1) of a local
linear stationary point is not violated then
f Xi >- -Bl and I Xi ' el ^ -£
and when the error expression
i '
we have our lower bound function
(10) gm(y°,y7) + 2(-x ) - g1(y°y0) - e^gm(Yl,y_)
1 2 i i 12 12
(Note that in the case where the functions are separable in
y., and y« as in our model, the convexity of the sub-problems
in the continuous variables is sufficient to ensure condition
(9) and the validity of the decomposition.)
The decomposition procedure then consists of alternately sol-
ving the relaxed sub-problems and the following
Master Problem
Find a feasible solution to,
i) gk(Y2) £ o
(incumbent)
3) g (y-,,Y9) + 2 (-x )g (y-,,y0) <1 -e
" peH p L 2
— 13 —
-------
where a constraint of type 2) is added when an optimal solu-
tion is obtained for the sub-problem and a constraint of type
3) is added when the sub-problem is inconsistent.
Convergence to an e-optimal solution in a finite number of
steps using this decomposition technique is easily demonstrated.
Assume we return from a sub-problem without achieving a bound
ed gain in the incumbent solution when we obtained a feasible
solution then
,.., m,o » . m,c cx
(ii) g (y1Jy2) 2 g (Y1^2)
and
(12) Z(-.x,.) . gi(y°,y2) > -e
if this set of dual variables is a repeat we contradict
,,_\ m,o > , „ , ,i,o \ , m/c c\
(is) g (y1,y2) + s(-xi)g (y-^y^ £ g (y^y^ ~ €
which has been established from the master problem. Since
there are only a finite number of dual solutions between bound-
ed gains in a finite number of steps we exceed any lower bound
or fail to solve the master. When we fail to solve the master
we have established from one of the lower bound functions that
we are within e of the optimum.
When we return from the sub-problems with an inconsistent
solution we have
(14) E (-x ) - gp(y° y ) > -gw(y° y ) - e
peH p *- * ^ z
Again if w and H are a repeat we contradict
(15) gW(yJ,Y2) + Z (-X )gP(y° y ) £ -e
peH F
which has been established from the master. So in at most a
finite number of steps we resolve the inconsistency or demon-
strate that the complete problem is inconsistent.
This type of decomposition has recently been employed in a
large warehouse location problem (see [3]) and proved quite
successful. It would seem to offer great promise in the area
of regional water quality control.
- 14 -
-------
IV. TRANSFER FUNCTIONS
Part 1: Theoretical Considerations
We begin by considering the two fundamental sets of mass bal-
ance differential equations. The most general form of these
equations is given by Thomann [11] or as equations (1) and
(2) by Pence, Jeglic and Thomann [8, pp. 383-83] . These basic
D.E.'s characterize a one-dimensional estuary divided into n
sections. The first set gives the transformation of polluter
BOD to in-stream BOD and can be written as
V.dL.
"alT = Qi-l,i[5i.1,iLi.1-Kl-gi.lji)LiI-ffli_lji(Li_1-Li)
where the sections are numbered from upstream to downstream.
LJ_ is the ultimate carbonaceous BOD in section i. QI_I i i-s
the net flow from section i-1 to i. Vi is the volume of sec-
tion i. §i-i i is the advection coefficient from section i-1
to i. E^_;L ±' is the eddy exchange coefficient from section
i-1 to i. d is the decay rate of BOD in all sections. K^ are
all direct sources of BOD discharged to section i.
Since we require dL^/dt = 0 in the steady state (SS) the above
system can be written as follows
where
and r. . = 0 if | i-j
•J
In matrix notation we have
RL + K = 0
where R is a square matrix of order n. Assuming R is non-
singular we have
L = -R"^
The second set of D.E.'s gives the transformation of stream
- 15 -
-------
BOD into D.O. and can be written as
dC.
-dLi)+Pi
where GJ is the D.O. concentration in section i. r^ is the
reaeration rate of D.O. in section i. Cs^ is the saturation
concentration of D.O. in section i. P^ is all other sources
and sinks of D.O. in section i. Since we require dC^/dt = 0
in the steady state the above system can be written as follows
0 =
where
s .
and s.. = 0 if
Let V be a diagonal matrix of order n with the V. as the diagonal
elements
V =
V, 0
0 V,
0
V
n
and let W be the column vector
W =
rics
rnCS
- 16 -
-------
In matrix notation, the above system is
SC = dVL - V W + P
where S is a square matrix of order n. Assuming S is non-
singular we have
C = -dS^V R-1K - S~1V W + S~LP
The essential idea in calculating the SS transfer matrix is
to inject 1 unit of BOD into each section and calculate the
corresponding D.O. changes in each section. A symmetric in-
terpretation is to remove 1 unit of BOD from each section.
This gives a sign reversal in the final result. Replacing
K above, by ej (unit vector) and noting that the 2nd and 3rd
terms on the right can be dropped we define
R~1e .
3
where (f) . . is the jth column of the transfer matrix. The trans-
fer matrix $ is then
5 = -dS'-'-V R"1! = -ds'V R"1 = -d(R V~1S)~1
The transfer matrix corresponding to a 1 unit removal in each
section is
$ = d(R V^S)"1
The mass balance D.E.'s hold for non- tidal streams as well as
estuaries. It is an interesting exercise to show that the
classical Streeter-Phelps equation is contained in these D.E.'s,
The Streeter-Phelps equation can be derived as follows:
In the absence of reaeration assume that the rate of change
of BOD at time t is proportional to the BOD at time t
or
where yt is the demand for oxygen (BOD) at time t and L is the
ultimate BOD at time zero. The general solution to this D.E.
is [1, Thm.l, p. 40]
(t) = ce~klfc
since L = (j)(0) = c, we have
= Le"1^!1^
Similarly in the absence of demand for oxygen assume that the
rate of change of D.O. deficit (DOD) is proportional to the
- 17 -
-------
DOD at time t.
diy
where Dfc is the D.O. deficit at time t, i.e., Dfc - Cs-Ct-
as before
Combining the above D.E.'s we have the Streeter-Phelps D.E.
-r-7— = k y -k0D, or D!+k0D. = k,y,
at 1 t 2 t tzt J_ t
The general solution is [ 1, Thm. 2, p. 41]
^"[^ «-» 4- \f rt \r •—T^" «""* "I"
/ _i_ \ JV / >— /* JX ,' ^.T T , JV / •—
(t) = e z / e * k,y dx + ce ^
n
0
using y = Le we have
1
, -
(t) = \\ f e2-l (v -k.)dx
k2~kl 0 2
klL r -kit -k2t, ^ ^ -
= Z t® ~ e + 'de
from the initial value problem we have D_ = () (0) = c and the
classic Streeter-Phelps equation is
In the Streeter-Phelps D.E.Q. let D, = C -C. giving
t w t
d(Cs-Ct)
dt klYt~ k2 (CS~Ct)
or
(16) -XT- = k2CS~ k2Ct~ klYt
vj.1— ^O ^L- J_l_
In the mass balance D.E. giving the transformation of stream
- 18 -
-------
BOD to D.O. let E,.= 0 and £. , .=1 which gives
i j • i— J.} i
dCi Qi-l iCi-l Qi i+ici
H?^ - = — ' - 3- , 3-TX 3. , _ _ -,T
(17) dt v± V.. + i Si i i " i
Note that in (16) and (17) the following quantities are equi-
valent: k2~r, k-j^-d, and yt~L. Comparing (16) and (17) we can
see that the terms of (16) are contained in (17).
In the D.E. giving the transformation of stream BOD into D.O.
we introducted the term P^ representing all other sources
and sinks of D.O. in section i. Now let P^ represent sources
of D.O. and Q^ sinks of D.O. in section i. From previous re-
sults we have
C = $K - S V W - S"1? + S-1Q
where P and Q are column vectors. The same basic idea used
to construct a SS transfer matrix for BOD can be used to con-
struct a SS transfer matrix for D.O. sources. We inject 1
unit of D.O. into each section and calculate the corresponding
D.O. changes in each section. Define
„-!
where ty-j is the jth column of the transfer matrix. The SS
transfer matrix Y is then
Y = -S"1! = -S"1
This SS transfer matrix can be used to account for D.O. changes
due to the following:
1) D.O. in the polluter discharges
2) D.O. in the regional plant discharges
3) in-stream aeration
4) flow augmentation.
For 1) and 2) let QJ-, and R, be the average D.O. concentration
of effluent from polluter i and regional plant i respectively.
Thus 1) and 2) are handled by the following modification to
the quality constraints
kt..(l-r.)p J kt..(l-r.)p J
AFJ + A 2 =H= ~ + C - AF J - A 2
- 19 -
-------
In-stream aeration 3) requires the introduction of a variable
specifying the amount of D.O. added in each section due to
aeration. Flow augmentation requires a conservation of flow
constraint for each reservoir and additional flow variables.
Part 2: Numerical Calculation of a Transfer Matrix
A general nonlinear programming problem for water quality con-
trol in a tidal estuary has been given by Graves., Hatfield
and Whins ton [4]. The flow chart [4, page 4] of the type of
operations involved in the overall approach is given below:
Basic Data from
Continuous
Monitoring
Basic Data
(costs,, removals,
etc.)
Compute Transfer
Matrix
Solve Nonlinear
Programming
Problem
Adjust Goals,,
Change Plant
Locations,, etc.
Acceptable
Solution
FIGURE 2
OPTIMIZATION PROCESS,, ASSUMING FIXED TRANSFER MATRIX
- 20 -
-------
In [4] the SS (Steady State) transfer coefficients were taken
as a matrix of fixed constants and the path labeled 1 -» 2 - 3
of Figure 2 was not investigated. There was no need to re-
compute the SS transfer coefficients since they were assumed
constant. Unfortuantely the SS transfer matrix A, is not com-
pletely independent of the variables of the nonlinear program-
ming problem. The calculation of the SS transfer matrix de-
pends, among other things, on the estuarial intersection flows.
Activating by-pass pipes and switching plant discharge points
can change the estuarial flow pattern. This suggests some type
of iterative process of: 1) solving the nonlinear programming
problem; 2) using the current solution to update the transfer
matrix; and 3) re-solving the nonlinear programming problem
using the updated transfer matrix. With these ideas, we can
modify the path 1-2 - 3 of Figure 2 as in Figure 3.
The purpose of this section is to develop the details of the
calculations implied by Figure 3. We begin by considering the
calculation of the transfer matrix. Let Qi-j^i be the net
flow from section i-1 to i. V± is the volume of section i.
F-_1 • is the advection coefficient from section i-1 to i.
Ei-l i is the eddy exchange coefficient from section i-1 to i.
d is*the decay rate of BOD in all sections. ri is the reaera-
tion rate of D.O. in section i. For computational purposes,
a convenient system of units is Qi-i i and Ei_]^i in km^/day,
V- in km3 and d and r± in I/day. Note that ?j_-i}i_ is dimen-
sionless. L
Compute Transfer
Matrix
Solve Nonlinear
Programming
Problem
Update
Transfer Matrix
Adjust Goals,
Change Plant
Locations, etc,
FIGURE 3
ITERATIVE PROCEDURE FOR FINDING STABLE TRANSFER MATRIX
- 21 -
-------
Let R and S be square matrices of order n (n = the number of
sections). If reR and seS and
then
r. . = c. . - V.d
11 11 i
s . . = c . . - V . r .
-
and rij = sij = ° if 'i - J I 2 2
Let V be a diagonal matrix of order n with the V^ as the dia-
gonal elements. Then the SS transfer matrix corresponding to
a 1 Ib/day BOD input into each section is given by
(20) A = -(4.536 • 10"7) d (R V~1S)~1
With the system of units chosen as above the elements of R V S
will be close to one and accuracy will be maintained in cal-
culating the inverse. The units of A are in mg/L/lb/day.
In order to calculate a SS transfer matrix for the Delaware
Estuary j we need numerical values for the Oi_i i, 5i-i i.>
Ei_i<;j_, Vj_j ri and d. These are in turn functions of the
physical characteristics of the estuary. The necessary data
on the physical properties of the estuary is given by Pence ,
Jeglic,, and Thomann [8, page 390] . The values of QI_I 1= Q(i)
in km3/day for a base flow of 3000 cfs are given in Ta6le 1.
Values of E-J__J_ i in km^/day are calculated from (see [8] equa-
tion 6)
(21) E = 1.57873 ^
where KI j_+± is the eddy-diffusion coefficient from section i
to i+1. 'AI^I+I is the cross-sectional area between section i
and i+1; and i^_ is the length of section i. Values of
Ki,,i+lj Ai,,i+l and "i are given in [8, page 390]. The EI_I 1=
E(i) are tabulated in Table 2. The advection coefficient is
taken to be PI- 1,1= ^(i) = 0.5 = i=5 , 6, . . . , 31. In sections one
to four^ the advective forces are stronger than the tidal action
22 -
-------
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
Q (i) =Q±_lfi
7.34055E-03
7.43108E-03
7.67577E-03
7.69779E-03
7.77119E-03
7.86907E-03
8.35844E-03
7.69779E-03
7.74672E-03
7.77119E-03
8.29727E-03
8.37067E-03
8.37067E-03
8.48323E-03
8.95057E-03
9.92932E-03
;
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
1.04040E-02
1.07123E-02
1.08591E-02
1.09081E-02
1.09081E-02
1.13778E-02
1.13803E-02
1.13925E-02
1.13925E-02
1.14047E-02
1.14170E-02
1.14170E-02
1.15393E-02
1.16127E-02
1.16861E-02
TABLE 1
NET FLOWS INTO SECTION i
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
E (i) = Ei_1 ±
9.77309E-04
2.43355E-03
3.37848E-03
3.88367E-03
4.49938E-03
5.38346E-03
8.71458E-03
1.56610E-02
2.02077E-02
2.18654E-02
2.40361E-02
2.56543E-02
2.60490E-02
2.86934E-02
2.32336E-02
1.93394E-02
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
2.07010E-02
2.23784E-02
2.49241E-02
2.81803E-02
3.95734E-02
5.99522E-02
6.42937E-02
6.97798E-02
1.00399E-01
1.04101E-01
1.08521E-01
1.18523E-01
8.65670E-02
7.03126E-02
8.49277E-02
TABLE 2
EDDY EXCHANGE COEFFICIENTS BETWEEN SECTIONS i-1 AND i
- 23 -
-------
and we take F ( 1 ) = 0.86686 5 ( 2 ) = 0.76969 F ( 3 ) = 0.61619,
and p(4 ) = 0.52767. Values of Vj_ in km3 are taken from [8]
pate 390 and are recorded in Table 3.
The reaeration rate r-j_ and the decay rate d are functions of
water temperature. This has the effect of introducing a de-
gree of freedom in the calculation of the SS transfer matrix
since each matrix will depend on the particular value of water
temperature chosen. Values of r-j_ in I/day are calculated from
(see [8] equations 8 and 9)
24 (8.1 • IP"5 u±)°'5 .018 (t _ 20)
ri = - Ti - e
where u^ is the mean tidal velocity in section i and h^ is the
mean depth in section i. Values of u^ and h^ are given in [8]
page 390. A value of t=20° C=68° F was chosen because:
(a) it simplifies the calculation of (22) ; and (b) it is
roughly the average water temperature from June to November.
The r-j_ are tabulated in Table 4. Values of d in I/day are
calculated from (see [8] equation 10)
d = 0.23 • 1.099(t "20)
With t=20° C we have d = 0.23/day-
The SS transfer matrix for the Delaware Estuary calculated
from (20) is given in Table 5. The SS transfer matrix supp-
lied to us by Prof. R.V. Thomann of Manhattan College , and
used in our earlier work is given in Table 6 for comparison.
Updating the Transfer Matrix:
Recall that there are five sets of variables (the F, P, T,
J and R matrices) in the nonlinear programming formulation.
It is relatively simple calculation to update the intersec-
tional flows Q^ i+i* given the current P and T matrices and
the original polluter discharges. If the base flow in the
stream is large, changing the polluter discharge points effects,
primarily, the intersectional flows. For the Delaware Estuary
problem we are assuming that all other secondary effects are
negligible. Some of the other quantities indirectly effected
by changing the flow pattern are: 1) the advection coeffi-
cient li^i+i/ and 2) the average depth h^ and consequently
the reaeration rate r^. This assumption is apparently justi-
fied for the Delaware Estuary; however, this is not always the
case, especially in situations where polluter discharges make
up a large percentage of the base flow. A detailed analysis
of this particular aspect of the problem is given in Section V.
- 24 -
-------
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
V.
i
6.85344E-03
1.03085E-02
1.30272E-02
1.50662E-02
1.80115E-02
2.14099E-02
1.28856E-02
1.42733E-02
1.50945E-02
1.64822E-02
1.78416E-02
1.85496E-02
1.96541E-02
2.27976E-02
5.26752E-02
5.74896E-02
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
6.18509E-02
6.78547E-02
7.62374E-02
8.30342E-02
4.28198E-02
4.45757E-02
4.80874E-02
5.07494E-02
5.23920E-02
5.44877E-02
5.81693E-02
6.36634E-02
1.38655E-01
1.59158E-01
TABLE 3
VOLUMES OF SECTION i
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
r .
i
1.00116E-01
1.19462E-01
9.23170E-02 .
1.00116E-01
1.15604E-01
1.25953E-01
1.06598E-01
9.87044E-02
1.15604E-01
8.55555E-02
8.94839E-02
5.71866E-02
4.56965E-02
8.39498E-02
8.94839E-02
1.12356E-01
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
1.30556E-01
1.12356E-01
1.12356E-01
1.12356E-01
1.21359E-01
9.93308E-02
1.21359E-01
1.21359E-01
1.41016E-01
1.95081E-01
1.78124E-01
1.78124E-01
1.29738E-01
1.20993E-01
TABLE 4
REAERATION RATE IN SECTION i
- 25 -
-------
I
to
TABLE 5
TRANSFER MATRIX
1
2
3
4
5
6
7
8
9
10
li
12
13
14
15
16
17
16
19
20
21
22
23
24
25
26
27
28
29
30
1
-1. 0180-05
-1.778D-05
-2.380D-05
-2.594D-C5
-2.4C1D-C5
- .9340-05
- .7210-05
- .4940-05
- .3C2D-05
- .1360-05
-1.CC4C-05
-8.9C60-C6
-7.633C-06
-6.267C-06
-4.581C-06
-3.033C-06
-1.9450-06
-1.262C-06
-8.0370-07
-5.0710-07
-3.5970-07
-2.7680-07
-2.0600-07
-1.4860-07
-1.1310-07
-8.3920-08
-6.213D-08
-4.6210-08
-2.8770-08
-1.245C-08
2
-1.5880-06
-1.2440-05
-2.0900-05
-2.4890-05
-2.4120-05
-1.9980-05
-1.8C20-05
-1. 5770-05
-1.381D-05
-1.2110-05
-1.0750-05
-9.5610-06
-8.2150-06
-6.7590-06
-4.9540-06
-3.2890-06
-2. 1130-06
-1.3730-06
-8.7580-07
-5.5310-07
-3.9250-07
-3.0220-07
-2.2490-07
-1.6230-07
-1.2360-07
-9.1730-08
-6.7930-08
-5. 0530-08
-3.1470-08
-1.3620-08
3
-1.113C-07
-1.0960-06
-1.4370-05
-2.2060-05
-2.3690-05
-2.0750-05
-1.9190-05
-1.704D-05
-1.509C-05
-1.333C-05
-1.1910-05
-1.0650-05
-9.1940-06
-7.5910-06
-5.5910-06
-3.7260-06
-2.4020-06
-1.566C-06
-1.0010-06
-6.33CO-07
-4.4970-07
-3. 4640-07
-2.5800-07
-1.8630-07
-1.4190-07
-1.054C-07
-7. 8050-08
-5. 8070-08
-3.6190-08
-1. 5670-08
4
-3.7840-09
-4.0820-08
-6.9190-07
-1.5C6D-05
-2.1270-05
-2.0870-05
-2.0200-05
-1.8380-05
-1.6550-05
-1. 4820-05
-1.3380-05
-1.2060-05
-1.0480-05
-8.7010-06
-6.4540-06
-4.3280-06
-2.8C40-06
-1.8350-06
-1.1760-06
-7.4610-07
-5.3070-07
-4.0920-07
-3.0500-07
-2.2040-07
-1. 6790-07
-1.2480-07
-9.2490-08
-6.8850-08
-4.2930-08
-1.86CD-08
5
-2.6160-10
-2.9550-09
-5.5450-08
-1.5780-06
-1.5420-05
-1. 9510-05
-2. 0470-05
-1. 9380-05
-1.7900-05
- .6350-C5
- .4990-05
- .3660-05
- .1980-05
- .0020-05
-7.5CCD-06
-5.0710-06
-3.3C60-06
-2.1740-06
-1.4CCD-06
-8.9050-07
-6.344D-C7
-4.8980-07
-3.655D-C7
-2.6430-07
-2.0160-07
-1.499D-C7
-1. 1120-07
-8.2790-08
-5.1670-08
-2.240C-C8
6
-3.347C-11
-3.8S4C-IC
-7.686C-09
-2.4350-07
-3. 0530-06
-1.52CC-05
-1.9UC-05
-1.949C-05
-i.eecc-05
- .773C-05
- .6630-05
- .541C-C5
- .3690-05
- .157C-05
-8.776C-06
-5.999C-C6
-3.944C-06
-2.6120-06
-1.69CC-C6
-1.C8CC-06
-7.7C9C-07
-5.96CC-C7
-4.453C-07
-3.225C-07
-2.462C-07
-1.6320-07
-1.36CC-07
-1.C14C-07
-6.3320-08
-2.747C-C8
7
-1.0350-11
-1.216D-1C
-2.465C-C9
-8.153C-C8
-1.107D-C6
-6.618D-Ct
- .496C-C5
- .77CD-05
- .838C-C5
- .819C-C5
- .7630-C5
- .672C-C5
- .5110-05
-1.294C-C5
-9.976C-C6
-6.911D-C6
-4.59CC-C6
-3.0630-06
-1. 995C-C6
-1.28CC-06
-9.163C-C7
-7.096C-C7
-5.31CO-C7
-3.85CO-07
-2.942C-C7
-2.1920-C7
-1.628C-C7
-1.215C-07
-7.598C-C6
-3.299C-C6
8
-5.3750-12
-6.367C-U
-1.3C9C-09
-4.4350-08
-6.2670-07
-4.C470-06
- .C650-C5
- .5360-05
- .7240-05
- .7850-05
- .7780-C5
- .718C-C5
- .574C-05
- .3620-05
-1.C63C-05
-7.4350-06
-4.973D-C6
-3. 3380-06
-2.1830-06
-1.4050-06
-1.CC8C-06
-7.8HD-C7
-5.8510-07
-4.2470-C7
-3.2470-07
-2.42CO-07
-1.7990-C7
-1.3430-07
-e.4CRD-OR
-3.652C-C8
9
-3.14SC-12
-3.752C-11
-7.791C-1C
-2.684C-C8
-3.896C-C7
-2.635C-C6
-7.484C-C6
- .167C-C5
- .516C-C5
- ,68tC-C=
- .75CC-C5
- .735C-C5
- .618C-C5
- .419C-C5
- .124C-C5
-7.961C-C6
-5.371C-C6
-3.629C-C6
-2.385C-C6
-1.541C-Ot
-1.1C7C-C6
-g.595C-C7
-6.445C-C7
-4.683C-C7
-3.584C-C7
-2.673C-C7
-1.969C-C7
-1.485C-C7
-9.3C7C-CP
-4.046C-C8
1C
-1.872C-12
-2.241011
-4.693C-1C
-1.64CC-C8
-2.43CC-C7
-1.7CIC-C6
-5.CE6C-C6
-6.3C6C-C6
-1.16CC-C5
-1.492C-C5
-1.661C-C5
-1.714C-C5
-1.642C-C5
-1.467C-C5
-i.ieec-c?
-6.544C-Ct
-5.83CC-C6
-3.972C-C6
-2.627C-C6
-1.7C6C-C6
-1.228C-C6
-9.5E1OC7
-7.173C-C7
-5.219C-C7
-3.998C-C7
-2.985C-C7
-2.222C-C7
-1.661C-C7
-1.C42C-C7
-4.533C-C8
-------
TABLE 5 Cont.
TRANSFER MATRIX
1
to
^]
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
11
-1.1390-12
-1. 3690-11
-2.8880-10
-1.0200-08
-1.5380-07
-1.1C50-06
-3.4240-06
-5.7590-06
-8.3990-06
-1.1600-05
-1.4780-05
-1.6290-05
-1.6260-05
-1.4910-05
-1.2430-05
-9.1290-06
-6.3180-06
-4.3500-06
-2.8990-06
-1.8930-06
-1.3670-06
-1.C65D-06
-8.0120-07
-5.8380-07
-4. 4770-07
-3.346C-C7
-2.4940-07
-1.8660-07
-1.1720-07
-5.1C40-08
12
-7.0230-13
-8.4710-12
-1.7980-10
-6.4110-09
-9. 7900-08
-7.1740-07
-2.2820-06
-3.9180-06
-5.8740-06
-8.4670-06
-1.1540-05
-1.4470-05
-1.5450-05
-1.4770-05
-1.2830-05
-9.6870-06
-6. 8290-06
-4.7640-06
-3.2050-06
-2.1070-06
-1.5270-06
-1.1920-06
-8.9880-07
-6.5620-07
-5.0400-07
-3. 7710-07
-2.8140-07
-2.1C8D-07
-1.3260-07
-5.7800-08
13
-4.2450-13
-5.1360-12
-1.0960-10
-3.939C-09
-6.0840-08
-4.5310-07
-1.4710-06
-2.56CC-06
-3.9250-06
-5.8250-06
-8.2840-06
-1.113C-05
-1.3630-05
-1.3990-05
-1.2960-05
-1.0180-05
-7.3590-06
-5.2220-06
-3.5550-06
-2.3570-06
-1.7150-06
-1.3430-06
-1.0150-06
-7.4290-07
-5.7140-07
-4.2830-07
-3.2010-07
-2.40CC-07
-1. 5120-07
-6.6030-08
14
-2.5900-13
-3.1420-12
-6.7320-11
-2.4360-09
-3. 7980-08
-2.866D-C7
-9.4580-07
-1.6680-06
-2.5900-06
-3. 9220-06
-5.7360-06
-8.0310-06
-1.0550-05
-1.2380-05
-1.2680-05
-1.0520-05
-7.8490-06
-5.6870-06
-3.9270-06
-2.6290-06
-1.9220-06
-1.5100-06
-1.1450-06
-8.3970-07
-6.4710-07
-4.8580-07
-3.6360-07
-2.7300-07
-1.7240-07
-7.5370-08
15
-1.2950-13
-1.5770-12
-3.3980-11
-1.2400-09
-1.9560-08
-1.5COD-07
-5.0460-07
-9. 0210-07
-1.4240-06
-2.2040-06
-3.3170-06
-4.8340-06
-6.7510-06
-8.7280-06
-1.1570-05
-1.0710-05
-8.4430-06
-6.3280-06
-4.4670-C6
-3.0350-06
-2.2350-06
-1.7640-06
-1.343C-06
-9.8860-07
-7.6380-07
-5.7480-07
-4.3130-07
-3.2450-07
-2.0550-07
-9.0010-08
16
-5.2980-14
-6.468C-13
-1.4CCO-11
-5.14SO-1C
-8.197C-09
-6.3670-08
-2.172C-07
-3.923C-07
-6.2650-07
-9.849C-07
-1.511C-06
-2. 2590-06
-3.268C-C6
-4.442C-C6
-6.529C-06
-9.7590-06
-9.C68C-06
-7.3S9C-06
-5.49CC-06
-3.851C-06
-2.878C-06
-2.2940-06
-1.760C-06
-1.3C50-06
-1.014C-06
-7.666C-07
-5.777C-C7
-4.3630-07
-2.7780-07
-1.222C-07
17
-2.1420-14
-2.6200-13
-5.693C-12
-2.104D-1C
-3.372C-CS
-2.642C-08
-9.1050-C8
-1.6550-C7
-2.664D-C7
-4.2290-07
-6.568D-C7
-9.9650-C7
-1.471C-C6
-2.054C-C6
-3.173C-C6
-5.5390-C6
-8.592D-C6
-8.291C-C6
-6.676D-C6
-A.906C-C6
-3.743C-06
-3.022C-C6
-2.343C-C6
-1.755D-C6
-1.372D-C6
-1.0440-06
-7.90SO-C7
-6.002C-C7
-3.8460-07
-1.7COO-07
18
-8. 6760-15
-1. 0630-13
-2.316D-12
-8. 5920-11
-1.3840-09
-1. 0920-08
-3.7890-08
-6.S210-06
-1.1200-07
-1.791D-C7
-2.8040-07
-4. 2980-07
-6.4320-07
-9.1310-07
-I. 4520-06
-2.7440-06
-5.C26D-06
-8.1260-06
-7.6740-06
-6.C85D-C6
-4.7870-06
-3.9370-06
-3.1COC-06
-2.3520-06
-1.R550-06
-1.4220-06
-1.C850-06
-8. 2070-07
-5.3550-07
-2.3810-07
IS
-3.5S9C-15
-4.415C-14
-•3.639C-13
-3.586C-11
-5.8CCC-1C
-4.598C-C9
-1.6C4C-C8
-2.941C-C8
-4.779C-C8
-7.t7SC-Ce
-1.21CC-C7
-1.868C-C7
-2.621C-C7
-4.C5CC-C7
-6.5ttC-C7
-1.3CCC-C6
-2.586C-C6
-4.916C-C6
-7.596C-C6
-7.C15C-C6
-5.817C-C6
-4.93CC-C6
-3.S72C-C6
-3.073C-C6
-2.454C-C6
-1.9C3C-C6
-1.467C-C6
-1.129C-C6
-7.381C-C7
-3.3C8C-C7
2C
-1.554C-15
-1.9C8C-14
-4.172C-13
-1.5560-11
-2.524C-1C
-2.CC8C-C9
-7.C36C-C9
-1.293C-C8
-2.1C8C-C8
-3.359C-C8
-5.376C-C8
-8.344C-C8
-1.268C-C7
-1.836C-C7
-3.C15C-C7
-6.152C-C7
-1.283C-C6
-2.635C-C6
-4.744C-C6
-6.S6SC-C6
-6.4690-Ct
-5.7S3C-C6
-4.857C-C6
-3.876C-C6
-3.158C-C6
-2.49CC-C6
-1.948C-C6
-i.siec-ce
-1.CC8C-C6
-4.566C-C7
-------
to
00
TABLE 5 Cont.
TRANSFER MATRIX
1
2
3
4
5
6
7
8
9
10
11
12
13
14
IS
16
17
18
19
20
21
22
23
24
25
26
27
26
29
30
21
-8.4490-16
-1.0380-14
-2.272D-13
-8.4820-12
-1.3780-10
-1.C98C-09
-3.857D-09
-7.C970-C9
-1.1590-08
-1.872D-08
-2. 9670-08
-4.6170-08
-7.0410-08
-1.C23C-07
-1.6SOO-C7
-3.4S8D-C7
-7.447C-07
-1.5780-06
-2. 9940-06
-4.9220-06
-6.CC5C-06
-5.974C-06
-5.3540-06
-4.4820-06
-3.757C-C6
-3.C31C-06
-2.418C-06
-1.9130-06
-1.2950-06
-5.945C-07
22
-5. 3980-16
-6.6350-15
-1.4530-13
-5.4280-12
-8.8270-11
-7.0460-10
-2.4770-09
-4.5630-09
-7.4570-09
-1.2C6D-08
-1.915D-08
-2.985D-08
-4.5620-08
-6.6440-08
-1. 1020-07
-2.3C2D-07
-4.9650-07
-1.0720-06
-2. 0950-06
-3.6370-06
-4.9300-06
-5.6350-06
-5.4300-06
-4.7600-06
-4.0940-06
-3.3690-06
-2.7310-06
-2.1870-06
-1.5C30-06
-6.9700-07
-3
-4
-9
-3
-5
-4
-I
-2
-4
-7
-1
-1
-2
-4
-7
-1
-3
-7
-1
-2
-3
-4
-5
-4
-4
-3
-3
-2
-1
-8
23
.4300-16
.2180-15
.24CC-14
.455C-12
.624C-11
.4940-10
.5820-09
.9160-09
.7700-09
.7240-09
.2280-08
.9170-08
.934C-08
.2840-08
.131C-08
.5CCC-07
.271C-07
.1720-07
.4340-06
.5880-06
.745C-06
.6000-06
.0990-06
.8230-06
.308C-06
.6460-06
.0200-06
.4580-06
.7220-06
.083C-07
24
-2.1620-16
-2.6590-15
-5.8280-14
-2.1810-12
-3.5520-11
-2.8420-10
-I. 0010-09
-1.8470-09
-3.0240-09
-4.9010-09
-7. 8000-09
-1.2190-08
-1.8700-08
-2.7340-08
-4.5650-08
-9.6690-08
-2.1280-07
-4.7230-07
-9.6160-07
-1.7880-06
-2. 7080-06
-3.4780-06
-4.1510-06
-4.5430-06
-4.3210-06
-3.8140-06
-3.2570-06
-2.7090-06
-1.9470-06
-9.2770-07
25
-1.5220-16
-1.8720-15
-4.1C5C-14
-1.5370-12
-2.5040-11
-2.CC50-10
-7.0710-10
-1.305C-09
-2.137C-09
-3.467D-C9
-5. 5210-09
-8.6370-C9
-1.3260-08
-1.9420-08
-3.2490-08
-6.91CD-08
-1.5300-C7
-3.4230-07
-7.0500-07
-1.335D-C6
-2.C75C-06
-2.7280-06
-3.3750-06
-3.9250-06
-4.0830-06
-3. 8020-06
-3.3650-06
-2.8670-06
-2. 1150-06
-1.0240-06
26
-I.C62C-16
-1.3C6C-15
-2.8t5C-14
-1.073C-12
-1.749C-11
-1.4C1C-10
-4.946C-1C
-9.131C-1C
-1.496C-09
-2.428C-09
-3.869C-09
-6.05EC-C9
-9.3C9C-C9
-1.365C-08
-2.2E7C-08
-4.ee3C-08
-1.C87C-07
-2.448C-07
-5.089C-07
-9.772C-C7
-1.55CU-06
-2.C73C-06
-2.63CC-06
-3.180C-06
-3.483C-C6
-3.584C-C6
-3.367C-06
-2.975C-06
-2.28CC-C6
-1.127D-06
27
C.O
-9. 0050-16
-1.975C-14
-7.399C-13
-1.2070-11
-9.675C-11
-3.416C-1C
-6.31CC-IC
-1.034C-C9
-1.679C-C9
-2.678C-C9
-4.1950-C9
-6.452C-CS
-9.468C-09
-1.5890-C8
-3.404C-C8
-7. 6070-08
-1.724C-C7
-3.612C-C7
-7.017D-C7
-1.131C-C6
-1.533C-C6
-1.9810-Ct
-2.463C-C6
-2.7880-C6
-3.041C-C6
-3.1840-C6
-2.984U-C6
-2.418C-C6
-1.23CD-C6
28
C.C
-6.2130-16
-1.3630-14
-5.1080-13
-6. 3350-12
-6.6840-11
-2.3610-10
-4.3630-1C
-7.1540-10
-1.1620-09
-1.854C-C9
-2.-5C60-09
-4. 4720-09
-6. 5680-09
-1.1C4D-08
-2.371C-C8
-5. 3170-08
-1.2100-07
-2.552C-07
-5.CC50-07
-8.165C-C7
-1.1180-06
-1.4650-C6
-1.8560-06
-2.149U-06
-2.4290-C6
-2.6S8D-06
-2.B110-06
-2.4800-06
-1.3120-06
29
C.C
-3.3K8C-16
-7.434C-15
-2.787C-13
-4.55CC-12
-3.651011
-1.291C-1C
-2.386C-1C
-3.914C-1C
-6.361C-1C
-1.C15C-CS
-1.593C-C9
-2.454C-C9
-3.6C8C-C9
-6.075C-C9
-1.31CC-C8
-2.952C-C8
-6.762C-C8
-1.43RC-C7
-2.856C-C7
-4.7??C-C7
-6.566C-C7
-8.745C-C7
-1.134C-Ct
-1.347C-C6
-1.581C-C6
-1.86CC-C6
-2.116C-C6
-2.372C-Cf
-1.373C-C6
3C
C.C
-1.234C-16
-2.7C9C-ie
-1.C16C-13
-1.659C-12
-1.332C-U
-4.712C-11
-8.712C-11
-1.43CC-1C
-2.325C-1C
-3.713C-1C
-5.829C-1C
-8.986C-1C
-1.322C-C9
-2.22SC-C9
-4.819C-C9
-i.cscc-ce
-2.51CC-C8
-5.37ic-C6
-1.C76C-C7
-1.8C3C-C7
-2.523C-C7
-3.398t-C7
-4.471C-C7
-5.3S3C-C7
-6.476C-C7
-7.859C-C7
-9.333C-C7
-1.148C-C4
-1.1C3C-C6
-------
TABLE 6
THOMANN'S TRANSFER MATRIX
to
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
2^
24
25
26
27
28
29
30
1
7.834E-06
1.3872-05
1.7712-05
1.806E-05
1.616E-05
1.159E-05
9.7612-06
8.196E-06
6.6332-06
5.3832-06
4.5712-06
3.o$3E-u6
3.3S42-06
2.779E-06
2.016E-06
1.382E-06
9.3522-07
6.2032-07
4.028E-07
2.577E-07
1.835E-0?
1.3852-07
9.9232-08
6.T71E-08
4.7502-08
3.218E-08
2.1872-08
1.5122-08
8.650E-09
3.7662-09
2
9.331E-07
1.0692-05
1.706S-05
1.8922-05
1.7682-05
1.3022-05
i.noE-05
9.3892-06
7.643E-06
6.2342-06
5.311E-06
4.6502-06
3.951E-06
3.2502-06
2.3622-06
1.6212-06
1.0982-06
7.2912-07
4.7362-07
3.0312-07
2.1592-07
1.6292-07
1.1673-07
7.9672-08
5.5892-08
3.786E-08
2.574E-08
1.7802-08
1.0192-08
4.4322-09
3
5.7672-08
8.3562-07
1.2402-05
1.7702-05
1.8342-05
1.428E-05
1.2482-05
1.0702-05
8.8082-06
7.2512-06
6.2172-06
5.4692-06
4.6642-06
3.848E-06
2.8062-06
1.9312-06
1.311E-06
8.7102-07
5.6632-07
3.6252-07
2.5832-07
1.9502-07
1.397E-07
9.537E-08
6.6912-08
4.533E-08
3.081E-68
2.131E-08
1.2202-08
5.308E-09
4
1.3432-09
2.1302-08
4.084E-07
1.2742-05
1.7502-05
1.528E-05
1.399E-05
1.2292-05
1.0302-06
8.6122-06
7.4592-06
6.611E-06
5.6722-06
4.701E-06
3.447E-o6
2.3812-06
1.62UE-06
1.079E-06
7.0232-07
4.4992-07
3.2072-07
2.4222-07
1.7362-07
1.1852-07
8.313E-08
5.633E-08
3.8292-08
2.6482-08
1.5172-08
6.599E-09
5
3.271E-11
5.439E-10
1.1572-08
4.770E-07
1.330E-05
1.5292-05
1.5232-05
1.393E-05
1.2032-05
1.0292-05
9.0472-06
8.1082-06
7.0152-06
5.850E-06
4.3212-06
3.0022-06
2.051E-06
1.3682-06
8.9212-07
5.722E-07
4.0812-07
3.082E-07
2.2092-07
1.5082-07
1.059E-07
7.1742-08
4.8782-08
3.3742-08
1.933E-08
8.4102-09
6
9.865E-1.3
1.685E-11
3.7782-10
1.7382-08
6.325E-07
1.2272-05
1.5192-05
1.510E-05
1.378E-06
1.2272-05
1.1052-05
1.0082-06
8.8332-06
7.4372-06
5.555E-06
3.8912-06
2.672E-06
1.7892-06
1.169E-06
7.509E-07
5.359E-07
4.049E-07
2.904E-07
1.983E-07
1.392E-07
9.436E-08
6.417E-08
4.439E-08
2.5442-08
1.107E-08
7
2.0572-13
3.557E-12
8.166E-11
3.921JS-09
1.548E-07
3.7292-06
1.199E-05
1.4352-05
1.4462-05
1.371E-05
1.2802-05
1.1962-05
1.0662-05
9.0912-06
6.8872-06
4.873E-06
3.368E-06
2.265E-06
1.484V -06
9.549E-07
6.8202-07
5.1562-07
3.6992-07
2.5272-07
1.7742-07
1.2032-07
8.1842-08
5.663E-08
3.246E-08
1.413E-08
8
8.947E-14
1.559E-12
3.6262-11
1.782E-09
7.3262-03
1.924E-06
7.4892-06
1.221E-05
1.3832-05
1.3962-05
1.3462-05
1.283E-05
1.161E-06
9.9942-06
7.6552-06
5.4582-06
3.7902-06
2.5572-06
1.679E-06
1.0822-06
7.731E-07
5.8462-07
4.196E-07
2.8672-07
2.014E-07
1.3662-07
9.2912-08
6.4302-08
3.6872-08
1.606E-08
9
3.855E-14
6. 758E-13
1.589E-n
7-953E-10
3.3692-08
9.3722-07
4.044E-06
7.4092-06
1.177E-05
1.3562-05
1.386E-05
1.369E-06
1.2672-05
1.108E-05
8.630E-06
6.224E-06
4.3532-06
2.9502-06
1.9432-06
1.254E-06
8.970E-07
6.788E-07
4.873E-07
3.3322-07
2.341E-07
1.5882-07
1.0802-07
7.4792-08
4.290E-08
1.8692-08
10
1.687E-ll»
2.973E-13
7.0572-12
3.5872-10
1.557E-08
4.524E-07
2.089E-06
4.08UE-06
7.418E-06
1.1912-05
1.356E-06
1.416E-05
1.3562-05
1.213E-05
9.6722-06
7.0832-06
5.0002-06
3.4082-06
2.2532-06
1.4582-06
1.044E»06
7.9052-07
5.679E-07
3.8852-07
2.730E-07
1.852E-07
1.2612-07
8.7322-08
5.011E-08
2.184E-08
-------
TABLE 6 Cont.
THOMANN'S TRANSFER MATRIX
CO
o
1
2
3
4
5
6
7
8
9
10
n
12
13
Ik
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
n
8.9162-15
1.5782-13
3.7712-12
1.938E-10
8.555E-09
2.55TE-OT>
1.229E-06
2.4862-06
4.8102-06
8.6342-06
1.213E-06
1.3822-05
1.3892-05
1.278E-05
1.049E-05
7.8212-06
5.5602-06
3.6292-06
2.5422-06
1.6492-06
1.1822-06
8.9602-07
6.4402-07
4.4082-07
3.0992-07
2.1032-07
1. If 322-07
9.9232-08
5.6972-08
2.484E-08
12
4.8462-15
8.6052-14
2.0692-12
1.073E-10
4.8002-09
1.4662-07
7.2552-07
1.5022-06
3.0232-06
5.7672-06
8.8722-06
1.2352-05
1.3562-05
1.3112-05
1.1242-05
8.6032-06
6.2312-06
4.315E-06
2.883JE-06
1.8772-06
1.3472-06
1.0222-06
7.3522-07
5.036E-07
3.5422-07
2.4062-07
1.6392-07
1.1362-07
6.5262-03
2.847E-08
13
2.5982-15
4.6272-14
1.1182-12
5.8432-11
2.6442-09
8.2202-08
4.163J3-07
8.767^-07
1.8142-06
3.6002-06
5.8342-06
8.9092-06
1.2042-05
1.2742-05
1.1752-05
9.361E-06
6.9282-06
4.8592-06
3.271E-06
2.1412-06
1.5402-06
1.1702-06
8.4252-07
5.7762-07
4.0662-07
2.7632-07
1.6042-07
1.3062-07
7.5122-08
3.2792-08
14
1.4232-15
2.5402-14
6.1612-13
3.2412-11
1.4802-09
.4.6662-08
2.4032-07
5.1282-07
1.0822-06
2.2052-06
3.6902-06
5.9302-06
8.7822-06
1.1302-05
1.1802-05
9.9652-06
7.5932-06
5.4152-06
3.683JE-06
2.4242-06
1.7492-06
1.3312-06
9.5972-07
6.5882-07
4.641E-07
3.1572-07
2.1542-07
1.4952-07
8.6062-08
3.7602-08
15
6.1222-16
1.0962-14
2.6742-13
1.4182-n
6.5502-10
2.10XE-08
1.1042-07
2.3902-07
5.1522-07
1.0802-06
1.8652-06
3.l4i2-o6
5.0062-06
7.2872-06
1.0962-05
1.0472-05
8.421E-06
6.1822-06
4.2742-06
2.8442-06
2.0612-06
1.5722-06
1.1362-06
7.8142-07
5.5122-07
3.75^2-07
2.5652-07
1.7822-07
1.028E-07
4.497E-08
16
2.0852-16
3-7442-15
9.1702-14
4.8942-12
2.2802-10
7.4062-09
3.9502-08
8.640E-08
1.8902-07
4.0322-07
7.1052-07
1.2292-06
2.0372-06
3.1392-06
5.3282-06
9.8172-06
9.43IE-06
7.5002-06
5.4132-06
3.6912-06
2.7032-06
2.0752-06
1.5072-06
1.041E-06
7.3682-07
5.0332-07
3.4492-07
2.4032-07
1.391E-07
6.106E-08
17
7.2022-17
1.2952-15
3.1822-14
1.7052-12
7.9922-11
2.6172-09
1.4092-08
S.ioiE-08
6.8422-08
1.4762-07
2.6302-07
4.6202-07
7.8152-07
1.2382-06
2.2142-06
4.8542-06
9.1082-06
8.643E-06
6.7472-06
4.7942-06
3.5682-06
2.7662-06
2.025E-06
1.4o82-06
i.
-------
TABLE 6 Cont.
THOMANN'S TRANSFER MATRIX
U)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
21
1.8292-18
3. 2982-17
8.1372-16
4.3902-14
2.0752-12
6.8812-11
3-7552-10
8.3402-10
1.8642-09
4.0792-09
7.3842-09
1.3232-08
2.2972-08
3.7602-08
7.1182-08
1.8202-07
4. 6602-07
1.1322-06
2.477E-06
4.6';5^-o6
6.0742-06
5.9082-06
5.0202-06
3.8992-06
2.9793-06
2.1YT2-06
1.5862-06
1.1642-06
7.2542007
3.360E-07
22
1. 1042-18
1.9912-17
4.9132-16
2.6512-14
1.2542-12
4.1602-11
2.2722-10
5.0462-10
1.129E-C9
2.4732-0?
4.4792-09
8.032E-0?
1.396E-OS
2.289E-0?
4.344E-Cc
1.117E-07
2.8882-07
7.1232-07
1.602E-C-5
3-1782-C5
4 . 6992 -C-S
5.51.42-cS
5.1192-0?
4.1962-3;
3.306E-C=
2.475E-C-S
I.84i2-Co
1.3732-Cc
8.7352-:?
4.103E-07
23
6.5692-19
1.1892-17
2.934E-16
1.58*12-14
7.4952-13
2.4872-11
1.3592-10
3.0212-10
6.7592-10
1.1) 812-09
2.6852-09
4.8182-09
8.3822-09
1.376E-08
2.6162-08
6.7582-08
1.7592-07
4.3902-07
1.0072-06
2.0742-06
3.2902-06
4.2042-06
4.7472-06
4.2832-06
3.5432-06
2.7502-06
2.1052-06
1.6052-06
1.0482-06
5.0092-07
24
3.871)2-19
6.9892-18
1.7252-16
9.3162-15
4.4092-13
1.46iiE-ll
8.0022-11
1.7802-10
3.9822-10
8.7322-10
1.5832-09
2.8432-09
4.9502-09
8.1332-09
1.5^92-08
4.0152-08
1.0512-07
2.6462-07
6.1592-07
1.302E-06
2.1622-06
2'. 9012 -06
3.5872-06
3.9982-06
3.606E-06
2.9622-06
2.3632-06
1.8562-06
1.2552-06
6.1282-07
25
2.5422-19
4.5862-18
1.1322-16
6.1142-15
2.8942-13
9.6132-12
5.2562-11
1.1692-10
2.6172-10
5.7392-10
1.0412-09
1.8702-09
3.2572-09
5.3542-09
1.0212-08
2.6522-08
6.9662-08
1.7632-07
4.1412-07
8.8902-07
1.5152-06
2.0842-06
2.6892-06
3.2472-06
3.3872-06
3.0142-06
2.5332-06
2.0592-06
1.4452-06
7.2122-07
26
1.6442-19
2.9662-18
7.3232-17
3-955E-15
1.8732-13
6.2202-12
3.4022-11
7.5672-11
1.6942-10
3.717E-10
6.743E-10
1.2122-09
2.1112-09
3.4722-09
6.6232-09
1.7242-08
4.5412-08
1.1542-07
2.7302-07
5.9312-07
1.0292-06
1.44 22-06
1.9132-06
2.425E-06
2.7222-06
2.8662-06
2.6362-06
2.2592-06
1.6722-06
8.5872-07
27
1.0482-19
1.8922-18
4.6712-17
2.5232-15
1.1952-13
3.9692-12
2.1712-11
4.8302-11
1.0822-10
2.3732-10
4.3072-10
7.7412-10
1.3492-09
2.2202-09
4.2372-09
1.1042-08
2.9162-08
7.4382-08
1.7692-07
3.8802-07
6.8332-07
9.7012-07
1.3142-06
1.7202-06
2.0202-06
2.3182-06
2.5632-06
2.3992-06
1.9172-06
1.0232-06
28
6.7342-20
1.2152-18
3.0012-17
1.6212-15
7.6762-14
2.5512-12
1.3952-11
3.1052-11
6.9532-11
1.5262-10
2.7692-10
4.9792-10
8.6602-10
1.4292-09
2.7282-09
7.1202-09
1.8832-08
4.8172-08
1.1512-07
2.5432-07
4.5282-07
6.4942-07
8.9242-07
1.1952-06
1.44 52-06
1.741)2-06
2.1062-06
2.3522-06
2.1232-06
1.193E-06
29
3.3162-20
5.9832-19
1.4782-17
7.9832-16
3.7812-14
1.2572-12
6.8752-12
1.5302-11
3.4272-11
7.5232-11
1.3362-10
2.4562-10
4.2832-10
7.0522-10
1.3462-09
3.5232-09
9.3422-09
2.3992-08
5.7672-08
1.2872-07
2.3262-07
3-3792-07
4.7312-07
6.5132-07
8.1432-07
1.0362-06
1.3562-06
1.7202-06
2.2472-06
1.4162-06
30
1.177E-20
2.1242-19
5.246E-16
2.8352-16
1.3432-14
4.4632-13
2.4422-12
5.4362-12
1.2182-11
2.6742-11
4.8542-11
8.7312-11
1.5232-10
2.5092-10
4.7972-10
1.2562-09
3.3362-09
8.5932-09
2.0752-08
4.6652-08
8.521)2-08
1.2502-07
1.7722-07
2.4852-07
3.1732-07
4.1652-07
5.6912-07
7.6552-07
1.1292-06
1.4152-06
-------
After updating the intersection flows to obtain a new set of
Q-i i+1 t^ie updated transfer matrix is obtained from (20) .
Steady State Stability of the Transfer Matrix;
Definition: Let a-j_j be any element of A and let a^j be the
corresponding element of A1 (updated transfer matrix) . Sta-
bility of A (or A1) implies that for all i, j [a-J^-a^j !<16/dmax
and
In the above definition f, is an arbitrary parameter chosen to
reflect the accuracy of calculation of DO changes and dmax is
some arbitrary BOD load. To motivate the definition let
€ij = aij-aij. Intuitively we would like stability to imply
that £j_j is so small that changes in DO in section i are negli-
gible when some maximum BOD load, dmax, is applied to section j.
Theref ore,, we want
(24) ^ij
Since dmax>0, 6^>0 ; and ^ij = aij~aij we have
(25) JCij| = !a!.-a..] 6/d for all
In practive we generally take dmax equal to the maximum pollu-
ter input in all sections.
Experimental Results:
The ability to calculate the transfer matrix when needed adds
quite a bit of power and flexibility to the water quality con-
trol algorithm. The first problem we investigated was to re-
solve the by-pass piping problem using the new transfer matrix
given in Table 5. Recall that this case has linear constraints
A solution of the by-pass piping problem using the transfer ma-
trix of Table 6 can be found in [5] (see Figure 8 page 31) .
The new solution, with no updating of the transfer matrix, is
plotted in Figure 4 in a manner similar to that of Figure 8 in
reference [5] . Upon comparison of the two solutions we note
that: a) the overall flow patterns are similar; and b) the
final costs after eliminating duplicate piping are almost the
same ($2,089,148. as compared to $2,070,150. see Figure 9 of
[5]). Note that as we are now applying the piping cost equa-
tion the unadjusted solution of $2,447,995. is more indicative
of the cost with consolidated piping.
The next problem we investigated was to update the transfer
matrix and continue the solution carrying each by-pass piping
problem to the optimum. The optimal solutions obtained were
- 32 -
-------
^ X/^wWPV^/ /
*—- x:—/s'1^ i '
L_ --IZII" jL-il*""_^//
Optimal Solution f2,447,99s.
FIGURE 4
SOLUTION OF BY-PASS PIPING PROBLEM
USING NEW FIXED TRANSFER MATRIX
- 33 -
-------
transfer matrix cost
(Table 5) $2,447,995.
1st update $2,436,462.
2nd update $2,436,573.
3rd update (stable)
The total computer running time was 132 seconds with one re-
start at 60 seconds (some calculating time was lost since the
restart occurred while solving a by-pass piping problem).
The final solution is plotted in Figure 5. We can conclude
that for the by-pass piping problem on the Delaware Estuary
updating the transfer matrix is of little or no importance.
Since pure by-pass piping is the most extreme case of re-
arranging the flow pattern the above conclusion holds for any
problem configuration.
- 34 -
-------
section D.o. goal
1 0
1 1
1
14
15
It
optiMl Solution ¥2,436,573.
FIGURE 5
SOLUTION OF BY PASS PIPING PROBLEM
USING NEW VARIABLE TRANSFER MATRIX
- 35 -
-------
V. ESTIMATION OF THE DEPTH AND VELOCITY OF A RIVER
FOR A GIVEN RATE OF FLOW
Main reference for this study has been: F.M. Henderson,
Open Channel Flow (see [6]). No attempt has been made to
make specific references and for further clarification re-
garding the physics of the problem, the reader is referred
to the above text or to its equivalent. The motivation for
this study stems from the general equation for the behavior
of dissolved oxygen in an estuary as suggested by Thomann
[10], One way to improve the DO level is by controlling the
flow in the estuary which in turn determines the depth and
velocity. The rate of flow appears explicitly in Thomann's
equation. The depths are implicit in the volumes of the vari-
ous sections. If during the solution process of the general
model,, large changes of flow are generated, the matrix of trans-
fer coefficients has to be updated to correspond to the new
flow variables which in turn requires estimation of the new
depths resulting from the new rate of flow.
Four cases to be considered in this chapter are:
A. No inflows or outflows.
B. There are spot inflows (pipes) .
C. There are spot outflows (pipes).
D. There is a continuous inflow (tributary).
The underlying assuptions in all these cases are:
1. Water level does not vary across any cross^section.
2. Velocities of flow are low (or equivalently—slopes
are small) and the steady state velocity is sub-
critical.
3. Velocity of flow is uniformly distributed in any
cross-section.
4. Rectangular cross-section.
Additional assumptions do not apply to all four cases and they
are listed in their appropriate sections.
A. Regular flow - fixed rate of flow in the whole sections
The Bernoulli equation:
H = y + z +-
applies to small slopes.
H = total head
y = depth of water in the particular cross^section
- 37 -
-------
v = velocity of flow (average over the cross-section)
z = height of river bed above reference line
g = gravitational constant
For two consecutive cross-sections:
H. = H. ,-, + i. . .,, 1=0,..., n and i=0 at the beginn-
1 1"T"J_ 1^ 1 i X
ing of the first section and i=n at the mouth of the river.
n = number of sections into which the river is divided.
!L^ i+1 = losses between cross-section i and i+1.
Applying the Bernoulli formula to H. and H. gives:
v2 2 1
Yi + Zi + 2? =
Defining: S^ = average slope of river bed at section
= length of (i+l)st section
leads to:
Substituting gives:,,
Eliminating v.^ and vi+1 by observing that
Q = vA = vby
whe re:
A = cross-sectional area
b = width of cross-section^
results in:
- 38 -
-------
or:
+ - 0.
The loss function:
2
Ai,i+l = (Sf} average ' ^i+1 + CL \2g) average
over over
section i+1 section i+1
where :
Sf = friction slope = (-^J
"R
n = Manning's resistance coefficient
R = hydraulic mean radius = A/P
P = wetted perimeter
A = Cross-sectional area
c = Eddy coefficient = 2— for angles between 90
LJ C
and 180°
r = radius of curvature of center-line
c
The arithmatic average used to calculate Sf is a crude approxi-
mation but this is also the nature of the whole loss expression.
The estimate of CL is biased upwards severely in most cases,
but in the calculation here, the model does not take into account
bridges or islands. In addition, the larger the river, the lar-
ger the Eddy Coefficient for the same b/rc.
Substituting gives:
29
TT^ — O
{—
R
2 '2g 2g
Substituting R = ^ = .^ and v =
- 39
-------
1 /fin ^ AV r (b,+2y.)4/3
— • I r*i I AX I J- J_ T~
— i VT /id' **"• • -i i- •*•
2 1.4y !+l in/-3
(by).
+ TTT-TT- ] + 1+1 \-—^ + Y~'
(by)^{3 4g {by)i (by)i+l
To approximate c_ for angles smaller than 90 :
- L b +b. , b.+b.^n
2b _ 2 . i i+l) ._0_ = i i+l -. ff
°-'- -ri+1( ^ ,0° ^ 90o
where: 3 = angle of curvature (in degrees).
T ( Qn )2-AX [ 1 (^_ + -L)
2 i+l b. y.
)4/3] -
.
by i+l b
(^) + UJ
i+1 i+1
2
4gri,i+l 9° bY i bY 1+1
Equations (26) and (27) define a recursive relation between
y^ and y^+i- The initial condition required to solve this
recursive set of equations is obtained by assuming that at
the mouth of the river
yn = CONST.
y has to be known to begin the solution process.
Equations (26) and (27) have to be solved by using a numeri-
cal method. To update the longitudinal profile after a vec-
tor of flow-changes has been induced, there are two possibi-
lites :
a. Identify the new flows in every section and re-
solve equations (26) and (27) .
b. Approximate the changes in the y values as
- 40 -
-------
Looking at equations (26) and (27) as components of a vector
function which implicitly defines
': >
and using the implicit function theorem^leads to:
Equation (26) F, [ (
Equation (27) F2[(y±+l> Q) > (yi^i, i+l) ] =0
or:
= 0 = F
so that:
BQ
Only the first row of the matrix is needed to obtain yi ,
since:
dy
i+1
dQ Byi+1 dQ
[F1
dQ
- 41 -
-------
= A =
-
= B = -y
*
2 1.
+
Q-n
1.49b.y.
y.
b.y.
i i
90
dF,
= 1
The determinant of this matrix:
A = A-l - BC. = A + y. C.
y *
F ' -1
Only the first row of [ y. ).. . , ] is required to obtain the
Yi 1} 1}
first row of f ' ( ) . This row equals:
~^1 1 1 Yi
( - i. ) . ii = r.i _i i
- 42 -
-------
is obtained from Ci by interchanging i and i+1 in the y
and b values in C..
BQ
= G =
30
D
E
E+yjro
A+y.C.
•* i i
These expressions may be used to calculate
^i = ^i .dyi+l +
dQ ay, ,-, dQ
This is a recursive relation which may generate a large round-
off error in addition to a high truncation error which accumu-
lates as i decreases.
It may be both simpler and more accurate to use method j. (see
p. 40) to update the profile. The same is true in the other
three cases examined, as will be shown later.
- 43 -
-------
B. The Longitudinal Profile with Spot Inflows (Pipes):
Additional assumptions for this case:
1. There is a negligible region where the change in heights
and velocities spreads across the width of the river and then
the cross-sectional water profile is independent of the lateral
position.
2. The pipe is perpendicular to the direction of the main
flow. This means that the incoming flow does not add momentum
in the direction of the main flow.
AX,^,
AX .
AX '
I I
v
I
v .
i
i
Q+AO
Vi+l
Q > 0
FIGURE 6
SCHEMATIC AND NOTATIONS FOR SPOT INFLOWS
The solving equations is obtained by postulating that the mo-
mentum is preserved in the direction of the main flow:
Q'V . , = (Q + AQ) •v •
1,1 1,2
y^ 2 may k"2 calculated recursively from the known situation in
cross-section i+1.
Q + AQ
v.
AQ is assumed to be fixed and positive.
Q
Substituting: y. =
ljl b. .v. ,
Q + AQ
- 44 -
-------
gives:
y.
Q-Q
b. , (Q+AQ) (Q+AQ)
•*-}•*- K \7
£
= CQ+AQ} " Yi,2
i.l =
dQ
! 2 2- AQ-VJ,2
I+^Q) [Q2(i+^Q)
Q Q
r2-AQ ^ _
Q(Q+AQ)
dy.
+ 1 i =
dQ J
dy.
i.2.
dQ
Consider now the case where AQ changes. The vector of flow
changes will define the new flow in every section. The new
dpeths and their first derivatives may be estimated as follows:
5
A(
i
T V V
O 1 1 , 1
I
\
i \
\
\
— Q0+AQ
Yl,2 ^2 3
'3
dyo dYl dyia dy1<2 dy2 dy3
dQ dQ dQ dQ dQ dQ
Vo vl Vl,l Vl,2 V2 V3
FIGURE 7
SCHEMATIC AND NOTATIONS FOR DERIVATIVES OF SPOT INFLOWS
The situation shown above describes the flow before the change.
Suppose a change of A(A.Q) occurs. Then for all the cross-sec-
tions downs tream, the new flow is Q0+AQ+A(AQ) (assuming as be-
fore that A.Q+A ( AQ)>°) - For all those downstream cross-sections,,
~ dy
The new dy/dQ may be estimated from the basic recursive rela-
tion. This process of updating will proceed from the mouth of
the river up to the cross section next to the transition re-
gion (yi2) -
- 45 -
-------
For the cross sections upstream from the change:
A(AQ)
d(AQ)
Again, since the profile is obtained from a recursive relation
starting at the river-mouth, all the sections downstream from
the change in flow do not distinguish between the possible
sources of the change in flow, so:
cLQ
and:
i,l
i,2
i 9
(Note that for computational procedure, the expression of the
derivative dy^ i/d(AQ) is similar to the expression dy^ ^/dQ. )
The new value of dy^_ i/dQ may be calculated from the expression
for this derivative, ' substituting all new values which are known,
The rest of the sections upstream may now be updated using y^_ ^,
dyj_ i/dQ and the basic recursive relation up to the place
where the next change occurs.
If there are more than one inflow (or outflow) into the same
section and they may not be lumped together, the section may
be further broken into sub-sections. A convenient notation
procedure may be : All odd second indices are upstream and all
even ones are downstream from the location of discharge.
- 46 -
-------
C. The Longitudinal Profile with Spot Outflows (Pipes)
v.
V
Axi+l
i +19
Q-AQ
AQ > 0
FIGURE 8
SCHEMATIC AND NOTATIONS FOR SPOT OUTFLOWS
In this case; the specific energy (energy per unit of weight,,
as represented by the Bernoulli equation) is preserved provided
some condition is met, as shown later.
t
Assuming that preservation of specific energy is possible:
v2
v2
vi,l
2gb
1 ,
2? <
2
2 ^s known from the recursive process
- 47 -
-------
It will be shown that there are two possible situations:
1 There are three real solutions to the cubic equation , one
of them negative. Since the flow is assumed to be sub-cri-
tical the larger positive solution is the appropriate one.
The sign properties of the roots of the equation may be ob-
tained directly from the form of the equation. Yet, since
it has to be solved anyhow, the nature of the solutions will
be demonstrated by evaluating them.
2. Critical conditions develop before the location of change.
To develop the two situations, the starting point is the solu-
tion of Situation 1.
The Cardan method [13] for solving the cubic equation is used.
Accordingly, the following substitution:
Yi,l = H +i tyi,2 + 2?
—} "
yields:
3 1 . r ., ,fiiAfi ,2 2 +JM£)2_J_[ + i (glAfi ,2,^0
u> ~ T LY- o + -I \-K,, i \ M> 9<-r *lv ?7L-*i 9 9a T-^
r** j. 1 ^ —— jiD y . ^ y J^ /.. / J-, ^i £* ^3
2g 1} £
or:
To obtain the solution:
a - +
A ' 21 ~ 2 + \l(2 27
c _ D _
27 2 2 27
or:
=^-^+^ 1 -
el _ D ci //n 27DX2
27 2 27
- 48 _
-------
To get all three roots to be real, it is required that A and
B be imaginary, or, the term inside the square-root sign has
to be non-positive:
2c
2c
or:
1 27D
i. r ;> o
2c
Inspection of the equation shows that:
D > 0
c > 0
27D
and so - - > 0 is satisfied.
2c
2.
2c
Substituting the original, values for D and c:
or, 2
Substituting, the problem is now transformed into:
_4_ 3 r -, . Vi.2 , 3 /" XX 2
27 Yi,2 L-L + - I > 1 U.;
- 49 _
-------
When this condition is met, Situation 1 describes the case.
Namely: Ignoring losses in a transition area, the specific
energy downstream from the change may also be maintained up-
stream. The diagram shows the depth-specific energy relation
for different levels of specific flow. (The specific flow is
the flow per unit of width.)
Curve for
specific
flow up-
stream
from out-
flow
specific energy
Curve for specific
flow downstream from outflow
FIGURE 9
CRITICAL,, SUPER-CRITICAL AND SUB-CRITICAL FLOWS
AS DEPICTED BY THE DEPTH-SPECIFIC ENERGY RELATION
If the change in flow is large enough and y^_ 9 is close enough
to the critical value, Situation 2 will prevail. In this case,
the specific energies cannot be equal because the specific flow
upstream from the location of the change is too high to be able
to ever be in the energy-level of the lower flow. Since y^ -i
and y^ 2 are determined freely (they are not forced by any don-
trolling device) , y.j_ 2 will be determined by the energy required
to move the specified flow and to overcome losses (namely, the
basic recursive relation) and y^ ^ will determine the lowest
energy level possible for the specific flow Q/b. This lowest
energy is the critical one. So, in Situation 2, critical con-
ditions will develop upstream from the change and there will
be energy dissipation at the transition area so that the condi-
tions defined by the basic recursive equation prevail after the
transition. The flow above the critical cross-section will a-
gain be sub-critical.
- 50 -
-------
Solution of Situation 1:
c3 ,, 27D, . c3 / n 27D.2
A = -TJ (1 3) + yy / (1 3) - 1
27 2c 2/ V 2c
c n 27D>, . c^ /~ M 27D.T
— U ; -t- i / J. - u 3)
Zl 2c ll 4 2c
c^ ,-. 27D, . cj / , ,, 27D,2
B = 27- (1 -- 3) - i 27 / 1 - (1 - - 3)
2/ 2cJ /7 V 2c
c3 rT! 27DN2 " 7! 27D,2 iff _
A = -^=- / (1 -- o) + 1 - (1 -- o) ' A
27 J J
27
Similarly:
B = — IicT
B SL
where: / 27D
/I - (1 -
2^L , !
a = arc tg [ ^ _ 27D ] = arc tg / — J^2 ~ X
_3 / ' 1 "" •*'
2c V 2c3
i ^ 27D
for 1 ^ 3
2c
cos -^
c a ,
— cos - +
f cos f
-r cos - - v-r c sin
= - — c sin (-| + p)
c_
B = arctan • 3 = arctan ^
"/Tc/3
- 51 -
-------
So:
2 . / a , TT\
= -- c sin (- + ^)
Similarly:
2
6'
To determine the behavior of the roots,, a is explored:
L-27D\2
= arctan
27D
2c
3J
1 -
27D
2c3
For - r = 1^ the value of CT is determined from the original
2c
expression for A and B. In that case:
a
A =
. c3 c3
X 27 = 27
B = -i
So^ CT = 2
Since -1 <
7 27
for
, 3
27D
2c3
0
the above expression for a shows that
0 < a £ TT
Substituting gives:
0 < f £ H-L < § c
0 > - <
f > 0
There is also ordering as follows
M-i ^ M-3 > M-2
- 52 -
-------
Substituting
+f
yields :
Since the flow is assumed to be subcritical , (y. ) is always
i , ± 1
chosen.
y±jl = f (2 cosf + 1) =
= I ^1,2
i,2
cos [ arct*n
- 1 ]+l]
where:
27D = 27. _1_
2c3 " 2 2^
2g
i.l = 1 rdc (2 cos £
dQ 3 ldQU C°S 3
_ ± . a
3 dQJ
i,2 2g
dy.
^
dQ " dQ
dc dyi, 2 _ (Q-AQ) r_-| , Q-AQ . dyi,2.
-« .~ ^ -i xv ^ I ^~ - - ^3 /^
dQ
- 53 -
-------
dQ i _1+1 " dQ [/ n_27D,2
(1 _ 27D.2
2c3
/i 27D.2 1 '
= (i - —3) • j • / r
27D.2 -1
o 3'
2c
/-, 27D. 2 1 J- , 0 x , / 27. d ,D .
= (1 - o) • ^- / i ^~2^ ± (-^~");T^ L3)
2c3' z / ^-^^T -1
_3 dD -2 dS
"
27
6
2 / 2 7D 2 c
J 2c3
27 •»• dD 3D dC
9^3* /In 27D.2 L dQ ~ C dQJ
Z(- / -L V-L- 0;
- has been found before.
dO
^a c _
d° gb2
Substituting these expressions gives dy^ j_/dQ.
As mentioned in Case B, changes in AQ will be regarded as
changes in Q in the sections downstream from the location of
change. The only variables affected by a change in AO are
Yi,l and dYi.,l/dQ • yi,l maY be updated using a first order
approximation, for which dy-j_ j_/d(AQ) is necessary:
.in
dyi.2
(Q-AQ) rl , Q-AQ . dy -,
9 L -1- + „ J-^ J
g(by )^ Yi,2
1}z
- 54 -
-------
Substituting:
d(AQ) dQ
yields for-r- ar* expression very similar to that for
Appropriate substitution yields:
dc = dc 2 (Q-AQ)
dQ
1 3 dD _ 2 dc
da = 21 . f d(AQ) "^ d(AQ) ,
d(AQ) 2 H ,n 27D,2 [ _6 ]
2c3
V 2c
dc
, . -^ has been found before.
d(AQ)
L dD - _dC_
d(AQ) C d(AQ)
d ( AQ)
So:
= 0
dCT _ 21 X r-3D dc_
/ ,x->\ ~ Q r O Tn T L /-• r!/AA\J
2c y i -
2c3
dy. 7
Substituting into gives a value for the first derivative
which may be used to update y^ 2 •
Situation 2 :
The critical depth is:
- 55 -
-------
dY,
d(AQ) *
Transition between Situations 1 and 2:
If the expression:
(28) (£[y1>2 + J, (g)2 ()]-<)J
zg i, £
changes sign as a result of the change in /\Q, then a transition
between the two situations occurs.
If expression (28) changes from negative to positive, then y^_ ^
and dy-j_ ±/dQ are calculated as in Situation 1 with the new tff~
If expression (28) changes from positive to negative,, then re-
gardless of yj_ 2> Yl l is critical and dy^_ ]_/d^O = 0. If AQ
becomes negative_, the outflow becomes an inflow. Then the up-
dating has to be done in two steps:
1. To the situation where /\Q = 0;
2. Using the routine described in Case B.
The procedure is similar in the case of changing an inflow
into an outflow.
D. The Longitudinal Profile with Continuous Inflow (Tributaries)
The notational convention is:
Capital letters - for main stream
Lower case letters - for tributary
Q^q = rates of flow
V,v = velocities
M = total momentum
M = component of total momentum in direction of main flow
X.
B,,b = widths of main flow and tributary
ot = angle between main flow and tributary at intersection.
- 56 -
-------
FIGURE 10
SCHEMATIC AND NOTATIONS FOR CONTINUOUS INFLOW
Assumptions:
1. Transition regions have negligible dimensions and effects.
2. The depth of the tributary equals that of the main flow
and it is approximated by y^ 2•
3. The component of momentum at the direction of the main
flow is maintained. (Actually this may work only for low
velocities. In reality, preserving the momentum would have
created a path of flow which is forced to follow the path
of the main-stream river-bed. In that case, there is a
bending phenomenon which may tend to preserve the velocity.
There is little theoretical or empirical knowledge describing
the behavior of the water in the junction. In the case of
high velocities the depth of the water is not uniform
across the river. In addition, there are energy losses
which cannot be approximated.)
4. Both flow velocities are sub-critical.
- 57 -
-------
M = QV + gv cos
q X
q = -r— sin
v = Const. = qO = qO
by by. _
i, z
V = BY
9 2
2 q X
2QY
dX
§- Q2 S , /V sin 2
dM
•*£ 71 / Q Q \
dx ~ M o ~ bf;
whe re:
A = BY = cross-sectional area of main flow corresponding
to point X
S = river-bed slope
S= frictional slope
n
4/3 1.49-B-Y BY
K
. (JL+ 1)4/3
V ;
where n is Manning's resistance coefficient
j = - -T— sin
dx b
Equating the two expressions for . x gives:
dx
B2y3So _ (_n2_)2 y (I + 1)V3 =
- 58 -
-------
= - 2QY Sin a - Q2
2yi
x j •
2 v/l . i}4/3 /_IOx2. B sin 2 » . V2
W n' "*" \v«/->' ••».. * JL —
dX V1.49' iVY ' B' T vbQ' 2y7 .
•i- • ^
t, o -, 2q sin a
-
-------
VI. TRANSIENT EFFECTS IN AN ESTUARY
A. Introduction
The transient behavior of the DO resulting from a change in the
system is not considered by the approach to DO control which
utilizes the transfer coefficient matrix. This study attempts
to justify the above approach by demonstrating that under "real
world" conditions, the system representing the Delaware Estuary
converges essentially monotonically without developing any
pathological behavior. It is also demonstrated that the length
of the transient period (i.e.,, the time required to get within
+ 1% of the steady state solution and stay there) in the slower
sections is approximately one month.
B. Representation of -the Transient
Thomann's model [10] consists essentially of two systems of
differential equations, one describing the behavior of B.O.D.
and the other one describing behavior of DO. The second sys-
tem couples with the first one through the B.O.D. levels.
Theoretically, an analytic solution to the system of differen-
tial equations is available. Its derivation in this case, how-
ever, involves inversion of a 30 x 30 matrix (for the 30 sec-
tions of the Delaware) with analytic expressions (as opposed
to a matrix with numerical entries) . This fact provides a con-
vincing argument for numerical integration. In addition, the
result of an analytic solution is again a 30 x 30 matrix which
does not have any intuitive meaning, necessitating numerical
evaluation of the transient behavior under assumed conditions
of loading.
Alternative approaches to the numerical integration are either
to investigate the stability of the system (Pole-Zero analysis)
or to load it with various patterns of Dirac Deltas (spikes).
Both possibilities have been explored by Thomann ([10], p.27,
[12], p.354). The main shortcoming of these approaches lies
in the fact that they do not provide a simple answer to the
question: What would be the consequences of a change in the
system and would the transient create undesirable effects?
The approach adopted in this study tries to provide the answer
by sampling numerically a few responses under conditions which
are plausible in reality. This approach provides also an indi-
cation on the sensitivity of the system to changes in certain
parameters.
C. Assumptions
The transient reaction of a linear system does not depend on
the forcing function (the transient is the homogeneous solution)
- 61 -
-------
However, since a numerical time-pattern of the system's behavior
was sought, a simple forcing function which would give intuitive
meaning to the result, had to be assumed.
The most obvious one would have been to remove all sources and
sinks of oxygen and to evaluate, under proper initial conditions,
the mathematical homogeneous solution. The disadvantage of this
procedure lies in the fact that the physical phenomenon of sat-
uration with DO is ignored. Since, as explained later (Part D),
the system may be allowed to oscillate between the zero-level
of DO and the saturation level of DO, the following alternative
was preferred:
By removing all sources and sinks of oxygen except for the re-
aeration from the atmosphere and by assuming that the water a-
bove the first section and below the last one is saturated with
oxygen, (together with an additional simplifying assumption ex-
plained in the sequel), the system is driven to a steady-state
solution at the DO-saturation level.
The additional assumption concerns the in-flowing tributaries:
The system of differential equations describes conservation of
matter. Accordingly, the amount of flow into a section has to
equal the amount of water flowing out. This is easily achieved
by equating these two quantities in the computer program, how-
ever, this implies that the conditions in the in-flowing tri-
butaries are the same as those prevailing in the next section
upstream. This is a somewhat artificial assumption but it is
the simplest one that could be made if only the main stream is
to be solved for.
It is easily seen that under these assumptions the DO satura-
tion level is the steady state solution for Thomann's equation.
D. Computation Method
As stated before, no analytic solution was sought for the sys-
tem. Instead, a numerical integration subroutine was used to
generate a time pattern of the behavior of DO. The integrat-
ing subroutine utilized Hamming's modified predictor-corrector
method. This subroutine checks the possible truncation error
at each step of the integration and automatically adjusts the
size of the interval of integration if the error exceeds a
specified value.
The integration period was deliberately chosen as 200 days,
which is much longer than the transient period (as it turned
out to be), in order to discover beats (dampening of the oscil-
lation and then a renewed increase in the amplitude) if any
would have occurred.
It should be stressed that Thomann's major assumption was the
- 62
-------
linearity of the main system. This assumption may not be valid
when the system approaches its boundary conditions, namely, no
DO or saturation, the reason being that at zero DO a different
physical process takes over, while the saturation phenomenon is
generated by a stochastic process of exchange of oxygen between
the atmosphere and the DO. These conditions are not explicit
in the system of equations, so there is nothing inherent in the
system of equations to prevent the solution from oscillating
above the saturation value or below zero if the numerical struc-
ture of the equations is appropriate. As a result, one of the
basic problems was the adequacy of Thomann's model for describ-
ing the transient. This is the reason why it was not clear
apriori if a beat would develop, or if the system would try to
oscillate beyond its boundary values.
The solution to the second problem was to bound the DO values
and their derivatives in case of a bad oscillation. As shown
in Part E, the actual runs were designed to investigate the
necessity of constraining the DO levels, with the result used
to indicate the adequacy of the model in estimating the transi-
ent behavior.
E. Computation and Results
1. The conditions assumed correspond to those prevailing
in the summer since this is the critical period.
2. The starting conditions for the integration of the DO
levels were the average values in the summer of 1964
(Table 3, Schaumburg [9]). The numbers are given in
Table 7.
3. The saturation value was calculated as:
c =14.652-0.410222T+0.0079910T2+0.00007777T3
s
(Source — Pence, Jeglic, Thomann [8]).
Substituting T=25°C=?77°F
yields
c =10.605976 mg/A
S
So the value used for the integration was 10,6.
4. Since the numerical integration is a sequential process
and at every step values are computed sequentially (by
a DO-loop) for all sections of the river, starting with
1 and ending at 30, the errors (truncation and round-
off) tend to accumulate towards the later points in
time and towards the sections indexed by higher numbers.
_ 63 _
-------
The integrating subroutine estimates the error for each
component of the solution vector generated at each step
of the integration. Then the sum of those error terms
is compared with a predetermined value. If the error-
sum turns out to be too high, the integration interval
is halved and the step is repeated.
An option is offered by the subroutine to assign weights
to the individual error terms, thus stressing the accur-
acy of certain components of the solution vector.
Utilizing this option, an attempt was made to improve
the results by re-weighting the error terms of the
individual sections. The purpose was to increase the
sensitivity of the error-sum to those components cor-
responding to sections of the estuary which were indexed
by higher numbers.
The re-weighting was done proportionally to r? and to
x/rf , where n is the section number. In the case of /rf ,
the error term corresponding to section n was multiplied
by /n/ 30 -/ri . In the case of n^ the weight was n^/ 30 n2
2j _ 2-i .,
n=l n=l
5. Since the system turned out to be insensitive to the
variations tried, they will be listed and the results
given for one of them.
Sensitivities checked:
i) increasing the saturation value to 12.6 mg/ 1. and
decreasing it to 8.2
ii) forcing as opposed to not forcing the DO levels
to be between 0 and the saturation value .
6. As stated before, the results were insensitive to the
above variations. Re-weighting as n^ improved a little
the numerical behavior. (For explanation, see comments
following Table 7.) Forcing the DO to be at most equal
to the saturation value and the derivative to be 0 or
lower in case the section is saturated, introduced addi-
tional error and the system converged to a somewhat
lower value than the saturation level. The runs made
without these constraints demonstrated almost no oscilla
tion above the saturation value, indicating that the
linearity assumption is adequate in this sense.
7. Lengths of time for the system to converge to + 1% of
the saturation value are listed in Table 7. In this
run the saturation value was 10.6 mg/£, the system was
not forced to remain below this value and the error
terms were re-weighted as n2 .
- 64 _
-------
Section
Number
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
-Transient
Period
in Days
4
4
6
10
14
17
21
24
25
27
29
30
31
32
33
33
33
33
33
32
31
30
31
29
28
25
22
23
21
16
Starting Error
Conditions Range
mg/ & mg/ &
7.5 nil
7.0
6.2
5.1
5.1
5.6
5.4
5.3
3.0
2.2
2.3
1.5
1
1
1
1
1
1
1
2
2
3
4
4
5
5
6
7
7
8
.2
.1
.1
.0
.0
.2
.4
.2
.9
.5
.1
.8
.5
.8
.2
.2
.6
.0
11
197
759
2386
7201
13230
15816
13395
6799
743
56
11
n
"
"
11
11
11
x
x
X
X
X
X
X
X
X
X
X
10
10~q
10 ^
10 ^
10 b
10
10~n
10~^
"Is
10 R
10 b
TABLE 7
TRANSIENT PERIOD OF SECTION i, CORRESPONDING
STARTING CONDITIONS AND COMPUTATIONAL ERROR
The error-range (last column in Table 7) is the maximal absolute
value of the difference between the following two numbers:
i. The saturation value (10.6 mg/jO .
ii. The numbers generated by the computer for the per-
iod between the end of the transient phenomenon
(column 2 of Table 7) and the end of the integra-
tion process (200 days) .
- 65 -
-------
These deviations represent the effect of the cummulative error
generated by the numerical integration process. Decreasing
these deviations was considered as "improving the numerical
behavior" (see E.6) and this was the reason for re-weighting
the error terms as explained in E.5.
The maximal deviations occur in sections 25, 26 and 27 and they
constitute an error of about 1.5%. The deviations in the rest
of the sections are well below 1%.
These observations may serve as indications on the accuracy of
the numbers generated by the integrating subroutine and on the
validity of the results of this section.
- 66 -
-------
VII. POLLUTION TAXES AND CHARGES IN PERSPECTIVE
Previous work [4], presented a model to find a minimal treat-
ment cost solution to a water quality problem. A solution to
such a problem gives an efficient solution in terms of cost
minimization. However, such a solution may be costly to imple-
ment and enforce unless each polluter accepts the regional solu-
tion and his share of its costs. We now wish to discuss whether
there are decentralized mechanisms so that decisions by indivi-
dual polluters correspond to minimization of regional costs.
The possibility of the decentralization of decision making is
important in economics since enforcement costs are minimized
when the social optimum can be achieved through individual de-
cision making. In economics, taxes or charges of various types
are often suggested as a decentralized means of making indivi-
dual behavior correspond to a social optimum. Several taxa-
tion and charge schemes will be discussed below.
Examination of the literature concerning the economics of
water quality reveals the assumption that society has a right
to a certain level of water quality, or there is a socially op-
timum level of water quality- By implication, the polluters
should either
1. be forced to treat waste material at a high enough
level to maintain the quality required by society
or 2. pay the costs society incurs to maintain the required
water quality.
In this view, two types of tax schemes have been discussed,
those relating to enforcing a quality standard and those re-
lating to raising revenue. The first type can be called an
enforcement tax. One example of such a tax would be to charge
a polluter an extremely high tax (or fine) on wastes dumped
above a certain level. Thus, above this level, the polluter
would rather cease polluting than pay the fine. (To do this,
he could either curtail production or build a treatment faci-
lity.) Below this level, the polluter is in effect given free
rights to pollute. Such a tax would not be intended to raise
any revenue.
Another enforcement-type tax would be to have a uniform tax
per unit of waste and the tax would be computed so that a pollu-
ter would dump no more than a certain amount. This is illus-
trated below
T
price percent
effluent
D
Eo
- 67 -
-------
DD represents a polluter's demand to dump effluent as a func-
tion of the tax. If E is the amount of effluent to be allow-
ed under a quality standard, then the per unit tax which would
result in Eo is given by the slope of T. If the polluter would
want to produce at a higher level than that corresponding to Eo,
then again he would have to build waste treatment facilities.
A related concept would be to view the service of waste disposal
via a river as a public resouce to be sold at a price. The
price and amount sold is set to ration use of the service cor-
responding to the quality goals. In both these cases, the tax
revenue collected is "gravy" since the tax collecting agency
does not construct treatment facilities; the taxes serve only
to force achievement of the quality standards. Any treatment
undertaken under these enforcement tax systems would be at an
individual level and would be individually paid for. Because
of the individual emphasis, such tax systems would not result
in minimizing regional treatment costs.
The second type of tax scheme would be for the purpose of
raising revenue for treatment purposes. For instance, a pollu-
ter could be charged for waste disposal according to the social
cost of his pollution and the revenues collected could be used
either to alleviate the damages or pay for treatment. Here
again, this may not lead to the least cost regional solution
since it is always easier and cheaper to treat wastes before
they enter a stream than to have to remove them from the stream
once they are there. Alternatively, revenue could be raised
through charging the polluter for his use of a public or "joint"
waste treatment facility. It would not be necessary for the
charge to be collected by a governmental agency in this case.
A contractual agreement could be formed among several polluters
to operate a joint facility such as a regional treatment plant
or reservoir. The charge for use of the facility would be
based on cost. The advantage of this system is that the least
cost regional treatment system could be obtained on a voluntary
basis. However, in this case, there would first have to be some
sort of enforcement mechansim, either legal or through enforce-
ment taxes, to give incentives to polluters to have wastes
treated in the joint facility rather than dumping them. If en-
forcement taxes were used as incentives, then the treatment
charges would have to be lower than the tax; otherwise polluters
would pay the tax rather than treat their wastes in the regional
facility. The charges would also have to be low enough for the
polluters to choose the joint treatment alternative rather than
treating their own wastes individually. The problem is compli-
cated by the need to cover full costs of building and operating
the joint facility from the charges.
It is possible that a proposed tax could fall into both of the
above two categories as it affects polluters. For example,
suppose a uniform effluent charge of $.05 per mg/d is to be
levied on effluent per unit B.O.D. concentration to pay for
treatment costs. Because of differing manufacturing processes,
- 68 -
-------
locations and other parameters the cost of removing waste mater-
ial may vary for polluters. For this reason one polluter may
have an on-site B.O.D. removal cost of $.04 per mg/d and another
on-site B.O.D. removal cost of $.06 mg/d. In this case, the
second polluter would pay his effluent charges since his margin-
al removal cost is more than the cost of dumping it in the river
and the first polluter would treat his own waste. This is not
inconsistent with the least cost regional solution since the
first polluter can remove his wastes himself more cheaply than
in a regional plant; such a situation can arise when a polluter
is so far from a regional plant that piping costs overcome possi-
ble economies of scale. Here we see that the tax has served
two purposes. In the case of the first polluter he has been
forced to treat his own wastes. In the case of the second pollu-
ter revenue has been raised to defer the cost of treatment under-
taken by the "government" or tax collection unit. However, it
is clear that unless treatment plant sizes and charges are based
on the least cost regional solution, in general the result of
such charges may not cover costs of the facilities. That is,
in the above example, if a treatment plant were designed to
treat wastes of both the first and second polluter and the charge
set accordingly at $.05 per mg/d, then the plant would in fact
be too large since the first polluter would do his own treatment
rather than use the plant at that charge.
The determination of charges for the purpose of quality enforce-
ment is generally based on two criteria:
1. Equity between polluters;
2. quality resulting from implementation of charges meets
quality goals.
Johnson in [7] determines several types of effluent charges
which meet the second criterion above. The methods used were
called least-cost, single effluent charge, and zoned effluent
charge.
The least-cost method proposed by Johnson seeks to find the set
of effluent charges which will result in the minimization of the
sum of the taxation costs incurred by the individual polluters
subject to the constraint that the quality goals are met. This
approach could possibly result in a different effluent charge
for each polluter.
The second approach seeks to impose a single effluent charge in
terms of $/lbs. to all polluters, again subject to the constraint
that the quality goals are met. Another approach, the zoned
method divides the polluters between regions and assigns a single
effluent charge to each region to minimize the total taxation
cost over the region. Actuarlly, all of these methods are vari-
ants of the zoned method since the single charge method treats
the entire basin as one zone and the least-cost method in effect
- 69 -
-------
assigns each polluter to a zone. In all cases the least-cost
solution is strived for given the constraints on the charges
and the quality requirements. Johnson's charges are all exam-
ples of enforcement taxes. His cost minimization criterion
is a way of determining how polluting rights are to be distri-
buted and is thus an equity criterion. It would also be possi-
ble to devise other methods of assigning effluent charges. For
example, instead of zones the polluters could be divided with
respect to type, each type being assigned an effluent charge.
Assigning the same charge for each type gives another equity
criterion.
As already mentioned, there is no reason to assume that the
revenue obtained from a uniform effluent charge would exactly
meet the cost of treatment. Also, an enforcement type tax
would not in general result in minimizing basic treatment cost.
For these reasons, recently the revenue collection type of tax
has become more interesting since the thrust has been to reduce
total basin costs significantly through the use of treatment
alternatives such as by-pass piping, regional plants and flow
augmentation. Each of these alternatives creates problems of
joint funding and operation of facilities. In such a case the
revenue necessary for the operations of these joint facilities
must be obtained from the polluters concerned. Charges of the
revenue collection variety are usually based on three criteria:
1. Equity between polluters;
2. quality goals maintained;
3. revenue sufficient to cover costs incurred by tax
collecting unit.
Graves, Hatfield, and Whinston [5] construct a pricing scheme
in their paper on by-pass piping which is of this type. Their
scheme is based on an equity criterion of charging polluters
according to damages. The scheme is not self-enforcing and so
would have to be preceded by some enforcement mechanism to give
incentives to meet the quality goal through by-pass piping.
The minimum cost piping arrangements is first solved for sub-
ject to meeting the quality goals and conservation of flow con-
straints. The dual variables associated with the flow conser-
vation constraints measure the marginal cost of extra effluent
flow from an individual polluter to the least cost basin solu-
tion. Suppose each polluter were charged according to the dual
variable per unit effluent. This would lead to a charge of
- 70 -
-------
denotes the optimal effluent flow from polluter k from solu-
tion of the least cost programming problem, TC is the least cost
solution, and the partial derivative is evaluated for the least
cost solution. Total revenue from this charge would be
y 2£!C £E = £
which in general would not equal TC. The expression in (29) i
only guaranteed equal to TC for all values of F|, k=l,n if the
least total cost function is linear homogeneous in the effluent
flows .
In the small 5 polluter example in the Graves, Hatfield, Whin-
s^ton paper [5], the least-cost solution is $55,577. The total
revenue collected using an effluent charge e.qual to the value
of the dual variables associated with the flow consistancy con-
straints is $179,499. In terms of the notation above
TC = $55,577
5
Z TP, = $179,499
k=l k
or
5
TC < Z TP, .
k=l K
In the Graves, Hatfield and Whins ton paper the revenue is forced
to equal cost by redefining the total payment of polluter j.
n
Tp* - (TP ./ 2 TP, ) TC
3 : k=l *
This automatically gives
Z TP* = TC .
k=l k
This cost allocation is based on achieving the quality goals
at least cost, costs are covered from the revenues, and there
is an equity criterion. The equity criterion is based on pollu
ters causing the same marginal damage paying the same per unit
- 71 -
-------
effluent charge, where marginal damage is defined to be the
increase in the basin solution cost due to an extra unit of
effluent.
Although the above is one type of equity criterion, it is clear
that there are other types,, such as the same type of pollutant,
or the same charge for polluters in the same zone. The notion
of what is meant by equity needs to be further examined. Equity
is included as a criterion for both revenue and enforcement
taxes. Its importance is due both to our notions of fairness
and to the difficulty of enforcing and collecting the tax or
charge when it is not felt to be equitable by those who will
be charged. In general, the notion of equity is based on two
principles:
1. Polluters who cause more damage should pay more.
2. Polluters who cause equal damage should pay the same.
The problem which is evident from these requirements is that
of interpreting "damage". The proposal to charge polluters
of the same type the same amount per unit weight makes the
implicit assumption that polluters who dump equal weights of
effluent of the same composition are causing equal damage.
However, this assumption may be false, for instance when the
flows (hence concentrations) of their effluents differ. The
location of a polluter also determines how much damage he will
cause. A polluter located in an area of high water quality
(and/or high flow) will do considerably less damage than the
same polluter located where the water quality is poor (and
stream flow is low). Damage also depends on what is located
downstream of a polluter. Damage will be much less if a pulp
mill is located downstream than if a municipal water plant is
located downstream from a polluter. Under closer examination
it is not at all clear that any measure of effluent composition
alone reflects pollution damage. In an economic sens.e, the
only measure of damage which is really relevant is the amount
of additional cost necessary to maintain desired quality goals
due to the existence of polluter k. If polluter k dumps a
certain effluent but causes no additional treatment cost to
the system (i.e., quality constraints are not violated), then
his existence is irrelevant as far as any treatment cost scheme
is concerned.
The Graves, Hatfield, and Whinston [5] paper recognized cost
in this way as a reflection of damage through their use of the
dual variables of the programming model to generate equitable
shares of the total basin costs. The problem associated with
this measurement of damage, as with any sort of marginal con-
cept, is that this information is only "local" in nature. Since
the relationship between the optimal basin solution, and the
flows from the individual polluters is non-linear, the marginal
cost at the optimum solution of an additional unit of flow from
- 72 -
-------
polluter i is no reflection of the added cost to meet the stan-
dards if effluent from a polluter were to be increased by more
than one unit. Thus, the equity criterion used in the Graves,
Hatfield, Whinston paper is not completely satisfactory.
Incremental Cost Allocation Scheme
In the previous sections of this paper we have indicated several
faults associated with the pricing and cost allocation systems
presented in the literature. These faults can be summarized as
follows:
1. Failure of revenue collected to meet total costs of
treatment.
2. Failure of some pricing systems to be based on the most
efficient (least-cost) combination of treatment methods.
3. Failure of pricing schemes to satisfy acceptable equity
criteria.
In this section we will present a cost allocation system which
we feel adequately deals with these problems. This approach takes
as given the present composition of polluters in the river basin
and the given quality goals and attempts to allocate the costs
of a basin wide system. Both the costs of individual polluters
and the joint costs of using regional treatment plants are con-
sidered in the allocation scheme.
After an optimal solution to a basin polluter problem, such as
that described earlier, has been obtained, it is not at all
clear what part of the cost of that system of dams, treatment
plants, and pipelines can be attributed to the fact that a par-
ticular polluter is in the basin. That is, when several pollu-
ters are present, it is not clear who has the responsibility
for violating the quality constraints and hence the exact damage
costs attributable to any one polluter is not well defined. If
we could measure the cost that an individual polluter causes to
the system, then we would have some basis to decide what his
contribution to the total cost of the abatement system should
be. We will now describe a way of allocating the total system
costs to polluters based on the idea of each polluter paying
the incremental cost he causes the system.
Consider a basin with only two polluters designated 1 and 2.
Denote the least-cost solution to the basin problem (29) by
C(l,2). The incremental cost caused by polluter 1 being in
the basin, given that polluter 2 is in the basin, is defined
to be
IC(1) = 0(1,2) - C(2)
where C(2) is the cost of polluter 2 treating his own waste.
- 73 -
-------
Both the values of C(lj2) and C(2) can be computed using the
programming model in the previous section.
A logical tax based on each polluter paying his social cost
would then appear to be to charge polluter 1 1C (1) and pollu-
ter 2 1C (2) , where
1C (2) = C(2) .
We note that
1C (1) + 1C (2) = C (1,2)
which means that the total cost of the basin project has been
covered by this tax scheme. However, the following question
now arises: What if the polluters are considered in reverse
order? Consider the following definitions of 1C* (1) and 1C* (2 )
1C* (2) = C(l,2) - C(l)
If polluter 1 is charged IC*(1) and polluter 2 is charged
1C* (2) the total cost of the basin system is covered as before,
but the charge to each polluter differs in these two schemes.
As an example, consider a case where if either of the polluters
are alone in the basin the quality goals are met without treat-
ment and if they are both in the basin the goals are violated.
Assume that the optimal solution when both are present is to
build a reservoir for flow augmentation at a cost of C(l,2) =
$100. The values of C(l) and C(2) are both zero since the stan-
dards are met if either polluter is absent. if the first cal-
culation of incremental cost is used the charges for polluters
1 and 2 are :
1C (1) = $100
1C (2) = $0
If the second calculation is used the charges are:
IC*(1) = $0
1C* (2) = $100
Polluter 1 clearly would prefer the second allocation over the
first cost allocation and by the same token polluter 2 would
prefer the first allocation. Since either ordering of the
polluters is equally plausible, we may average all the possible
incremental costs which each polluter could cause to obtain a
cost allocation of
p(l) = & IC(1) + & IC*(1)
= 4 (100) + £ (0) = 50
- 74 -
-------
and
p(2) = | 1C (2) + | 1C* (2)
= £ (0) + i (100) = 50
The sum of the payments from the polluters covers the total
basin cost of 0(1,2) = $100. However, the equity problem has
been taken care of by considering all possible orders of the
polluters, in this case (1 2) and (2 1), and averaging over
possible definitions of incremental cost for each polluter.
This approach is easily extended to a case of three polluters
where there are six possible orders for calculating incremental
costs for each polluter. For example, considering polluter 1,
the following orders are possible.
Order 1: (123)]
V polluter 1 before polluters 2 and 3
Order 2: (132) )
Order 3: (213)} polluter 1 after polluter 2
Order 4: (312) } polluter 1 after polluter 3
Order 5: (231)]
V polluter 1 after polluters 2 and 3
Order 6: (321) I
Corresponding to each of these orders is a way of defining in-
cremental cost for polluter 1 (the subscript denote s the order
under consideration) :
= C(l,2) - 0(2)
= C(l,3) - C(3)
= C(l,2,3) - C(2,3)
= C(l,2,3) - C(2,3)
For polluter 2, the following incremental cost definitions are
possible :
IC1(2) = C(l,2) - 0(1)
IC2(2) = 0(1,2,3) - 0(1,3)
- 75 -
-------
IC3(2) = C(2)
IC4(2) = C(l,2,3) - C(
IC5(2) = C(2)
IC6(2) = C(2,3) - C(3)
And for polluter 3,
= C(l,2,3) - C(
IC2(3) = C(l,3) -
IC3(3) = C(l,2,3) - C(l,2)
IC4(3) = C(3)
IC5 (3) = C(2,3) - C(2)
IC6(3) = C(3)
Giving each of the possible orders equal weight, the average
incremental cost for each polluter is obtained:
p(i) = j- 2 1C. (i) , 1=1,2,3
6 j=l ^
Note that
3
Z p(i) = C(l,2,3)
1=1
So that total costs would be covered by charging each polluter
the amount p (i) .
Generalizing to the case of n polluters where there would be
n! orderings, the cost allocation formula is given by
(30) p(i) =i 2*IC. (i) , 1=1, ...,n
n- j=i 3
where ICj (i) is the incremental cost due to i in ordering j.
Also,
n
2 p(i) = C(l,2, ...,n)
1=1
This cost allocation scheme is discussed in [14]. Also in [15]
- 76 -
-------
there is a theoretical development of the cost allocation me-
thod described. It is demonstrated that formula (30) is a
unique formula which satisfies the following axioms:
1. Charges for the use of the basin treatment system must
cover costs.
2. A user's charges are based only on the incremental
cost caused by the user.
3. The user's charge is independent of the ordering or
labeling of the users.
4. If a user's incremental costs increase by \, then his
charges will increase by \.
The equity properties of this scheme are evident from the above
axioms 2-4.
The following example illustrates cost allocation for the five-
polluter case. In the report [5] the solution of a small scale
pure by-pass piping problem is given on page 90. The dual solu-
tion for this problem is
XL = $3285/MGD xg = -$1, 929, 780 .MG/L
X2 = $11970/MGD X10= °
x3 = $6445/MGD x^ 0
X4 = $2865/MGD
X = $2381/MGD
The extremal value can be calculated as follows
Q± = 19.7 Q1x1 = 64,714 K-jXg = -147,107
Q2 =
Q3 =
Q4 =
Q5 =
7.0
3.0
6.1
4.0
Q2x2 = 83,790
Q3x3 = 19,335
Q4x4 = 17,476
Q5X5 = 9,524
194,839
-147,107
Kx = .07623 194,839 $ 47,732
A pricing scheme suggested in [5] was to use the values of the
xiQi to give the proportion of the optimal solution to be paid
by each polluter. This gives
- 77 -
-------
Polluter % Paid Payment/Year
1 33.2 $15,643
2 43.0 20,261
3 9.9 4,664
4 8.9 4,146
5 4.8 2,261
It would be interesting to compare the costs from the global
scheme to those given above. This is economically feasible
since we only have five polluters. The coalition costs were
obtained from the optimal solutions of 31 L.P. problems and
are given as
C(l,2,3,4,5) = $47,119 €(3,4,5) = 0
C(l,2,3,5) = 30,763 C(l,2,5) =$25,238
C(l,2,4,5) = 29,365 C(2,3,4) = 16,690
C(l,3,4,5) = 0 C(l,3,4) = 0
C(2,3,4,5) = 18,939 C(l,2,4) = 27,115
0(1,2,3,4) = 37,903 C(l,2,3) = 27,844
C(l,3,5) = 0 C(l,2) = 22,988
C(l,4,5) = 0 C(l,3) = 0
C(2,3,5) = 14,812 C(2,3) = 12,563
C(2,4,5) = 14,084 C(l,4) = 0
C(2,4) = 11,902 C(5) = 0
C(3,4) = 0 C(4) = 0
C(l,5) = 0 C(3) = 0
C(2,5) = $11,210 C(2) = $10,381
C(3,5) = 0 C(l) = 0
C(4,5) = 0
The polluters not in the coalition were allowed to return their
wastes to the estuary at no cost which in effect made them non-
existent. The global scheme gave the following cost allocation:
- 78 -
-------
Polluter Local Method Global Method
1 $15,643 $10,372
2 20,261 24,474
3 4,664 5,159
4 4,146 4,633
5 2,261 2,478
Application of Pricing Scheme
In the last subdivision we presented an alternate cost alloca-
tion mechanism to be used for regional treatment systems. In
this subdivision we apply the scheme to be the Delaware Estuary
Mode1.
We now wish to divide the total cost of the regional system,
in this case the cost of the plants and pipes, between the
polluters using the cost allocation method explained in the
last subdivision. This necessitates the calculation of C(G),
G <^ N where G ranges over possible subsets of polluters. The
interpretation of C(G) in terms of river basin model is not
straightforward. For example, how do we compute C(j) and what
do we assume about polluters I^j. Are they non-existent, are
they treating at some uniform level? Are they dumping their
effluent with no treatment? An infinite number of assumptions
could be made concerning the behavior of the polluters not in
the coalition to be examined.
The only restraint on the behavior of polluters outside of the
group being considered is that they are treating in a manner
which will allow the polluters to gain feasibility of the basin
problem (i.e., the quality constraints can be satisfied for
some level of treatment of wastes). For example, if the pollu-
ters outside of the coalition are assumed to be dumping all
wastes directly into the stream it is very likely that several
sections of the river will violate the quality constraints no
matter what the polluters do. One procedure which can guaran-
tee the feasibility of the problem is to find the least cost
solution to the basin problem, C(N), and assume that all pollu-
ters outside of G are required to treat their effluent to the
same level as they did in the solution for C(N).
The term C(N) can now be interpreted as the least-cost the
polluters in G must pay to maintain the quality goals given
that the other polluters in the basin are treating on site at
the same level required for the least cost base solution. The
polluters in G can make use of regional plants and piping. In
the special case where G contains only one polluter piping is
- 79 -
-------
the only alternative method available. In some problems no
savings over on-site treatment may result for some coalitions
of polluters. For example, if polluters in the coalition are
so far apart that piping costs overwhelm any economies of scale
gained from a regional plant, the best possible alternative for
these polluters is to treat on site. In such a case, the cost
to the coalition would just be the sum of the individual on-
site costs, or
C(G) = Z C(i)
ieG
- 80 -
-------
VIII. REFERENCES
[1] Coddington, E. A., "An Introduction to Ordinary Differential
Equations/" Prentice-Hall, 1961.
[2] Geoffrion, A. M., "Generalized Benders Decomposition," Publica-
tion forthcoming in Journal of Optimization Theory and Applica-
tion, 1972.
[3] , Graves, G. W., "Multicommodity Distribution System
Design by Benders Decomposition," Forthcoming in Working Paper
No. 181, Western Management Science Institute, October 1971.
[4] Graves, G. W., Hatfield, G. B. and Whinston, A., "Mathematical
Programming for Regional Water Quality Management," Water Pollu-
tion Control Research Series, U.S. Department of the Interior,
Federal Water Quality Administration, August 1970.
[5] , Hatfield, G. B., Whinston, A., "Water Pollution Con-
rol Using By-Pass Piping," Water Resources Research, Vol. 5, No.
1, February 1969.
[6] Henderson, F. M., "Open,Channel Flow," Macmillan Series in Ci-
vil Engineering, 1966.
[7] Johnson, E. L., "A Study in the Economics of Water Quality Man-
agement," Water Resources Research, Vol. 3, No. 2, 1967.
[8] Pence, G. D., Jeglic, J. M., Thomann, R. V., "Time-Varying
Dissolved-Oxygen Model, " ^Journal of the Sanitary Engineering
Division Proceedings of the American Society of Civil Engineers.
Vol. 94, No. 5915, April 1968.
[9] Schaumburg, G. W., "Water Pollution Control in the Delaware
Estuary," Submitted to the Department of Applied Mathematics
in partial fulfillment of the requirements for the degree of
B.A. with honors, Harvard College, May 1967.
[10] Thomann, R. V., "The Use of Systems Analysis to Describe the
Time Variation of Dissolved Oxygen in a Tidal Stream," A disser-
tation in the Department of Meteorology and Oceanography sub-
mitted in partial fulfillment of the requirements for the degree
of Ph.D. at New York University, February 1963.
[11] , "Mathematical Model for Dissolved Oxygen," Journal
of the Sanitary Engineering Division Proceedings of the American
Society of Civil Engineers. No. 3680, October 1963.
[12] , "Recent Results from a Mathematical Model of Water
Pollution Control in the Delaware Estuary," Water Resources
Research, Vol. 3, pp. 349-359, 1965.
- 81 -
-------
[13] Uspensky,, J. V., Theory of Equations. Me Graw-Hill, 1948.
[14] Loehmarij E., Whinston., A.,, "A New Theory of Pricing and
Decision-Making in Public Investment," Forthcoming in Bell
Telephone Journal of Economics and Management, Fall 1971.
[15] , , "An Axiomatic Approach to Cost Alloca-
tion for Public Investment^" Forthcoming.
_ 82 _
-------
1
Accession Number
w
5
2
Subject Field & tjroup
06 A
SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Organization
California University, Los Angeles
Title
Extensions of Mathematical Programming for Regional Water
Quality Management
10
Authors)
Glenn Graves
16
21
Project Designation
16110 EGQ
Note
22
Citation
Water Pollution Control Research Series, EPA, 16110 04/72 EGQ
23
Descriptors (Starred First)
*Regions, *Computer Models, *Mathematical Models, Networks,
Water Quality Control, Water Pollution Treatment, Optimum
Development Plans
25
Identifiers (Starred First)
*0ptimal Configuration of Regional Treatment Plans,
Non-linear Programming
27
Abstract
This report extends the work contained in the earlier work, "Mathematical
Programming for Regional Water Quality Management." A mixed integer, con-
tinuous variable non-linear programming model is developed which promises
to be much more realistic and effective in selecting an optimal configura-
tion of regional treatment plants and pipe juncture nodes. A new Bender's
type decomposition for the non-linear model is also presented as an appro-
priate and feasible solution technique. In the area of the physical response
of the environment, the earlier work linked the dissolved oxygen profile and
the control options through a constant transfer matrix obtained as a steady
state solution to a system of differential equations. This work shows how
the transfer coefficients can be recalculated periodically in the course of
the system optimization. Full scale computations are carried out demonstrating
this capability. The theoretical framework for including two additional treat-
ment methods, mechanical re-aeration and flow augmentation via reservoirs is
also developed.
Abstract
Institution
WR:102 (REV. JULY 1969)
WRSIC
SEND WITH COPY OF DOCUMENT. TO: WATER RESOURCES SCIENTIFIC I N F ORM A T I ON G E NT ER
U.S. DEPARTMENT OF THE INTERIOR
WASHINGTON, D. C. 20240
* OPO: 1970-389-930
------- |