WATER POLLUTION CONTROL RESEARCH SERIES
16110 FPX 08/70
Mathematical Programming
for
Regional Water Quality Management
U.S. DEPARTMENT OF THE INTERIOR • FEDERAL WATER QUALITY ADMINISTRATION
-------
Mathematical Programming
for
Regional Water Quality Management
by
Graduate School of Business Administration
University of California
Los Angeles, California 90024
for the
FEDERAL WATER QUALITY ADMINISTRATION
DEPARTMENT OF THE INTERIOR
Program #16110 FPX
Formerly WP-01210
August, 1970
For sale by the Superintendent of Documents, U.S. Government Printing Office
Washington, D.C. 20402 - Price $1.29
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Reports describe
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, develop-
ment, and demonstration activities in the Federal Water
Quality Administration, in the U. S. Department of the
Interior, through inhouse research and grants and con-
tracts with Federal, State, and local agencies, research
institutions, and industrial organizations.
Inquires pertaining to Water Pollution Control Research
Reports should be directed to the Head, Project Reports
System Planning and Resources Office, Office of Research
and Development, Department of the Interior, Federal
Water Quality Administration, Room 1108, Washington, D. C,
202H2.
-------
ABSTRACT
This report is an application of a non-linear programming
algorithm to the problem of optimal water quality control in an
estuary. The mathematical model that was developed gives the
solution to the general mixed case of at-source treatment, re-
gional treatment plants, and by-pass piping. The non-linear
algorithm is developed in considerable detail and a sample pro-
blem is worked out. Detailed results are presented for a pilot
problem to illustrate the method of solution. Actual data from
the Delaware Estuary is used to solve a large scale problem and
the solution is given. The results indicate that a regional
treatment system for the Delaware Estuary is superior, in terms,
of total cost, to other proposed schemes. This report was sub-
mitted in fulfillment of Grant No. 16110 FPX between the Federal
Water Quality Administration and the University of California,
Los Angeles, directed by Dr. Glenn W. Graves, assisted by Andrew B.
Whinston and Gordon B. Hatfield.
-------
Section I.
Section II.
Section III.
Section IV.
Section V.
Section VI.
Section VII.
Section VIII.
TABLE OF CONTENTS
List of Tables
List of Figures
Conclusions and Recommendations
Introduc tion
Mathematical Development*
Derivation of the Non-linear Problem
Generation of the Underlying Array of the Local L.P.*
Development of Treatment Costs
A Small Scale Problem
Cost Analysis for the Delaware Estuary
Economic Considerations
Appendix A
Appendix B
Indicates mathematical chapters.
-------
LIST OF TABLES
TABLE TITLE PAGE
1 Flows and Concentrations for Example Estuary 41
2 Distance Polluter to Section for Example Estuary 42
3 Distance Polluter to Plant for Example Estuary 42
4 Distance Plant to Section for Example Estuary 42
5 Transfer Coefficients for Example Estuary 43
6 DO Changes for Example Estuary 43
7 Example of Priority Class 52
8 Matrix Derivatives 54
9 Polluter Cost and Discharge Data 66
10 Comparisons of Estimated Plant Cost 78
11 Priority Class Configuration 1 84
12 Polluter Data for Example 85
13 Solution Data for Example Configuration 1 86
and Priority Class 3
14 Solution Data for Example Configuration 1 87
and Priority Class 3
15 Priority Class Configuration 2 88
16 Solution Data for Example Configuration 2 89
and Priority Class 1
17 Solution Data for Example Configuration 2 89
and Priority Class 3
18 Priority Class Configuration 3 90
19 Solution Data for Example Configuration 3 91
and Priority Class 1
20 DO Standard of 3 mg/1 94
21 Cost Comparison of Treatment Schemes 95
22 Transfer Matrix for the Delaware 97
23 Mileages from Polluters to Sections on the 100
Delaware Estuary
24 Mileages from Polluters to Regional Plants on the 102
Delaware Estuary
25 Mileages from Regional Plants to Sections on the 103
Delaware Estuary
26 Solution for Treatment at the Polluter for the 104
Delaware Estuary
27 Data for Optimal Solution 108
-------
LIST OP FIGURES
FIGURE TITLE PAGE
1 Dynamics of Using the Model 4
2 Graphic Solution of Non-linear Example 24
3 Graphic Solution of First L.L.P. 28
4 Graphic Solution of Second L.L.P. 31
5 Graphic Solution of Second L.L.P. 37
After Parametric Adjustment
6 Physical Layout for Example 41
7 Piecewise Linear Cost Curve 69
8 Smooth Curve Fit of Cost Data 70
9 One Price Break 72
10 Two Price Breaks 72
11 Three Price Breaks 72
12 Plant Costs 75
13 Solution of Example Configuration 1 and 86
Priority Class 1
14 Solution of Example Configuration 1 and 87
Priority Class 3
15 Solution of Example Configuration 2 and 89
Priority Class 3
16 Solution of Example Configuration 3 and 90
Priority Class 1
17 Map of Delaware Estuary 93
18 Optimal Solution (map) 107
-------
CONCLUSIONS AND RECOMMENDATIONS
The Delaware Estuary provided an excellent test for our
general non-linear mathematical programming model. At the op-
timal solution the three alternative treatment modes, treatment
at the polluter, by-pass piping, and treatment at regional plants
were all utilized. This would indicate that none of the treat-
ment modes dominates the others in a heavily polluted estuary.
Our final solution annual cost of 2.291 million dollars, to
achieve a 3 mg/1 dissolved oxygen standard throughout the estuary,
compares very favorably with the 4.1 million annual cost of
optimal treatment at the polluter, the next best known solu-
tion. The currently favored policy of requiring "uniform
treatment" at all polluters would require an annual expendi-
ture of 8.153 million to achieve the 3 mg/1 quality standard.
The solution time of about ten minutes for the complete model
makes it practical to repeatedly solve the problem with vari-
ous plant and pipe configurations, various quality goals and
different cost estimates. It would seem, therefore, in terms
of both utility and practicality, the non-linear programming
model and solution techniques developed here provide an effec-
tive management tool in pollution control.
The difficulties encountered in solving the non-linear pro-
gramming model were not due to the sheer size of over 2,000 varia-
bles and 80 constraints alone. The range of values of the trans-
fer coefficients made the scaling of the numbers for single pre-
cision (6 decimal digits) calculation quite sensitive and com-
-------
pounded the usual numerical analysis difficulty of round-off error
propagation. The economies of scale (mathematical concavity) in-
herent in piping costs and treatment plant costs also proved
troublesome. In fact the opening and closing of regional plants
and pipes hinges on very delicate physical considerations made in
evaluating indeterminate cases in the first partial derivatives
of the cost function. Many innovations in non-linear programming
solution technique were developed and incorporated in the machine
code to achieve the results and solution time reported here. The
technique of only generating columns of the underlying local linear
programming problems as needed, in addition to substantially reduc-
ing computation time as a direct effect, indirectly contributed
to the effectiveness of the other innovations. The use of priority
classes for the variables helps convergence, cuts computations,
permits the solution of specialized sub-problems such as pure treat-
ment at the polluter, and makes possible the opening of plants and
pipes in spite of their concave cost functions. The use of logi-
cal techniques for range restricted variables and piecewise linear
cost functions instead of expanding the size of the problem again
substantially cuts computation time. Perhaps, the most signifi-
cant innovation, however, was the use of a quantified "relaxation"
technique. Instead of obtaining a complete optimal solution to
each local linear programming problem in the non-linear step-wise
procedure, a control parameter determines the degree of optimiza-
tion. The most significant effect of this procedure is to lift
the classical curse of "zig-zagging" in gradient optimization pro-
-------
cedures. In addition, it permits decisive control of round-off
error propagation.
In view of the successes achieved to date, we would recommend
proceeding with further development of the approach. In particular,
development of a "global" approach to deal with the combinator-
ial problem of treatment plant location and the non-convexity
of the cost functions. Also a coupling of the present techni-
ques with an up-dating procedure for the transfer coefficients
which describe the physical quality response of the estuary
would contribute to the realism of the model.
-------
I. INTRODUCTION
It has become apparent that major changes are desirable
in the institutional structure for water quality management
in the United States. The desirable structure would be a
regional or basin-wide water quality management authority
[1, P. 75].
A single agency would control all discharges (industrial
and municipal) and operate all treatment plants in a region.
It would construct new regional treatment plants in optimum
locations and control the distribution and re-distribution
of treated and partially treated wastes. The authority would
be responsible for finding and implementing the least cost
solution of meeting the stream (or estuary) water quality
goals. The question of setting water quality goals relates
to social and economic desires and needs. This question
need not be resolved by the authority but it could make
valuable contributions to the rationality of the process. When
a change in goals is proposed, the authority will determine
the optimum solution to meet the changes in addition to the
cost. Thus, an informed public should be better able to
decide what quality of water it is willing to pay for.
It is the purpose of this work to present a planning tool
(algorithm) to provide optimal solutions for the complex
choices involved in balancing alternative methods for attaining
water quality goals. As one might suspect, there are a tremen-
dous number of alternatives that would achieve these desired
-------
2.
goals in a body of water. One of the most commonly proposed
solutions is for the polluters to increase their levels of
treatment. This is also one of the most expensive solutions.
Within the framework of a proposed regional water quality
management authority, we are going to investigate this problem
with some powerful tools from non-linear programming and
control theory. Note that for the Delaware estuary, which
we are going to study, an authority such as proposed above
(namely the Delaware River Basin Commission) already exists [2].
It appears that Thomann [3] was the first to recognize
the significance and usefulness of characterizing the dissolved
oxygen changes in a body of water as a linear feedback and
control system. Other models for predicting water quality in
an estuary by solving the basic differential equations subject
to boundary conditions have been presented [4]. However, the
flexibility of a linear control approach appears to be far
superior. This is due in part to the availability of an
existing body of knowledge and in part to the tremendous
amount of detail that can be brought to bear on the problem
with no loss of computational effectiveness.
The concept of a transfer function is basic to every
feedback and control system [5, 6]. Thomann has developed the
transfer functions for a general body of water and in particu-
lar a one-dimensional estuary. In conventional feedback control
system design, lead, lag, or an amplifier is introduced into
the feedforward loop to modify the gain and re-orient the poles
-------
of the open loop transfer function. For the purposes of
water quality control, the important thing is to obtain the
transfer functions of the body of water. Control will be
exercised not by modifying the transfer functions but by
changing the input pattern to the control system.
The implementation of regional water quality control
systems for the major rivers and estuaries of this country
should occur sometime in the next 10 to 30 years. The im-
plementation will be preceded, of course, by planning,
analysis, and design of the system. At this point it is
appropriate to give some guidelines to the governmental or-
ganizations or other agencies which will undertake these
projects.
The most important consideration of the agency in charge
of the project will be the overall approach to the problem.
The methods used must be flexible enough and general enough so
that in analyzing the physical system all relevant inputs are
considered and accurate outputs are produced. Reality should
not be compromised due to inadequate methods. Likewise, in the
choice of control measures, the mathematical model must be of
sufficient generality so that all possible alternatives are
considered. It follows that the mathematical algorithm must
be flexible enough to select from the tremendous number of possi-
ble solutions the "best" one in whatever sense "best" is taken.
At the present time, the approach which seems to adequately
meet all these requirements is a combination of a linear feed-
back and control system and a nonlinear programming problem.
The approach is general, flexible, and comprehensive. Any polluted
-------
4.
body of water can be analyzed and any particular pollutant can be con-
sidered. Thomann gives extensive results for the variation of dis-
solved oxygen, DO, in a one-dimensional estuary. The choice of con-
trol measures includes treatment at polluter, treatment at regional
plants, and by-pass piping taken singularly or in any combination.
A rough flow chart of the type of operations involved in imple-
menting the overall approach is given in Figure j . The remainder of
this section will discuss the details involved in each of these oper-
ations. FIGURE 1
Dynamics of Using the Model
Basic Data
from continuous
monitoring
I
Compute
Transfer
Functions
;
Basic Data
(cost, removals,
etc.)
!
Solve Nonlinear
Programming
Problem
r
Adjust goals,
Change plant
locations, etc.
Acceptable
Solution
-------
5.
The considerations involved in determining the boundaries
of the particular river or estuary are as follows: 1) all
polluted sections of the body of water should be included; and
2) all sections, polluted or unpolluted, into which domestic
or industrial waste is discharged should be included. Generally,
the upstream boundary should be chosen at a location relatively
unpolluted and the downstream boundary at a location where the
recovery from pollution is complete (if such a location exists).
For estuaries these locations might be the head of tide and the
mouth of the bay. In the following discussion, we will consi-
der the operations of Figure 1 that are relevant to an estuary
where the water quality variable of interest is the dissolved
oxygen.
The success of a linear systems approach depends on the
availability of continuously recorded time series data. Water
temperature, stream BOD, and dissolved oxygen should be continu-
ously recorded at a number of monitoring stations along the es-
tuary. This data is used in constructing the forcing functions
and in evaluating the theoretical DO time series calculated from
the transfer function. At the least, 1 year of record should
be available before computing the transfer functions-
To calculate the transfer coefficients, the estuary is
segmented into a finite number of sections. The data required
for each section consists of:
1) rate of flow
2) volume and length of the section
3) the diffusion coefficient
4) reaeration coefficient
-------
Also required are the flows and BOD concentrations for
each of the polluters. This data is also basic data for the
nonlinear programming problem. Specific equations for calcu-
lating the transfer functions of an estuary are given by Tho-
mann.
The input to the nonlinear programming problem is the
transfer functions and certain basic data. The data required
for each polluter consists of:
1) rate of waste flow
2) BOD concentration
3) % removal of BOD
4) cost of removing BOD (linear piecewise costs)
The data required for a regional plant consists of:
1) performance curves (BOD removal efficiencies
for different types of wastes separated and
mixed)
2) cost curves for removal of BOD
Also required are pipeline costs and a preliminary set of DO
goals.
Given a preliminary set of dissolved oxygen goals,
initial computer runs will be made (using rough cost equa-
tions which we will supply) to determine the general pattern
of the solution. This will identify preliminary regional
treatment plant locations and preliminary .pipeline locations.
Surveys can be made in these1 areas to determine accurate pi-
ping costs and small-scale pilot plants can be built to de-
termine the feasibility of mixing wastes, and to establish
treatment efficiencies and costs. As this accurate data be-
comes available, further computer runs can be made which may
-------
change the solution and identify the need for better data in
other areas. In addition., some re-calculation of the transfer
functions may be required. This can possibly occur when large
quantities of flow are shifted to different discharge points.
Thus, we envision the interaction of lab teams, field teams,
planners, and goal setters whose activities will be governed
by the information generated by the program.
The entire process can be considered as a learning device.
Each cycle of analyzing output from a computer run followed by
changing DO goals and regional treatment plant locations, or re-
stricting certain classes of variables to be active or non-ac-
tive yields a certain amount of information about the response
of the system. This information, properly interpreted, gener-
ates further actions and the cycle is repeated. This process
is continued until an acceptable solution is obtained.
Considering the size and complexity of the problem, the
cost of computer time, and assuming a reasonable number of
computer runs, one of the prime goals of this work was to develop
an algorithm that gives efficient solutions in a reasonable
time. This became a formidable undertaking because of the many
undesirable features inherent in our non-linear programming
model of the physical situation. It has over 2000 variables
and over 80 constraints. Some of the first partial derivatives
are discontinuous and the transfer functions have such a wide
range that scaling is extremely difficult.
-------
8,
Most of the techniques to be presented here have not been
available to researchers attempting to' solve large scale water
pollution control problems. This is due to the rapid changes
in fields outside water resources engineering and the limited
contact with these fields by researchers. On this account,
we present considerable detail on our solution techniques,
particularly in the area of non-linear programming. It is,
of course, not necessary to master the mathematics of these
techniques as presented especially in Sections II and IV to
evaluate the realism of the general mathematical model and
the conclusions drawn from the model.
-------
II. MATHEMATICAL DEVELOPMENT
The mathematical development of this section will con-
sist of a proof of convergence of the general non-linear pro-
gramming algorithm followed by a sample problem worked out in
complete detail. A previous paper by the authors [7, pages
36-47] gave a complete discussion of the solution of large
scale linear programming problems. The idea given there was
to reduce the number of calculations by carrying a truncated
tableau and generating elements as needed. This approach is
appropriate for problems where the number of variables is much
larger than the number of constraints, or vice versa. It was
shown that, by carrying only the matrix inverse of the currently
binding constraints, any element (or column) of the basic tab-
leau could be generated from a few simple equations. This,
of course, assumes that the elements of the underlying array
are available by generation or from combinatorial considera-
tions. Thus up-dating the truncated tableau instead of the
complete basic tableau greatly reduces the number of calcula-
tions.
In the discussion that follows, it will be shown that
for every non-linear programming problem we can define an as-
sociated local linear programming problem. In fact most of
the burden of calculation in solving a non-linear problem is
-------
10.
in solving the series of local linear problems as we move to
successive points. Thus we can appreciate the efficiencies
obtained by using the truncated tableau to solve the local
L.P.'s of the non-linear problem.
Development of a non-linear programming algorithm
The algorithm employed for solving the non-linear
programming problems of this paper is a general purpose
algorithm which solves problems of the form:
Subject to: g1(y) £ 0 i = 1,...,m-l
m U)
Min gm(y)
y
where y is a vector in En and g1(y), i = 1,.,.,m, are
differentiable functions. It is a "local", "gradient",
"stepwise" correction descent algorithm. By a "stepwise"
procedure we mean given a point y in the domain of the
functions, a "correction" vector Ay is determined and a
new point y = y + k Ay is used for the successor "step".
It is a "local" method because the correction direction
Ay and its length (determined by the scalar k) are obtained
from the behavior of the system in a "sufficiently" small
neighborhood of the current point y°. It is a "gradient"
technique inasmuch as the gradients of the function g1(y)
are principally used to obtain the correction direction.
-------
11,
The algorithm is a "local-gradient" technique because of
reliance on the approximation of the functions g (y) by
the well known:
Approximation Theorem
Let geC1 in an open region D and let E be a closed
bounded subset of D. Let
Vg(y°) «
°
acr(v)
P
m
be the "gradient" of g at the point y°eE. Then,
g(y° + Ay) - g(y°) + ?g(y°)-Ay + RUy) (2)
where
llm. &L&L „ o
uniformly for y°eE.
The direction of improvement Ay in the stepwise
procedure is obtained from the function approximation in
(2) by estimating the Rx(Ay) and solving the associated
local linear programming problem:
• ^fc * •
Vg (y )-Ay' £ -g1(y°) - k r1 (3)
min vgm(y°)'Ay
The r are an arbitrary set of positive constants
• t
which serve as surrogates for the Rx(Ay), where r1 = 0 for
the strictly linear functions, of course the R^Ay) are
-------
12
truly functions of the direction Ay being sought, and sub-
stituting positive constants causes the solution of the local
linear programming problem to be an approximation. In prac-
• •
tice, the r are taken to be the R (Ay) obtained from the
previous Ay correction vector. The k is used to parametri-
cally adjust the solutions of the local linear programming
problems. In general as k decreases, it becomes easier to
achieve feasibility in the local linear programs but unfor-
tunately the gain that can be assured in the convergence of
the non-linear problem also decreases. As will be shown in
the convergence proof it is necessary to have k ^ s > 0 for
some fixed s to prevent convergence to a non-optimal limit
point.
NATURE OF THE SOLUTION
In order to discuss the convergence of the algorithm
we must develop the concept of an e -solution of the problem.
Given any e > 0, it is first necessary to construct a suitable
"feasible" set. This consists of exhibiting an e-j_ such that
0 < e £ e and defining the e-feasible set P = (y | y € E and
g±(y) £ «T i=l, .. .,m-l). The solutions y* are then char-
acterized as:
Global e-minimums
For all yeF
gm(y*) £ gm(y) + e
-------
13.
or
Local e-minimums
For some 6 > 0 and all yeF such that
jy-y*l < 6
gm(y*) £ g(y) + «
Obviously the concept of an e-solution is necessary
when treating real valued functions. The solutions in general
will be real valued numbers and an answer can be determined only
when the level of accuracy desired is specified. A problem is
considered completely solved when for an arbitrary positive e
an e-solution can be guaranteed in a finite number of steps.
We also require the following well known fundamental
results from the theory of linear programming. A general mixed
linear programming problem can always be stated in both a primal
and dual form.
-------
14,
Primal Form
Subject to the constraints
* aiq yq *= ain; I *• i ^ ml
f aiq ^ * ain'* ^ < i < m (6)
SI
Yq >. 0 nx < q < n
Minimize
2 a y
q mq q
Dual Form
Subject to the constraints
xp ^ 0 m-j^ < p < m
Maximize
2 x a
p P pn
The fundamental inequality of linear programming states:
Given any feasible solution y of the constraints
q
in (6) and any feasible solution x of the con-
P
straints in (7), then
2 a y ^ Z x a (8)
-q mq *q ^ p P pn l '
The solvability or unsolvability of the linear program-
ming problem as well as the relationship between the primal and
dual problems can be summarized in a general theorem.
-------
15
Linear Programming Theorem
Given a linear programming problem, exactly one of the
following conditions holds:
1. There exists a finite value v and feasible vectors y*
and x* of (6) and (7) such that
P
2 a y* = 2 x* a = v (9)
^q p pn
(and by the fundamental inequality (8) , they must then
be optimal feasible vectors) .
2. The constraints for the primal problem are inconsistent,
and the constraints of the dual problem are inconsistent
or the dual extremal function is unbounded.
3. The primal extremal function is unbounded and the con-
straints for the dual problem are inconsistent.
There are three bounds required in our development.
Bound 1
S (-x ) £ B1 (10)
P
The bound on the magnitude of the dual variables encountered
in solving the local linear programming problems has to be
imposed as an additional condition on the problem. This condi-
tion merely acknowledges a limit to our computational ability
to distinguish non-singularity in inverting the basis associated
with the local linear programming problem.
-------
16.
Bound 2
£B (ID
Since E is a bounded set, the vectors y of interest are
bounded. By the triangle inequality
|Ay| = |y - y
and the correction vector is also bounded. Therefore, ve can
always superimpose constraints of the form L1 £ y. <[. U1, which
can be treated algebraically without explicit introduction into
the constraint set. These bounds will always exclude alternative
3, an unbounded solution in solving the local linear programming
problems .
Bound 3
£ B, (i = I,...,!!!) (12)
Transposing (2),
R^Ay) - gi(y) - gi(y°) - vgi(y°)-Ay.
But since g (y)eC', both g (y) and Vg (y ) • Ay are continuous
functions and the difference of continuous functions is a
continuous function. Hence R (Ay) is a continuous function
over a compact set E, and is therefore bounded. Taking the
•
largest bound over the finite set of bounds for the R^tAy),
we obtain an over all bound for the remainder terms.
When the constraints of the primal problem are incon-
sistent in solving a local linear programming problem, there
still exists a set of dual variables to a restricted subproblem
of the original problem. These are exhibited naturally in the
-------
17.
linear programming algorithm employed [see "A Complete Construc-
tive Algorithm for the General Mixed Linear Programming Problem, "
Naval Research Logistics Quarterly. Vol. 12, No. 1, March 1965].
Although not all algorithms naturally exhibit these dual vari-
ables, their existence is easily demonstrated. Let H index a
consistent subset
H = [ijVg^y^-Ay £ (-g^y0) - k r1) for some Ay),
(where choosing Ay = L. the set H contains at least the bounding
constraints of the form L. <1 Ay. <£ u.), and w is an index such that
w , Ov . . w , o, r- w
vg (y ) • Ay ;> -g (y ) - k r .
Solve the subproblem
Subject to
vg1(y°)-Ay £ -g1 (y°) - k r1, ieH
Minimize
?gw(y°)-Ay (13)
In this subproblem only the first alternative can occur
since we know the problem to be consistent and bounded. There-
fore, we v/ill always have a w and a set of dual variables x £ 0
such that
vgw{y°) = 2 x vgp(y°) (14)
peH p
fcince the Ay are of arbitrary sign and the dual variables solve
the dual constraints), and
vgw(y°)-Ay= 2 x (-gp (y°) - k~ rp j (15)
peH p
feince only the first alternative of the Linear Programming
Theorem can occur).
-------
18
For a subproblerr. where w < m, to be consistent at its
conclusion and employing (15), k is bounded above since,
T W
Vgw(y°)-Ay = 2 x (-gp(y°) - k rp) £ -gw(y°) - k r
peH p
or
k£-(gw(y°) + 2 (-x )-gp(y°)/(rw + 2 (~x)rp). (16)
peH peH p
In actual practice the k is parametrically adjusted
downwards to a lower bound as required in the course of the
computation. A bound of
* > e/(B3 + B1-B3) (17)
will be sufficient to insure a gain in the original non-linear
problem. Using (16) and (17) 3 the condition for a violation
of this bound is,
2 <-x J-gfy) ^-g(y) - € (18)
peH p
In general, from (15), after reaching non-linear
feasibility,
Vgm(y°)-Ay = 2 (-xjgp(y°) + k -2 (-xjrp (19)
where
2 (-x )-gp(y°) £ 0
P P
and (20)
2 (-x )rp > 0
P P
Hence, a small k improves the rate of gain but decreases the
admissible step size. The art of constructing an efficient
algorithm requires the balancing of the rate of gain with the
admissible step size.
-------
19.
We can now give the conditions for a stationary point
of the algorithm.
Local Linear Stationary Point
There exists a set of dual variables x such that
P
1) 2 (-x ) > B, (insufficient resolution)
P P
There exists a w and H such that
2) 2 (-x )-gp(y°) ^ -gw(y°) - e (inconsistency)
peH F
There exists a yQeF such that
m-1 .
3) 2 (-x^-g (y ) ^ -e (optimality)
i=l
An interpretation of the consequences of termination of the
algorithm with these stopping criteria will be provided after
a demonstration that they are the only conclusions.
Convergence Theorem
I.t is always possible to choose k_ and k in such
a. manner that in at most a. finite number of steps
iL local linear stationary point is reached.
Proof
Choose k = e/(B3+
In order to establish k, we must dominate the
remainder terms R1 (Ay) and require the following con-
sequence of the Approximation Theorem (2) :
For any y eE, there exists a 51 > 0 such that for
R(k-AY)A' |Ay| < ^]f- , i = l,...,m-l (21)
(22)
-------
20.
Recalling that |&y| <£ B , and taking 6 = min 51,
2 i
then for k = 6/B ,
< k" r1, i = 1, ...,,m-l (24)
(25)
Let
supg = max (g1(y°) )ga (y°) > 0]
i
or (26)
supg = 0 if g1(y°) £ 0, i = l,...,m-l
The feasibility of the constraints is then achieved V*
or maintained when k is choosen to satisfy the
constraints,
+ vg1(y°) •&•* + R1(k-Ay) ^ a
with 0 ^ d < 1, and as close to zero as possible when
supg > 0.
Assume we do not encounter a local linear stationary
point. At the conclusion of the solution of the local
linear programming problem (3) we cannot be inconsistent
without violating stationary condition 2) . Thus,
vgNyVayk <1 [-g1^0) - k r1] -k. (28)
Substituting this into (27) it is sufficient to choose
k such that the stronger condition holds:
g1(y°) + [-g^y0)-^ r1] -k + R1(k-Ay) £ d-supg (29)
or
v r-1 / d supg - (1 - k) •qi(y°} R1 (k . Av)
-k r £ k k
With an appropriate choice of d, the term
-------
21
1°
d-supg - (1 - k).g(y) ^ 0 f or o < k < 1
JC
Choose d = 0 when supg = 0, inasmuch as supg = 0 implies
g1 (y°) <^ 0. Choose d - (1 - k) when supg > 0, and recall"
supg > g1 (y°) .
Now with our choice of k = 6/B and invoking (24) , the
condition (30) is satisfied at every step. But then
when d = (1 - k) = (1 - 5/B2) , the infeasibility supg
must be reduced by at least 6/B • supg at each step.
Given e^> 0 in at most a finite number of steps , supg < e,
and remains so in all subsequent steps.
Further, at the conclusion of the solution of the
local linear programming problems, we must have,
2 (-x )-gp(y°) < - e (31)
P P
or we would violate condition 3 in our definition of
a stationary point. Restating (19) for convenience,
Vgm(y°)-AY = 2 (-x)-gp(y°) +k~Z (-x ) -rp. (32)
P P
Employing these results, we can demonstrate a bounded
reduction in the extremal function as follows:
gm(y) = gm(y°) + vgm(y°)-Ay-k + Rm(k^y)
[by (32)]
= gm(y°) + [2(-x_) -gp(y°) + k-z(-x ) -rp] -k
P P P P
[by (31), our choice of k and k, and (25)]
B^ " 2-d+BjT
-------
22
[by algebraic simplification]
B, ,
m / o, £• 6 r i 1 1
= g (y ) ~~ "0— [ i ~ ~7i 5—T" ~~ o—M 1
2 11
Vi'* D x (33)
Inasmuch as gm(y) is reduced'by at least e-6/2*B2(1+B,)
at each step, we can contradict any bound on g (y) in a
finite number of steps. Hence we conclude that a sta-
tionary point is reached in at most a finite number of
steps.
In order to interpret a stationary point, we need an
additional result. Consider any y = y + Ay that is an
e-feasible solution of the constraints of (1), the original
problem. By the Approximation Theorem
gp(y) = gp(y°) + vgp(y0)-Ay + Rp(Ay) £ ^
or
Ov T^yO\ P/\i / ** A \
y )*Ay £ -g (y ) ~ R (Ay) + e, (34)
Now recall that for any solution x £ 0 of the dual problem
as stated in (14) ,
2 x ?gp(y°) = vgm(y°). (35)
P P
Using these two results, (34) and (35), we can obtain a lower
bound on the extremal function as follows:
gm(y) = gm(y°) + 7gm(y°).Ay + Rm(Ay)
= gm(y°) + [2 x vgp(y°)]-Ay + Rm(Ay)
P p
- gm(y°) + 2 x [vgp(y°)-Ay] + Rm(Ay)
P p
+2 (-x ) gp(y°) +2 (-x)-Rp(Ay)
p p p
e
P
x -e (36)
-------
23.
e
Choose e-i = s~ an<^ assume that condition 1) of a local
1 Bl
linear stationary is not violated then
2 x^ ^ -B, and 2 x 'e, ^ -e .
p P x P P
Using this fact and recalling that global or local convexity
for all functions would imply R (Ay) ^ 0. Formula (36)
would then imply.,
gm(y) :> gm(y°) + s (-x ) gp(y°) -e (37)
p p
When condition 3) in the definition of a stationary point
obtains, we conclude
gm(Y°) £ 9m(y) + 2e for all y, (38)
and therefore by definitions (4) and (5) we have achieved a
global or local e-minimum. The identical argument can be
applied to a subproblem when condition 2) in the definition
of a stationary point obtains, to demonstrate e-inconsistency.
In fact, if the stronger condition
Z (-x ).gp(y°) > -gw(y°) (39)
peH p
occurs at termination, the argument would demonstrate absolute
inconsistency of the constraints. Without convexity it is
theoretically possible to terminate at a "saddlepoint", although
this is virtually impossible as a practical matter. The algor-
ithm can be extended to a second order method by explicit intro-
duction of first order necessary conditions as additional con-
straints. It can also be extended to treat equality constraints,
-------
24,
An Example Problem
As an illustration of some of these ideas consider the
following non- linear programming problem.
Subject to:
y - 9 o
g1 (y) = Y
g (y) =-y - Y2 + 1
o
min g(y) =
(Y± ^ 0)
- 3)
- 2)
Graphically the problem becomes
Y2 FIGURE 2
Graphic Solution of Non-linear Example
2--
-------
25
The optimum occurs at the intersection of the line,
- 3 = 0
connecting the centers of the circles g (y) and g (y) and the
binding constraint circle g (y) = 0.
Eliminating y, from the above linear equation
y, = 3/2 y2
and substituting in g (y) = 0.,
9/4 y2 + y2 = 9
y2 = 6/VU
plugging this value back into the expression for y..
yL = 3/2 • 6/VTJ = 9/
-------
26,
which are handled in the internal logic and not introduced
explicitly. The associated local linear programming problems
are:
s.t. vgi(y°).AY^ -gi(y°) - k r1
min vgm(y )-Ay. (see pg. 11 )
In our problem,
vg1(y) = [2y1, 2y2J
vg2(y) = [-1 , -1]
vg3(y) = [2(y1 - 3) , 2(y2 - 2)]
In general, r is taken as the error term R (Ay) of the pre-
vious step and k = 1. The error terms R (Ay) reflect the length
I AY I of the previous step and may cause the local linear pro-
gramming problems to be inconsistent or prevent a local linear
gain. However, as given on pages 17 and 18 , the Tc can then
be adjusted down to as low as
k~ = e/(B3 + B1-B3)
parametrically in the course of the computation to gain consis-
tency and a local linear gain. More specifically, if an incon-
sistency is encountered in the local linear programming problem
there will exist an index w such that
(-gW(y°)-krw) — ygw(y°)-Ay < 0.
We employ (15) page 17 and calculate Tc to eliminate this vio-
lation or,
(-gW(y°)-krW) - ( Z x (-gP
-------
27
k = -(gw(y°) + 2 .(-x ) -gp (y°) )/(rW+ 2 (-x. )rp)
peH P peH p
-(gW(Y°) + 2 <-x )gP(y°)) < €
peH p
then re-arranging the terms,
2 (-x )-9p(y°) > -gw(y°) - €
peH p
and we terminate with the knowledge the non-linear problem is
inconsistent inasmuch as we have satisfied condition 2) page 19
of a local linear stationary point. Note if
-(gw(y°) + 2 (-x).9p(y0)) ^ €
peH p
then by our definition of the bounds B^ B^ and B^>
k ^(B + B-B
and we guarantee a movement of sufficient magnitude to achieve
convergence. In order to insure a local linear gain, we must
have
vgm(y°)-Ay < 0
Using the equivalent expression given in (19) page 18 , yields
2(-x )gp(y°) + k-Z(-x )rp < 0
P P P P
or,
k < -Z(-x )gp(y°)/2(-x J
P p P p
The dual variables x required for these parametric adjustments
are readily available with our L.P. code.
-------
28,
This parametric adjustment will be illustrated in the course
of solving the present problem.
Suppose we start at the point y = (0,0) and since we
1 2
have no previous step estimates, take r = r =0. The initial
linear programming problem is :
s.t.
mn
The implicit range restrictions on the variables are
0 £ &! £ 5 and 0 £ t&2 1 5
FIGURE 3
Graphic Solution of First L.L.P.
1\ 2 34
V
-------
29,
The optimal solution is Ay = (5,5) . The new point will be of
the form.
y = y +
0
0
5
5
where we must determine k £ 1 to reduce SUPG = 1 or achieve
feasibility. The error term R (Ay) is obtained from
R1(Ay)
by a single additional function evaluation of g (y) inasmuch
as g (y°) is already available and vg (y°) is given by the
local linear programming solution. Carrying this out,
R1(Ay) = g1(y) -
- vg1(y°)-Ay
= 41 (-9) -
Using the quadratic approximation,
0
= 50.
g1(y) =
(vg1(y°) • Ay) -k + R(Ay)-k2
(which happens to be exact in this case) to solve for the
k yields
-9 + 0-k + 50-k2 = 0
k = 3/V50" = .424264
Now R2 (Ay) = 0 since the second constraint is linear and
g2 (y) = g2 (y°) + vg2 (y°) - Ay-k
1 10 -k = 0
at k = .1
-------
30,
Hence, we are feasible for any value of k in the interval
.1 <£ k £ .424264. We take k = .424264 which yields the best
gain in g (y) .
The new point becomes
y° = (2.121321, 2.121321)
and in setting up the next linear programming problem we have,
1(y°) = o. , g2(y°) = -3.242641
g
= so , R2(Ay°) = o
The second linear program becomes
s. t.
4.24264 &1 + 4.24264-Ay2 £ -50
-1.0 Ay-L - 1.0 Ay2 £ 3.24264
min -1.75736 Ayj + .242641- Ay2
The range restrictions are
0 £ y° + Ay £ 5
or
-2.121321 £ Ay £ 2.878679
-2.121321 £ AY2 ^ 2.878679
Now, if we take both Ay^ and Ay2 at their lower bounds, the
first linearized constraint becomes,
4.24264 x (-2.121321) + 4.24264 x (-2.121321)
= -18 > -50
-------
31,
and the problem is clearly inconsistent. In the course of the
solution of the linear programming problem, the initial k = 1
is parametrically adjusted down to k = .275148.
The right hand side of the adjusted linear programming pro-
blem is
-gl(y°) - k
= -0 - k-50 = -13.7574
•0 = 3.24264
FIGURE 4
Graphic Solution of Second L.L.P.
-g2(y°) - k R2Uy°) = -92
-------
32
The two linearized constraints were adjusted until they coin-
cided to give a feasible region. The unrefined optimal solu-
tion is
ty1 = -1.12133
Ay 2 = -2.121321
The computer program at this point makes a post-optimal adjust-
ment of k to further improve the solution of the local linear
program. The optimal solution of a linear program always occurs
at a basic solution. The basis consists of n linearly independ-
ent row vectors or constraints where for theoretical purposes
we include the range constraints of the form
e j Ay
-e Ay
- y
y - L..
at the present optimal solution the basis is:
B =
4.24264
4.24264
-1
This consists of the first constraint and the lower bound con-
straint on Ay 2
optimal solution is
"1
= B" (-g - kr) .
In this particular instance,
B
-1
.235702 1
0 -1
-------
33.
and,
.235702
+ k
-1
2.121321
-50
0
2.121321
-2.121321
-11.7851
+ k
With our current adjusted k = .275148 we obtain,
~-l. 12133
Ay =
-2.121321
Now observe,
3/o
)-&y= (-1.75736,.242641)
2.121321
-2.121321
+k
-11.7851
= -4.24264 + k-20.71066
With our current adjusted k = .275148
vg3(y°)-Ay = 1.455856,
and since this is positive, we can not even obtain or achieve a
local gain. In general, as k—yO, we obtain the best "rate of
gain", but the distance we can move approaches zero when we are
on a boundary.
-------
34.
Distance Moved
i o
Suppose one of the constraining functions g (y ) =0 (on bound-
ary1) and we return to this boundary after movement,
1) g1(y) = vgi(Y°).Ay.k + R1 (kAy) = o
From the local linear program the linear gain is,
2) vgi(y°).Ay = -k r1
Substituting expression 2) into expression 1) and using the
approximation: R1(kAy) = k •R1(Ay) yields,
k-f-k-r1 -f k-R1(Ay)) = 0
and we are at the boundary for k = 0 (no movement) and for,
k = k
" Rx(Ay)
Now we neglect changes in the error term due to direction (with
circles there is no error since they have uniform curvature) .
* •
Since we always take r1 = R1(Ay°) this assumption translates to
r1 m Ri(Av°) = AV°TAV°
Rx(Ay) R1(Ay) AyTAy
and the distance moved becomes,
_ T
_ k (AV° &v°)
- m
(Ay Ay)
Assuming we move this distance, the objective is to minimize
gm(y) = gm(y°)+ vgm(y°) -Ay-k + R(Ay)-k2
-------
35
and employing our assumption that
AyAy
as well as our expression for k yields
_ T
m, * m, o m, o. t k ( Av &v )
g (y) = g (y ) + vg (y )-Ay — "»* — **—*-
AY Ay
T
Av°)2
i
T2 T v/
(Ay Ay) hyo ^yo
gm(y°) + (vc*m (v°} ' AV'^+Rm( Ay0) '*2 ) ' Av° Av°
Ay Ay
Now defining
d. = B-1g and £. = B^r1
we have
Ay = d. + k £
and
T T *
Ay Ay = d. id+2^
Further defining dm = Vgm(y0)-d, pm = Vgm(y°) -p_
vgm(y ) -Ay = vgm(y°) -d + k.vgm(y°) -jg, = dm + kpm
With the additional definition of
p m r = p m + Rm(Ay°)
we can rewrite the expression for gm(y) as:
9m(y) = gm(y°)
T
d d+
-------
36
To obtain the minimum value of g (y), we set the derivative
, m, o.
d g (y ) _
d Tc
and solve for k.
Working through the algebra
-B +7B2 _ AC
k =
where,
A = 2- (c!Tp_) -pmr - dm- (£Tp)
B = (dTd)-pmr
C = (dTd)-dm
In our example,
(;p_Tp_) = 138.88882
(dT£) = -25.
(dTd) = 9.
pmr = 70.710663
dm = -4.24263
A = -2946.286
B - 636.39746
C = -38.183762
k* = .0324352
-------
37,
—*
Using this optimal k yields
AY =
2.121321
-2.121321
+ .0324352
-11.7851
0
=
1.739072
-2.121321
and the adjusted right hand side of the underlying linear
program becomes
-g1^0) - k*R1(Ay°) = 0 - k-50 = -1.62175
-g2(Y°) - k*R2Uy°) = g2(y°) - k-0 = 3.24264
FIGURE 5
Graphic Solution of Second L.L.P.
After Parametric Adjustment
-------
38.
Again using the formula,
we calculate
= g1(y) - g(y°)
= 7.524391
Using the quadratic approximation,
g(y) =
(which is exact in this case) to solve for k yields
0 - 1.621761-k + 7.524391 k2 = 0
k = *
7.524391
o
Now R (Ay) = 0 since the second constraint is linear and
2 , « 2 / o\ 2 / o\ i
g (y) = g (y ) + vg (y )-Ayk
= -3.242641 + .3822479-k = 0
for k £ 8.483084.
Hence, we are feasible for
0 <1 k £ .2155338
and the optimum value of g (y) on this interval occurs for
k = .2155338.
Using this value of k, we calculate
y = y0 + k-
2.121321
2.121321
2.496150
1.664104
.215534
1.739072
-2.121321
-------
39,
the optimum solution. The optimum value is
g3(y) = .3666912.
It should be observed that in making the parametric adjustments
to k, we may transfer basis in the linear programming tableau.
This requires calculating the range of adjustment of k to pre-
serve the current basis and making additional pivots if we ad-
just beyond these limits.
-------
III. DERIVATION OF THE NON-LINEAR PROBLEM
In this section,, we will derive the constraints and
objective function of a comprehensive non-linear program-
ming problem for the regional management of water pollu-
tion control. In order to facilitate the derivation we will
set up a small scale problem. In Figure 6 we have divided
a length of estuary into three sections numbered 1, 2, 3.
We have located five polluters numbered 1,2,3,4,5 along
the length of the estuary discharging waste water into the
section denoted by the arrow. Finally we have fixed the
location of three regional treatment plants numbered 1, 2, 3.
The regional plants can be located from political or engineer-
ing considerations. Note that to consider the location of the
regional treatment plants as variable introduces a combinator-
ial programming problem mapping polluters into regional plants
and regional plants into sections. This combinatorial pro-
gramming problem will be dealt with in subsequent work. In
this study, we take the plant locations as fixed. Finally,
observe that if no waste water is piped to a particular re-
gional plant, then that regional plant is not opened.
Suppose that we have obtained from each polluter the
current flow f^. in (MGD), where f... is the flow from
polluter j to section i, and the current BOD concentration 3.
in Ib/MG of effluent discharged to the estuary (see Table 1 )
-------
41,
FIGURE 6
Physical Layout for Example
Direction
of
Flow
TABLE 1
Flows and Concentrations for Example Estuary
polluter
Numbe r
1
2
3
4
5
Current Flow f . .
(MGD) 3
19.7
7.0
3.0
6.1
4.0
Current BOD J.
(Ib/MG) 1
155.
1801.
666.
278.
557.
In order to compute our piping costs we have obtained the
following distances all measured in miles: (1) d..., the
distance from polluter j to the center of section i (See
Table 2); (2) d.., the distance from polluter j to regional
treatment plant i with » indicating a very large number to
-------
42.
prohibit piping wastes across the estuary (See Table 3); and
o
(3) dr., the distance from regional treatment plant j to sec-
tion i (See Table 4).
TABLE 2
Distance Polluter to Section for Example Estuary
Polluter
1
Section 1
2
3
0
3
5
2
0
2
5
3
2
0
3
4
2
0
3
5
6
3
0
TABLE 3
Distance Polluter to Plant for Example Estuary
Regional 1
Plant
2
Polluter
12345
2 03 1 00 4
00 1 CO 2 00
4 00 3 00 1
TABLE 4
Distance Plant to Section for Example Estuary
Section
Regional Plant
123
1
2
3
1 0.1 5
0.1 2. 2
4. 6. 0.1
-------
43.
Lastly we have a 3 x 3 steady state response matrix of
transfer coefficients a.. in mg/1 per pound per day (mg/1/lb/ctay)
(See Table 5 ). The transfer coefficient is defined as the
change in DO in section i of the estuary due to a change in
input of BOD in section j of the estuary. Thus, in Table 5
an input into section 1 of 100,000 Ibs/day will produce a DO
decrease in section 3 of 0.8421 mg/1.
TABLE 5
Transfer Coefficients for Example Estuary
Section 1
Number
2
3
Section Number
123
1.096 x 10~5 5.328 x 10~6 2.214 x 10~6
1.047 x 10"5 9.817 x 10"6 4.854 x 10"6
8.421 x 10~6 9.431 x 10"6 9.108 x 10~6
The five waste discharges into the estuary will produce
a DO profile down the length of the estuary for steady state
conditions. Now suppose some regulatory agency has decided
to change this profile and has given us a vector of DO
changes, ci, in (mg/i) (See Table 6 ). Thus our problem can
be stated as follows: we must achieve the desired DO changes
in each section at a minimum overall cost.
TABLE 6
DO Changes for Example Estuary
DO Change Vector
mg/1
Section 1
Number
2
+ 0.12
0.00
-0.12
-------
44,
In deriving the mathematical programming formulation, we
will need the following additional notation:
f.. the flow from polluter j to section i in (MGD)
p.. the flow from polluter j to treatment plant i in (MGD)
t.. the flow from treatment plant j to section i in (MGD)
j. the BOD concentration of effluent from polluter i
(Ib/MG)
r. the fractional removal of BOD at treatment plant i
The corresponding matrices (F, Pf T, J, and R) of these five
quantities completely specify the variables of the problem.
Finally we will adopt the convention that a variable with a
bar, such as f.., denotes the current value of this variable.
The mathematical programming formulation of this problem has
three sets of constraints.
Set Is Material Balance at a Polluter
In matrix notation we have
FT-1 + PT-1 = FT-1 + PT-1
which expands to
fll f21 f31
f!2 f22 f32
f!3 f23 f33
f!4 f24 f34
f!5 f25 f35
1
1
1
+
Pll P21 P31
P12 P22 P32
P13 P23 P33
P14 P24 P34
P15 P25 P35
1
1
1
=
f 11 ° °
?12 o o
o f23 o
0 f O
24
o o f35
1
1
1
0
1
1
1
Thus we see that the matrix multiplications give the desired
set of constraints.
-------
45
Set 2: Material Balance at a Regional Plant
In matrix notation we have
P*l = T
which expands to
P21 P22 P23 P24 P25
P31 P32 P33 P34 P35
V
1
1
1
1
=
32
Set 3: Quality Constraints
The third set of constraints, called the quality con-
straints, is more complicated. The quality constraint for
section 1 is:
Change in BOD
input to section 1
L12
113
Change in BOD
input to section 2
— - ••
_~_ ^^~~
Change in BOD
input to section 3
From the first term we need the change in BOD input to section 1
in Ibs/day. There are two components of this change: (1) the
direct discharge (polluter to section), and (2) the indirect
discharge (treatment plant to section). For the first term
the direct effect is
-------
46,
The indirect effect is:
11 W21 W31
t, ,
_ _
12+t22+t32
P25J5)(t
P3555)( .
^^
2 3^3 3
We can combine the direct and indirect effects and simplify
to give:
(1-r.) (p, J)t..
O
— — — — % 1 ~1 • 11 1
cL ^^l-*»^^i-* i £ ^**J*^.^i-i _w • "*~ "J ~™ ™* r * Tvn ™"^^™~^^ •• "^^^"™"
T—
T T
-------
47
The second and third terms of the quality constraint are
(l-r )(g 3)E^ (1-g )(Pl.J)tal
-
f24:i4 ' f2-J
t22
T— T1
*-2 X fc-2
(l-?3) (g3.j)t23 (l-r3)(p3.J)t2
lTt , " lTt . -1
3
f3-J
(l-r2)(p2.J)t32 (l-r2)(p2.J)t32
^-2
(l-r3)(p3.J)t33
,Tr
Now combine the three terms and multiply through by (-1) to
give:
Jjt.., (l-ra)(p2.J)t.2 d-r3)(p3,J)t.
T
-2
)(p J)t '
C;L - AI.FJ -
Cl-r3)(p3.J)t.3-J
-------
48
Generally we will have k treatment plants which gives the
generalized expression:
t.^l-r^ (Pi.J)
.jd-Sj) (P±. J)
i=l
-i
Now we generalize for all three sections to obtain the
desired result:
k
AFJ + A
.u, irt
i
ytji ^ (Pi~'
>, A.,
k
C - AF J - A
!Xt ,
0 ,
The two sets of material balance equations and the set of
quality constraints constitute the explicit constraints of this
problem. For this problem^ the three sets of explicit constraints
and the five blocks of variables can be conveniently represented
in a schematic as follows.
Variables
Constraints X 15
Material Balance
at a Polluter
Material Balance
at a Regional
Plant
Quality
Constraints
8
11
30
39
44
47
ij
o
ij
Pij
Pij
Pij
o
fcij
ij
O
o
*i
©
0
r .
-------
49,
Thus the local linear programming problem has 11 con-
straints and 47 variables in addition to the upper and lower
bound constraints. By eliminating certain blocks of vari-
ables (and, or blocks of constraints), we can consider a
number of different problem combinations. These combinations
are:
(1) treatment at polluter only;
(2) by-pass piping only;
(3) treatment at regional plants only;
(4) treatment at polluter and by-pass piping;
(5) treatment at polluter and treatment at regional
plants;
(6) by-pass piping and treatment at regional plants; and
(7) by-pass piping, treatment at polluter, and
treatment at regional plants.
The combination we are specifically interested in is num-
ber 7. Combinations 1 and 2 can be reduced to linear pro-
gramming problems and have already been solved. The linear
programming problem for treatment at the polluter (case 1)
has been solved by Thomann and Marks [8]. Their formula-
tion can be obtained from the overall problem very easily.
Note from the schematic that the variable j^ appears only in
the quality constraints. The general quality constraints ares
k k
*» . —._ — ..— —
v t (1—r.) (p • J)
AFJ + A L -^ ST — + C - AFJ - A
-------
50,
which reduce to
C (where AJ = J - J)
inasmuch as P = o, T = o, and F = F without treatment plants
or by-pass piping. Recall that j^ the current concentration,
is used to reflect additional treatment at the polluter. Im-
posing a policy that polluters shall not decrease their level
of treatment (or treat negative amounts) leads to the range
restrictions:
0 < j • < U . /f . where j . = U . /r" .
•^ Ji -^ r i Ji r i
The objective functipn is
2 B±W± = 2 s.. (j^ - j/f^) = S (s/f^-aj^
where s. ($/lb/day) is the cost of removing an additional pound
per day. The possibility of multiple price breaks is treated
in detail in Section V.
The linear programming problem for by-pass piping (case 2)
has been solved by Graves, Hatfield, and Whinston {7], For
this formulation the variables are the f . . and the j. are fixed
at their current values j . . The first set of constraints in
the overall problem (material balance at a polluter) reduce to
m _m
F -1 = F -1.
There are no material balance at a regional plant constraints
and the quality constraints reduce to
AFJ - AFJ C.
-------
51,
Let J = J which gives A(F - F)J ^ C, Multiply through by (-1)
and re-arrange to give
AFJ £ -C + AFJo
The objective function is given by
2 c. .f. .
ij n ^
where the piping cost coefficient c.. is given by (c..=1865-
,1 -.402, 1 1D
Vij '•
To complete this section,, we introduce the
concept of a priority class for the variables. Basi-
cally, the idea is to solve the problem with only a
small subset of the variables active (priority class 1)
and the rest set to zero. The problem is then resolved
with the variables in priority classes 1 and 2 active
and the remainder of the variables set to zero. Note
that we can have as many priority classes as we want.
By putting the more important variables in lower pri-
ority classes, computational efficiency is gained as
we step up through the priority classes. This is be-
cause the values of the critical variables will tend
to stabilize and the number of pivots to achieve fea-
sibility and optimality will be reduced. As an exam-
ple, we might want to arrange the variables in the
following five priority classes.
-------
52
TABLE 7
Example of Priority Class
1 (J1, J2, J3, J4, J5, fix, f12, f23J f24> f35^
2 (pi;L, P22, P13, P24, P35, t21, t12, t33)
3 (P12, P14, P15, P21, P23, P25, P3l, P32, P33,
i- t- f- t- f- +• 1
11J 13' 22* 23J 31* 32J
_I_ JL J- ~J fc* fc* £• ^J «J «i. +J £~
4f-F -F -F f -F -F F -F -F
L 22J 13J 14J r!5J t21J 34' 25' T3l> T32>
^ f T" T* T1 T
5 lrl' r2J 3J
Note that the solution to priority class 1 gives us the
solution to case 1 (treatment at the polluter only). Opening
up priority class 2 allows each polluter to pipe to the
nearest regional treatment plant and each regional plant
discharges to the nearest section. Priority class 3 opens up
all other combinations of polluter to regional plant and
regional plant to section. Priority class 4 opens up all
by-pass pipes and priority class 5 allows for variable
treatment efficiencies at the regional plants.
The development of the costs of treatment at the polluters
and the regional plants is given in Section V. In Section VI,
we present the solution to the small scale problem of this
chapter for a number of different problem combinations.
-------
IV. GENERATION OF THE UNDERLYING ARRAY OF THE LOCAL L.P.
In Section II we reviewed the truncated tableau and
generation of elements technique for solving linear program-
ming problems . We required the underlying array to be
available from secondary storage or by generation. The
underlying array of the local L.P. can be conveniently
represented by a schematic as follows:
Material Balance
at a Polluter
Material Balance
at a Regional
Plant
Quality
Constraints
9*13
..
&G:
as!
Note that each column of this array corresponds to a particu-
lar variable. Thus, given any variable, the problem is to
generate the corresponding column of partial derivatives.
Mathematically, we need the matrix derivatives with respect
to the particular variable. Systematic procedures for obtain-
ing matrix derivatives are given by Dwyer and Macphail [ 9].
At this point we present the matrix derivatives that we
will need in the development that follows. Let Y be a matrix
whose elements are functions of the elements x of X and
let A, B, and C be matrices whose elements are not functions
of x . Define K as the matrix having the dimensions of X
-------
54,
with all elements zero except for a unit element in the m
row and the n column. The formulas for the derivatives
are given in the following table.
TABLE 8
Matrix Derivatives
th
(1)
(2)
(3)
(4)
FORM
Y = XB
T
Y = X B
Y = AX
Y = AXC
?>Y
9xm,n
K B
ra,n
T
Kn,n,B
*Vn
*Vnc
The material balance at a polluter constraints are
-1 - PT'l = 0.
G1 = FT-1 + PT-1 -
"^" BG T
The ^j— has the form (2) and we write ^ = Kn m«l where
m,n m,n *
K has the dimension of F. Similarly the
SPn
has the
form (2) and we write
sp
= K •! where K has the
m,n
dimension of P.
-------
55,
The material balance at a regional plant type constraints
are
G2 = P-l - TT-1 = 0.
") 2
The dG /dp has the form (1) and we write dG /9Pm n =
K •! where K has the dimension of P. Similarly, the
dG2/dt has the form (2) and we write dG2/at = -KT •!
m^n m,n n,m
where K has the dimension of T. The quality constraints
can be written
G3 = AFJ + A Y t*i(1"'ri) (pi«J) +
t m
i=l lTt.±
C - AFJ - A YE.i(1-;i)(Pi-5) £ 0.
Block 1:
The dG /df has the form (4) and we write dG3/df
m,n ' I
AK J where K has the dimension of P.
m,n m,n
Block 2:
The dG /dPm n has the form (4) which gives
iTt.m
where K is a row vector with a 1 in the n colxunn. The
m«
row has the same dimension as p> .
-------
56,
Block 3:
For SG /Bt n we have
(lTt.n)2
^T"\
where Te is the negative of t^ except for the m row where
T
= It - t
*n
Block 4
The SG /dj has the form (3) which gives
3 (1-r>
i=l -i
where K is a column vector with a 1 in the m row. The
*m
column has the same dimension as J.
Block 5:
For 5G /Sr we have
m
lTt
•m
At
These rules hold for all cases except for blocks 2f 3, 4, and
5 when the numerator and denominator of the fraction goes to
zero. For this case we need L' Hospital's Rule which can be
stated as follows [10., page 110].
L' Hospital's Rule:
Let f (x) and g (x) represent functions which are
differentiable at each point of a neighborhood N£ (a) .
Assume f (a) = g (a) =0 and g' (x) ^0 for all xeN (a) .
-------
57
If
lim f ' (x)
x-«a g1 (x)
exists, then
lim f (x)
g(x)
also exists and
lim f (x) _ lim f (x)
g(x) x-*a g1 (x)
We will need the following relationship
T
Pi-1
Block 2
3 (1-r )
Now consider
r (1"ri)
~ L ~^~ [At-ijj
i L 1 t 1 3
Let f(tk±) =
and g(tki) = tk±.
Note that when t^- 0, f (0) = g(0) = 0,
Now apply L'Hospital's Rule to give
f' ft ) = (1-r.)i.a i
t ^ki; ^ ^a.'-^ -k
g1(tki) =1
and
_ n-r- ^-i a
- (1~r) a
7) ~ vx • i'Jj -k
-------
Block 3:
58.
Now consider
and
3P
..) (p...J)A(-t...)
(1Tt-j) (1Tt>j)J
.. l-r.(p..J)a.
L P.l J
JK
Taking (B) first, let
f(PjK) = (l-r^
and g(pjK) = PjR.
Then
and
f • »w / **V **, «4 ' T
•— I r> 1 T J
-Kj L Pj.X lTt.^ J
Let PjK = fcKj =
e snce
lim
g(€)
-------
59,
We have
lira
= (i-r H
( rJ
and
lira
"
Combining (A) and (B) gives
Block 4:
^gi ^ y d-r±)
a^m *m jT-^ 1 t.j_
Now consider
Here we have to consider each term separately.
Let f(tji) = d-r^P^a.jt^
and g(t.^ = t^.
Apply L1 Hospital's Rule to give
= 0
9' (tj±) = 1
-------
60
Block 5:
B r T
ri It...
Now consider
r- Pj.j*t.i-.
. L J
Let f (t^) = -
-------
where
N is the variable number
N is the number of polluters
Thus,
i =
4 ~
+1=1 and j =4- [0x5] =4
fll
f21
f31
f!2
f22
f32
f!3
f23
f33
f!4
f24
f34
f!5
f25
f35
Applying the rules we have
61.
000
000
000
100
000
^14J
a21
a31
1
1
1
H «•
=
a22 a23
S32 a33
all J4
a2l j
a31 j
4.
4
0
0
0
1
0
= o =
oooio
00000
0 0 0 00
0
0
0
the completed column is
-------
62.
0
0
0
1
0
0
0
0
21
where k,, is the cost coefficient. In this example note that
in calculating the dG /3f14 we multiplied the A matrix
against K,4J. Essentially we are taking the inner product
of each row of the A matrix against the column vector K-.J.
In fact, every block of the quality constraints can be
reduced to a multiplication of the A matrix against a column
vector! When we consider the size of the underlying array
of the Delaware Estuary problem, and the number of columns
generated to solve one local L.P., we can get some idea of
the tremendous number of inner products computed in solving a
problem of this size. A scheme for computing inner products
more efficiently has been given by Winograd [11].
-------
63,
Let X = (xx, x2, x3, x4, x5, xg)
and let Y = [y^ y2, y3, y4, yg, y6l .
The standard computation of X-Y is
X-Y = xiyi + x2y2 + x3y3 + x4y4 + x5y5 + x6y6
which has six multiplications and 5 additions. Now suppose
we calculate X-Y as follows
X-Y = (xx + y2) (x2 + yx) + (x3 + y4> (x4 + y3) +
(x5 + y6) (x6 + y5) - (Xlx2 + x3x4 + x5x6) -
(y!Y2 + y3Y4 + y5y6)
At first glance it appears that we have lost ground. But
if the last two terms have been pre-calculated and the
results stored, we have only three multiplications and ten
additions. In terms of machine computation, we have gained
ground because we halved the multiplications while doubling
the additions. The formulas used in this procedure are as
follows.
Let X = (xif...txn)j Y= (yx, ...,yn)
n/2 n/2
X2j-l X2j 1 = E Y2j-l y2j
n/2
Z
3=1
j} (X2j
-------
64,
The inner product X-Y is then given by
Y - § - 1] if n is even
-^-T] + xy ifnis odd
* ' n n
We make use of this procedure by noting that in the expres-
sions for the partial derivatives of the quality constraints
we can factor the A matrix out of each expression. Thus we
pre-calculate all the row products for the A matrix and store
the sum. The remainder of the expression is reduced to a
column vector and the column products are calculated and
the sum is saved. Then we apply the above formulas to
get the inner product.
We can also make use of this inner product routine in
evaluating the quality constraints, G . Recall from
Section II that in setting up the local L.P.j we evaluate
G (y°) for the right hand sides. The quality constraints
are ,
3 * t (l-ri)(p J)
G = AFJ + A ) ~ AFJ + C £ 0
i=i * fc-i
Factoring out A gives
3 r £ t.i(1-ri)(Pi.J) - T
G = A[ FJ + 2, —-—|—i FJ J + c <: o.
i=l 1 t-i
Note that the expression in brackets reduces to a column
vector.
-------
V. DEVELOPMENT OF TREATMENT COSTS
Treatment Costs at the Polluter:
The costs to polluters to reduce BOD discharges
were determined by the Delaware Estuary. Comprehensive Study
Group through direct questioning of the officials of the
industrial firms and municipal treatment plants. They
were asked to determine how much it would cost to reduce
BOD discharges by the year 1975 [12, page 50].
The data collected included present discharges and
the cost to reduce BOD discharges (nearly) to zero. Thus,
each polluter provided either one, two, or three points on
his waste reduction cost curve. That is, each polluter
gave the cost of reducing waste discharges by one, two,
or three different amounts.
The cost estimates provided by the polluters are the
best cost estimates available, yet it _s important to recog-
nize their limitations. Engineering estimates of costs
cannot be made with complete confidence, and as a result
there may be significant errors in the cost data. It may
also be the case that some of the polluters biased their
estimates in an effort to protect their interests. Since
process changes were considered as an alternative to waste
treatment in the estimation of costs, the inherent diffi-
culty of estimating the cost of process changes is involved,
The data we used in this study is given in Table 9 .
-------
66,
Polluter
TABLE 9
Cost and Discharge Data
Polluter
Number
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Slope s .
$/lb/day
5980
3769
2455
1843
10702
1749
105
4809
59
1051
191
2735
33
192
66
91
149
1452
199
285
14
498
1836
1471
250
4157
Increment
Bound B .
Ib/day
2040
1133
1833
346
228
2093
1333
445
5808
1161
892
892
50100
35070
42507
8501
9712
1942
130854
18319
12238
2448
1780
214
880
587
Flow f.
MGD
19.7
6.1
250.
15.2
4.6
3.0
2.0
4.0
144.7
22.4
7.0
107.6
3.5
10.6
4.9
Present _
Discharge U .
Ib/day
3060
1700
2750
990
3140
2000
7550
2230
125250
59510
12605
170110
15910
3560
1760
-------
67
TABLE 9 CON'T.
Polluter
Number
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
Slope s.
$/lb/day
247
258
3416
559
6230
135
4377
232
291
123
1828
346
724
6488
177
324
324
198
1108
3235
3474
3403
126
3630
348
Increment
Bound B.
Ib/day
20672
9923
1985
3346
305
1552
345
117150
18223
1335
267
9471
8287
1973
5091
1389
1389
27160
4000
2000
2107
717
1350
270
5365
987 3745
2317 2000
3791 1400
252 7275
4192 1455
Flow f .
i
MGD
28.2
51.0
6.8
2.0
118.0
0.7
38.6
0.3
80.4
11.6
4.0
0.7
6.6
14.8
8.6
Present _
Discharge U.
Ib/day
21925
12900
3955
2070
156200
1870
25650
8100
34160
3160
1075
1890
7750
7745
10185
-------
TABLE 9 CON'T.
68,
Polluter
Number
31
32
33
34
35
36
37
38
39
40
41
42
43
44
Slope s.
$/lb/day
263
872
2912
487
107
3286
198
410
796
384
858
2900
803
7197
171
359
140
3067
3040
5055
469
686
381
107
2340
458
15844
116
2535
Increment
Bound B .
i
Ib/day
1219
1273
831
2855
1400
280
19590
2857
3265
8044
16089
3448
1445
289
53732
8059
1073
238
2500
2750
74000
30500
5066
1402
312
1923
385
1298
288
Flow f .
i
MGD
112.2
5.4
1.4
1.1
106.9
33.3
60.0
1.0
4.8
103.3
10.1
1.2
278.6
1.3
Present
Discharge U.
Ib/day
3600
3145
1820
32650
28730
2890
85970
1430
8480
110000
6755
1870
2500
1730
-------
69,
As an example of the meaning of this data, consider
the polluter number 23 which supplied three points on the
cost curve. This polluter is currently discharging 8100
Ibs/day of BOD. For the meaning of the slope and upper
bound, see Figure 7. The costs in Figure 7 are the pre-
sent value of the construction and operating costs. They
were computed using a present value factor of 13 by the
following formula:
(Total cost) = cap. cost + 13(op. cost) 1
To make these costs compatible with our piping costs which
are annual costs, we make the following conversion:
(Total
= cap.^cost + op. cost
FIGURE 7
Piecewise Linear Cost Curve
Present
Value
Cost^
$xlO
5091 6480 7869
Ibs/day removed
See [12] pg. 52. The present value factor of 13 is based on
a discounting period of 20 years with a discount rate between
3% and 5% per annum.
-------
70.
For our purposes, the above data is not in a desir-
able form . We need an expression for the cost as a func-
tion of the variable j. Initially, we used the general ex-
pression,
K = c-^UB-j) 2 where K = cost,
to fit the data points supplied by that polluter (See
Figure 8).
FIGURE 8
Smooth Curve Fit of Cost Data
Annual
Cost K
$x!06
LB
10"3lb
MG
However, this expression was found to be undesirable be-
cause the fitting introduced local minimums near the upper
bound. Thus, most of the polluters moved slightly off their
upper bounds and a large amount of calculating time was wasted
making slight improvements which were meaningless.
-------
71,
The most satisfactory cost function was found to be of
the simple piecewise linear type. No external constraints
were added and all the algebra and logic were done internally
in the machinery of the algorithm.
Using a piecewise linear cost function, we define for
each polluter with three price breaks:
U - U - Bl -
Jl = - x 10 J, J2 = x 10 J
f f
U - Bl - B2 .. U - Bl - B2 - B3
J3 - _ x 10"-% J4 = x 10~J
x 10~3/13, k2=s2f x 10~3/13, k3=s3f x 10~3/13
If we have two breaks we simply omit J4 and k_ and if we have
one price break we omit J4, J3, k_ and k2. The cost function
becomes
K = (Jl - j)•kI J2 < j £ Jl
K « (Jl - J2) -kx + (J2 - j) -k2 J3 < j £ J2
K = (Jl-J2)-k1 + (J2-J3)-k2 + (J3-j)-k3 J4 ^ j £ J3
The cases with one, two and three price breaks are plotted
in Figures 9, 10, 11.
-------
72
Annual
Cost K
$x!06
FIGURE 9
One Price Break
1 Bound
10 3lb
MG
Annual
Cost K
$x!06
FIGURE 10
Two Price Breaks
2 Bounds
10 3lb
MG
Annual
Cost K
$x!06
FIGURE 11
Three Price Breaks
3 Bounds
10"3lb
MG
-------
73,
The partial derivatives of the cost function for each
of the three cases are as follows:
Case 1: J2 £ j. £J1 = -s^F - (lo"3/13)
-I- o J -• J- J-
Case 2: J3 £ j, £ J2 = ~s?f. ' (10"3/13)
j- °-'
Case 3: J4 £ j . £ J3 7 = "s?F-i ' (10~3/13)
-•• J-
Cost of treatment at a regional plant:
The variables corresponding to regional treatment plants
are the r^ (the fractional removal of BOD at the ith regional
plant). The costs of BOD removal as a function of the percent
removed have been given by [13, page 10]. Costs were given
in terms of total annual costs of sewage treatment per MGD
capacity based on 20 year, 4-|% municipal bonds covering the
total initial capital investment plus the annual operation
and maintenance costs of the treatment plant. Costs are
representative of plants designed and operated in humid areas
where water scarcity is not a problem. These costs are
given in Figure 12 with the design capacity as a parameter
for four different sized plants. The costs were fitted to
an equation of the form
K = a(w - .5)3 + b
-------
74,
The equation was parameterized to yield the following general
expression
49.22
Kc =8-°W - -53 + l (40)
where
K is the total annual cost per MGD in $1000
w is the fractional removal of BOD
Q is the design capacity in MGD.
As an example of the economies of scale obtained by
using regional treatment, consider the following situation.
Suppose four polluters, each having a flow of exactly 2.5MGD
decided to build treatment plants to obtain 50% removal.
From Figure 12 we see that the total annual cost would be
38, 000 x 2.5 MGD x 4 = $380,000
If the four polluters were located close enough so that
piping costs would not be excessive, they might consider build-
ing only one plant and sharing the cost. From Figure 12, the
total annual cost would be
i x 10 MGD = $260,000
which is a $120,000 savings.
Theoretically, the cost curves go to infinity as w-»1.0.
To cope with this physical reality and its mathematical
consequences, we restrict w ^ .98. The bound can be chosen
to reflect the accuracy of the cost data. Before presenting
the scheme for computing costs at the regional plants, it
-------
80
70
60
o
o
2 50
4*
40
oc
UJ
Q.
(0
O
O
30
< 20
z
z
<
< «0
I-
o
I-
% BOD REMOVAL
FIGURE 12
Plant Costs
10. MGD
2.5 MGD
Ke»304(V-.3)3* 38
Ke «2I6(W- .5) •»• 27
50 MGD ,
Ke ' 152(W- .5) * 19
100 MGD
K, =I28(W- .5)3+ 16
-J
cn
-------
76
might be interesting to compare annual costs at the polluters
with those at the treatment plants. Notice that equation (40)
is an efficient tool for making this comparison.
At polluter 1 we have:
f-j^ = 19.7 MGD BI = 2040 Ibs/day
j-L = 18.6 mg/l s^ = 5980 $/lb/day
U, = 3060 Ibs/day
qj_ = -85 (fractional removal at polluter 1)
If we removed all 2040 Ib/day, the cost would be
•^H^ x 204° = $938,400
with U-L = 3060 - 2040 = 1020 Ib/day
^ - 1020 ' _ ._
Dl ~ 19.7-8.33 = 6'2 mg/1
and
- i _ 6.2 _ Q_
ql •"• 124 " 'y5
~ 18 6
Note that j-^ = ^1_ Q^^ = 124 mg/l (untreated concentration)
Now evaluate (40) with Q = 19.7 and w = .85 and .95.
w = .85
K~ =- 49'i/4[8.0(.85-.5)3 + l] = 31.376
C 19,
and the cost is 31,376-19.7 = $618,111.
w = .95
49 221— "} —i
Kr, = '~TTK 8.0(.95-.5) J + 1 = 40.394
19.7'
and the cost is 40,394*19.7 = $795,767.
-------
77.
The first cost $618,111 is the present annual cost to
maintain treatment at the 85% level. The second cost $795,767
is the annual cost of constructing and operating a new plant
that would give 95% removal. The estimate given by polluter 1
is $938,400. We conlude from the magnitude of this figure
and its closeness to the theoretical value fo $795,767 that this
cost is not for additional units. It must be the cost of a
complete new plant to give 95% removal.
This same type of comparison was made for all 44 polluters
and the results are given in Table 10. Column two gives the
classification of the type of waste, i.e., domestic or indus-
trial. Column three gives the untreated concentration in
(mg/1). Columns four, five and six give the linear costs ver-
sus the costs computed from (40) for different levels of per
cent removal at the polluter. Note that the comparison is made
for the 22 industrial wastes although equation (40) is not valid
for plants treating industrial wastes. The conclusion from
Table 10 is that the costs are mixed, i.e., some costs are for
new plants. Thus the cost data are inconsistent in the sense
that two identical polluters can have widely different costs to
achieve the same increase in treatment.
Note here that we are using a translated set of w's in
the sense that the fractional removal depends on the point of
view. In the above example, we could say that since the treat-
ment level was increased from 85% to 95% the increase in treat-
ment was 10%. However, we didn't know this fact we
say that since 2040 of the 3060 Ibs/day were removed.
-------
78,
TABLE 10
Comparisons of Estimated Plant Cost
POLLUTER
NUMBER
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
TYPE
WASTE
D
D
I
I
I
D
I
D
D
D
D
D
I
I
D
I
I
mg/1
124.
223".
9.
8.
546.
178.
697.
268.
416.
456.
333.
292.
840.
40.
144.
267.
47.
PER CENT
BOD REMOVAL
0.85
0.95
0.85
0.95
0.85
0.95
0.0
0.35
0.58
0.85
0.95
0.55
0.85
0.95
0.35
0.85
0.95
0.75
0.85
0.95
0.75
0.85
0.92
0.30
0.80
0.90
0.35
0.85
0.95
0.35
0.85
0.92
0.35
0.85
0.95
0.0
0.50
0.56
0.70
0.85
0.95
0.65
0.98
0.35
0.85
0.95
COST
EQN. 1
618111.
795767.
256575.
330292.
4155968.
5350188.
0.
368565.
380439.
207628.
267291.
112309.
150656.
194016.
80543.
111177.
143115.
156617
186966.
240703.
2310178.
2757838.
3270604.
474354.
616253.
766251.
206100.
284487.
366243.
1599974.
2208395.
2618988.
122548.
169142.
217761.
0.
289148.
289651.
172476.
217703.
280319.
618584.
1135208.
913970.
1261511.
1624174.
COST EST.
OF POLLUTER
938400.
328483.
346155.
49052.
236749.
281589.
10767.
175382.
26359.
120222.
13106.
200769.
127177 .
645134.
215805.
275312.
111314.
328221.
2003072.
2404680.
13179.
106957.
251391.
275606.
16923.
204628.
392768.
196933.
718530.
-------
TABLE 10 CON'T.
79,
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
I
D
D
D
I
I
I
D
D
I
D
I
D
I
I
D
I
107.
207.
265.
458.
123.
9261.
78.
218.
215.
463.
217.
157.
203.
6.
108.
240.
4454.
0.35
0.90
0.95
0.40
0.85
0.95
0.40
0.85
0.92
201667.
313312.
358397.
82116.
111135.
143878.
290043.
16117.
143065. i 132276.
1748091. !
2366618. 2090676.
2806635. 2498590.
0.30 35257.
0.80
0.90
0.35
0.59
0.80
0.85
0.65
0.87
0.93
0.99
0.35
0.87
0.94
0.98
0.85
0.95
0.85
0.95
0.30
0.80
0.90
0.35
0.80
0.60
0.79
0.90
0.97
0.30
0.80
0.90
0.35
0.57
0.80
0.95
0.35
0.94
0.35
0.85
0.95
0.20
0.68
0.75
0.83
45782. 12631.
56907. ! 50175.
741643.
766668. ; 252074.
926876. 713596.
1023678. 1698275.
20491.
28036. 69316.
32642. 103934.
38732. 138552.
1285866.
1843309. 413668.
2240168. 754591.
2497884. 1252283.
415490. I
534932. 563055.
186966.
240734. 187689.
35257.
45803 13085.
56953. 88477.
197202. !
246438. ! 143617.
374367.
446450. 284332.
556894. 640793.
677931.
231361.
300572.
373738.
1651005.
1701495.
2063128.
2933687.
169648.
293226.
61638.
85077.
109530.
41448.
55334.
59476.
68067.
1049054.
141023.
610205.
24661.
110050.
296194.
106953.
11523.
82298.
298371.
388476.
588395.
-------
TABLE 10 CON'T.
80.
35
36
37
38
39
40
41
42
43
44
I
I
D
D
D
I
I
D
I
D
129.
10.
215.
286.
212.
197.
401.
312.
2.
266.
0.75
0.82
0.96
0.99
0.0
0.50
0.60
0.20
0.70
0.78
0.40
0.85
0.95
0.0
0.29
0.62
0.35
0.79
0.97
0.80
0.95
0.40
0.85
0.95
0.35
0.85
0.95
0.40
0.85
0.95
1840888.
2065289.
2910532.
3176476.
0.
682299.
687757.
831899.
1129011.
1237634.
48826.
66133.
85118.
0.
148583.
161772.
1551777.
1897312.
2898458.
339091.
482134.
55981.
75762.
97557.
3265802.
4507479.
5804563.
59444.
80508.
103625.
237607.
1299481.
2068650.
89257.
249251.
706782.
929335.
11555.
67705.
584615.
1653941.
2669692.
4279152.
148473.
11540.
67700.
67749.
536975.
11582.
67742.
-------
81.
the treatment increase was 66%. The relation between the
translated variables and the r variables is
w. - w.
r. = -^ ^ (41)
1 1 - w.
To compute costs at a regional, plant, we will make two evalu-
ations of equation (40) and take the difference as the cost.
At the ith regional plant, we note that Q = Zp.. in (40) and
_ 3
we calculate a v. and a v. as follows
51 T"> ( ~1 — ~1 1
and v. £ v. £ .98. Thus, we have
V> ~ -1 ~ D 5 (42)
If Zp. . =0, (42) is indeterminate. In this case, the value
j *•!
of v. is irrelevant since the cost equation will automati-
cally go to zero. From (41) we have
vi = ri + ^i " ri^i
The cost equation is now written as
k
Cost($xl06) = Y .04922(Zp. .)3/4[ (8.0(v.-.5)3 +
l-> -i -LJ J-
- (8.0(V;L - .5)3 + 1)]
-------
82
which can be simplified to
k
Cost($xl0) = .39376(Zp. .
(44)
with v. and v. given by (42) and (43).
Computing costs in this manner assumes that at a regional
plant we build only the additional units that we need to achi-
eve the desired per cent removal. With the economy of scale
effect from combining flows and using only the cost of addi-
tional units , the costs at the regional plants should be com-
petitive with those at the polluters.
The respective partial derivatives of the cost function
are as follows :
= 1.18128 (Spj,)^4 (v^.5)2 (1-v^
?- = 1.18128
(Zp. .)
13
3/4. Pik
2
L-ri>
1.18128 (Zp. .)
Vj-T1
~1/4
OSPij)
j
3/4
-------
83.
If site costs are known for each of the regional plants,
we can include the relative differences of these costs in (44)
as follows:
k
Cost($xl06) = Y . 39376 (Ep..)3/4[ (v.-.5)3 - (v.-.5)3]-S
L~J • J- J -L _L
i=l D
where S is greater or less than one depending upon the ratio
of site cost to some standard site cost. This effect can
differentiate between two plants located extremely close
together.
By-pass piping costs:
The piping costs associated with each flow variable are
given by the general expression [7 } page 23]
MGD - yr
> = 1865
where
Q. . is the flow in MGD and corresponds to a particu-
f±j, P±r or
The cost function for this problem is made up of piping costs
and treatment costs at the polluters and at the regional plants,
The entire cost function can now be assembled and is given as
follows.
.2. "xij + -39376 f
piecewise linear costs at polluters
-(Vi-.5) ]
-------
VI. A SMALL SCALE PROBLEM
Recall that in Section III, we set up a small scale pro-
blem with 5 polluters, 3 sections, and 3 potential regional
treatment plants (See Figure 6 and Table 12). Recall that
the j variables are bounded above and below (See Section V).
The r variables are absolutely bounded, of course, by 0.0 and
1. However, in order to control the magnitudes of the v. (See
Section V) and maintain stability of the cost function, we will
restrict r so that 0 £ r. £ 0.70. Note that from (41) in the
previous section, the real removal rate would be
w. = r. • (1 - w.) + w. .
When we have 95% removal at a polluter followed by 70% removal
at the regional plant, the total removal would be
w. = .70- (1 - .95) + .95 = .985.
The first configuration that we will investigate has the
variables arranged in the following priority classes:
TABLE 1 1
Priority Class
Configuration 1
!2' 23' 24> 35' 11J 13'
P31' P333 P35' fc!2J fc21' t33' rlj
2 (jx, J2, J3, J4, J5)
3 ^fcllj t!3' fc22' fc23' t31J> t32^
4 tf!3' f!4> f!5' f21' f22' f25' f31J
5 (P, P, P, P
-------
85
TABLE 12
Polluter Data for Example
Polluter
Number
1
2
3
4
5
Number
(from Table 9)
1
11
6
2
8
fij
(MGD)
19.7
7.0
3.0
6.1
4.0
Ib/Ac
155
1801
666
278
557
% removal
BOD
85
35
55
85
75
Note that in priority class !_, treatment is restricted
to the regional plants. Activating class 2 forces us to
choose between treatment at the polluter and treatment at
the regional plants. Class 3 allows the regional plants to
choose discharge sections and class 4 opens up all the by-
pass pipes.
From Table 6 (See Section III), we see that there is a sub-
stantial negative goal in section 3. Thus, we might guess
that when the by-pass pipes are opened up, there will be a
dramatic shift in the solution. However, it should be re-
called our cost function for the regional plants and by-
pass pipes is non-convex and we are using a local method.
This confines us to local optima and makes our solutions de-
pendent on the order in which we introduce the variables.
This is illus-trated by the slight discrepancy in the solu-
tions to Configuration 1 and 2 Priority Class 3 which we
would expect to be identical.
-------
86
The solution to priority class 1 is given in Figure 13
and Table 13.
FIGURE 13
Solution of Example Configuration 1 and Priority Class 1
Priority class / 1
Cost = $93.738
*MGD
TABLE 13
Solution Data for Example Configuration 1
and Priority Class 1
regional 1
plant
2
% removal BOD
0
69.4
8.69
Note from Table 13 that regional plant 1 is not opened and
is nothing more than a junction point. Opening up priority
class/y2 gives no further gains.
-------
87
The solution to priority class 3 is given in Figure 14
and Table 14.
FIGURE 14
Solution of Example Configuration 1 and Priority Class 3
Priority class / 3
©
\
.031
Cost = $51,369
TABLE 14
Solution Data for Example Configuration 1 and Priority Class 3
regional
plant
% removal
BOD
1 0
2 36.9
3 —
Polluter
Numbe r
1
2
3
4
5
Lower
Bound
Ib/MG
52
136
74
93
111
Concen-
tration
Ib/MG
155
1801
666
278
557
Upper
Bound
Ib/MG
155
1801
666
278
557
% removal
BOD
85
35
55
85
75
-------
88
Notice that in this solution plant 1 does not treat and
becomes a junction point for by-pass piping. This occurs be-
cause the negative goal in section 3 makes by-pass piping
attractive and the by-pass flow variables remain in priority
class 4. The best solution opening all priority classes was
$44,997.
The second configuration we investigated had the variables
arranged in the following priority classes:
TABLE 15
Priority Class
Configuration 2
13, P15
P14,
The solution to priority class 1 is given in Table 16.
Notice that this solution is pure treatment at the polluter
with a cost of $180,843. Opening up class 2, of course, had
no effect. The solution to class 3 is given in Figure 15 and
Table 17.
-------
89.
TABLE 16
Solution Data for Example Configuration 2 and Priority Class 1
Polluter
Number
1
2
3
4
5
Lower
Bound
Ib/MG
52
136
74
93
111
Concen-
tration
Ib/MG
155
358
222
278
334
Upper
Bound
Ib/MG
155
1801
666
278
557
% removal
BOD
85
87
85
85
85
FIGURE 15 Cost = $180>843
Solution of Example Configuration 2 and Priority Class 3
Priority class # 3
Cost = $49,416
TABLE 17
Solution Data for Example Configuration 2 and Priority Class 3
regional 1
plant
2
3
% removal Polluter
BOD Number
— 1
38.9 2
— 3
4
5
Lower
Bound
Ib/MG
52
136
74
93
111
Concen-
tration
Ib/MG
155
1801
666
278
557
Upper
Bound
Ib/MG
155
1801
666
278
557
%removal
BOD
85
35
55
85
75
<
-------
90.
Notice by comparing Tables 14 and 17 that in this solution
we have a slightly higher treatment level at plant 2 than in the
solution to Configuration 1 Priority Class 3. This enables pollu-
ters 1 and 4 to discharge their total flow directly into the near-
est section and gives a slightly improved cost. The best solu-
tion opening all priority classes was $45,406.
The final configuration we investigated was the following:
TABLE 18
Priority Class
Configuration 3
1
2
3
4
*fll' f!2' f!3' f!4' f!5' f21' f22' f23' f24' f25'
f31' f32' f33' f34' f 35 *
^1' J2' J3' J4' J5^
fpll' P13' P15' P22' P24' P31' P33' P35' Hi' fc!2 '
fc!3' fc21' fc22' fc23' fc31' fc32 ' fc33' rl' r2 ' r3^
{P12' P14' P21' P23' P25' P32' P34]
The solution to class 1 is given in Figure 16 and Table 19
FIGURE 16
Solution of Example Configuration 3 and Priority Class 1
Priority class
Cost = $47,119
-------
91.
TABLE
Solution Data for Example Configuration 3 and Priority Class 1
Polluter Lower
Number Bound
Ib/MG
1 52
2 136
3 74
4 93
5 111
Concen-
tration
Ib/MG
155
1801
666
278
557
Upper
Bound
Ib/MG
155
1801
666
278
557
% removal
BOD
85
35
55
85
75
This is the pure by-pass piping solution. The annual cost
of $47,119 is considerably better than the annual cost of
$180,843 for pure treatment at the polluter given in Table 16
and the annual cost of $93,738 for pure treatment at regional
plants given in Figure 13 and Table 13. These comparisons
however are misleading for general full scale problems. The
very short distances obtainable with three sections and the
negative goal in section 3 account for the distortion. in
general all three treatment modes can be utilized efficiently
to contribute to the optimal solution. When seeking the gen-
eral mixed treatment mode optimum, it is usually best to opti-
mize first in Configuration 1 Priority Class 1. This forces
the regional plants to operate at high enough levels to achi-
eve their economies of scale before introducing other cost
comparisons.
-------
VII. COST ANALYSIS OF WATER POLLUTION CONTROL SCHEMES
FOR THE DELAWARE ESTUARY
The Delaware Estuary is bordered by a large metropolitan
and industrial complex, including one of the largest oil re-
fining and chemical areas in the United States. Hence a num-
ber of substantial and complex wastes enter the estuary along
the major part of its 84 mile length. Figure 17 shows the es-
tuary segmented into 30 sections of 10,000 feet or 20,000 feet
in length. Also, there are located 44 major waste dischargers
and their discharge section along with 9 potential regional
treatment plant sites. Note that the influences of the tide
extend to Trenton, New Jersey.
A report by Schaumburg [12 ] gives a detailed comparison
of the total costs of different water pollution control schemes
for the Delaware Estuary. He took the average dissolved oxygen
concentration in the Delaware Estuary in the summer of 1964 as
a baseline. Thus given a water quality standard (say 3 mg/l)
we calculate a set of D.O. goals by subtraction. A standard of
3 mg/1 means that the D.O. in each of the 30 sections must be
at least 3 mg/l. Since he gives extensive results for this stan-
dard, we will use it too. The corresponding D.O. goals are given
in Table 20.
-------
FIGURE 17
Map of Delaware Estuary
93
O = Major Waste Discharger
A = Potential Regional Plant Site
-------
94.
TABLE 20
D.O. standard of 3 mg/1
Section Goal
i _
2
^ «
4
5
6
7
8
9
10 0.8
Section
11
12
13
14
15
16
17
18
19
20
Goal
0.7
1.5
1.8
1.9
1.9
2.0
2.0
1.8
1.6
0.8
Section Goal
21 0.1
22
23
24
25
26
27
28
29
30
In the next few paragraphs we will give a short description
of each of the schemes he considered along with their costs;
1) Required secondary: no specific standard is given
and all polluters are required to use a secondary
treatment process. Those already using a secondary
process do nothing.
2) Uniform treatment: given a set of quality constraints
all polluters are required to reduce BOD discharges
by the same percentage.
3) Least cost method: identical to treatment at the
polluter (see formulation on page 51 )•
4) 7-zone uniform: the polluters are divided into 7 zones
according to their location along the estuary. Within
each zone all polluters reduce BOD discharges by the
same percentage.
-------
95.
5) 3-zone geographic: the estuary is divided into 3 zones
of equal length and within each zone all polluters re-
duce BOD discharges by the same percentage.
6) 2-zone uniform: zone 1 contains all municipal polluters
and zone 2 all industrial polluters. Within each zone
all polluters reduce BOD discharges by the same percentage
7) 3-zone percentage: zone 1 consists of all polluters
currently treating over 80%. Zone 2 consists of all
polluters treating between 40 and 80% and zone 3 con-
sists of all polluters treating less than 40%. Again
within each zone all polluters reduce BOD discharges
by the same percentage.
8) Effluent charge: polluters are charged for each pound
of BOD discharged into the estuary, and in accordance
with marginal cost analysis they reduce discharges un-
til the marginal cost of reducing discharges equals the
amount of the effluent charge.
The total annual cost of each of these schemes for a quality
standard of 3 mg/1 is given in Table 21. Schaumburg gives
these costs as a present value over a discounting period of
20 years with a present value factor of 13.
TABLE ?,].
Cost Comparison of Treatment Schemes
Treatment Scheme
1) required secondary
2) uniform treatment
3) least cost
4) 7-zone uniform
5) 3-zone geographic
6) 2-zone uniform
7) 3-zone percentage
8) effluent charge
Total Cost ($ x 10°)
10.331
8.153
4.0
5.138
5.715
7.084
7.292
4.692
-------
96.
Uniform treatment and required secondary are favorite
schemes of politicians because they are "fair". From Tab-
le 21, we see that the cost of these "fair" schemes are
"unfair" to whomever pays for them.
Part of the input data for the Delaware Estuary pro-
blem is given in Section V (Tables 9 and 10). The transfer
matrix is given in Table 22 and the three distance tables
are given in Tables 23, 24 and 25. In Table 24, the infinity
sign eo indicates large distances and effectively rules out
piping across the estuary where our piping costs would not
be accurate.
We began by solving the linear case of treatment at
the polluter. Our solution gave a total annual cost of
4.115 million dollars. We believe the additional $115,000
can be attributed to further refinements made in our trans-
fer matrix. Unfortunately, Schaumburg does not report the
transfer matrix, flows, or concentrations which he used.
Schaumburg's data as well as ours was supplied by the Dela-
ware Estuary Comprehensive Study Group; and, we feel confi-
dent in comparing our results with his results. Also, it
appears that he rounded his final total costs. The details
of our solution are given in Table 26. In Table 26, the
values in column 2 are effluent BOD concentrations at the
optimal solution. It is interesting to note that of the
44 polluters on the Delaware Estuary, only 7 have 80% or
better BOD removal. From Table 26, we see that 22 (or half)
of the polluters have 80% or better BOD removal in the op-
timal treatment at the polluter solution.
-------
TABLE 22
Transfer Matrix for the Delaware Estuary
1
2
3
k
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
1
7.834B-06
1.387E-05
1.771E-05
1.806E-05
1.616E-05
1.159E-05
9.7612-06
8.1962-06
6.6332-06
5.3832-06
4.5712-06
3.384E-06
2.779E-06
2.016E-06
1.382E-06
9.352E-07
6.203E-07
4.028E-07
2.577E-07
1.835E-07
1.385E-07
9-9232-08
4.750E-08
3.218E-08
2.187E-08
1.512E-08
8.65CE-09
3.766E-09
2
9-331E-07
1.069E-05
1.706E-05
1.892E-05
1.7682-05
1.302E-05
i.noE-05
9.3892-06
7. 6432-06
6.234E-06
5.3112-06
4.650E-06
3.9512-06
3.2502-06
2.3622-06
1.6212-06
1.0982-06
7.2912-07
4.7362-07
3.0312-07
2.159E-07
1.6292-07
1.1672-07
7.9672-08
5.5892-08
3.786E-OG
2.574E-08
1.7802-08
1.0192-oS
3
5.767E-08
8.356E-07
1.240E-05
1.770E-05
1.834E-05
1.4282-05
1.248E-05
1.070E-05
8.8082-06
7.251E-06
6.217E-06
5.469E-06
4.664E-06
3.848E-06
2.8062-06
1.9312-06
1.311E-06
8.710E-07
5.663E-07
3.625E-07
2.5832-07
1.950E-07
1.397E-07
9.537E-08
6.691E-08
4.533E-08
3.081E-08
2.13LE-08
1.220E-08
5.3082-09
4
1.343E-09
2.1302-08
1.2742-05
1.7502-05
1.5282-05
1.399E-05
1.229E-05
1.0302-06
8.612E-06
7.4592-06
6.61UE-06
5.672E-06
4.701E-06
3.447E-06
2.3812-06
1.621E-06
1.079E-06
7.0232-07
4.4992-07
3.207E-07
2.4222-07
1.736E-07
1.185E-07
8.3132-08
5.6332-08
3.829E-08
2. 6482-08
1.517E-08
6.5992-09
5
3-271E-11
5.439E-10
1.157E-08
4.770E-07
1.330E-05
1.5292-05
1.5232-05
1.393E-05
1.203E-05
1.0292-05
9.0472-06
8.1082-06
7.015E-06
5.8502-06
4.323E-06
3-002E-06
2.051E-06
1.368E-06
8.023JE-07
5.722E-07
4.0812-07
3.082E-07
2.2092-07
1.508E-07
1.0592-07
7.174E-08
4.8782-08
3-3742-08
1.933E-08
8.410E-09
6
9.8652-13
1.685E-11
3.778E-10
1.7382-08
6.3252-07
1.2272-05
1.519E-05
1.510E-05
1.3782-06
1.2272-05
1.1052-05
1.0082-06
8.8332-06
7.4372-06
5.5552-06
3.8912-06
2.6722-06
1.7892-06
1.169E-06
7-509E-07
5.3592-07
4.049E-07
2.904E-07
1.9832-07
1.392E-07
9.436E-08
6.417E-08
4.439E-08
2.544E-08
1.107E-08
7
2.057E-13
3.557E-12
8.166E-11
3.921E-09
1.548E-07
3.729E-06
1.199E-05
1.4352-05
1.446E-05
1.3712-05
1.280E-05
1.196E-05
1.066E-05
9.0912-06
6.887E-06
4.873E-06
3.368E-06
2.265E-06
l.464v-o6
9-549E-07
6.820E-07
5.156E-07
3.699E-07
2.5272-07
1.7742-07
1.203E-07
8.1842-08
5.663E-08
3. 2462-08
1.413E-08
8
8.947E-1'*
1.5592-12
3.626E-H
1.782E-09
7.326E-08
1.924E-06
7.4892-06
1.221E-05
1.383E-05
1.3962-05
1.3462-05
1.2832-05
1.161E-06
9.9942-06
7-655E-06
5.458E-06
3.790E-06
2.557E-06
1.679E-06
1.082E-06
7.7312-07
5.846E-07
4.196E-07
2.8672-07
2.014E-07
1.366E-07
9.2912-08
6.4302-08
3.687E-08
1.6o6E-o8
9
3.855E-14
6.758E-13
1.589E-11
7.9532-10
3-369E-08
9-3722-07
4.044E-06
7.409E-O6
1.1772-05
1.356E-05
1.386E-05
1.369E-06
1.267E-05
1.1082-05
8.6302-06
6.224E-06
4.3532-06
2.950E-06
1.943E-06
1.254E-06
8.9702-07
6.7882-07
4.873E-07
3.3322-07
2.34LE-07
1.5882-07
1.080E-07
7.479E-08
4.2902-08
1. 869E-08
10
1.68TE-1&
2.973E-13
7.0572-12
3.587E-10
1.5572-08
4.5242-07
2.089E-06
4.0812-06
7. 4182-06
1.191E-05
1.356E-06
1.416E-05
1.356E-05
1.2132-05
9.672E-06
7.0832-06
5.000E-06
3.4082-06
2.253E-06
1.4582-06
1.044E-06
7.905E-07
5.6792-07
3.885E-07
2.730E-07
1.8522-07
1.261E-07
8.732E-08
5.0112-08
2.1842-08
-------
transfer Matrix (Continued)
1
2
3
k
6
7
8
9
10
n
12
13
15
16
17
18
19
20
21
22
23
2k
25
26
27
28
29
30
n
8.916E-15
1.57CE-13
3-771E-12
1.938E-10
8.555E-09
12
1.229E-06
2.1f86E-06
8.63^-06
1.213E-06
1.382E-05
1.38915-05
1.2703-05
l.o'tyE-05
7.82IE-06
5.580E-06
3.629E-06
2.5lf2E-06
1.1G2E-06
8.960E-07
6.MtOE-07
l».l»OuE-07
3.099E-07
2.103E-07
l.lf32E-07
9.9232-08
5.697E-08
2.069E-12
1.0732-10
lf.800E-09
7.255E-07
1.502E-06
3.0233-06
5.7672-06
8.872E-06
1.235E-05
1.3562-05
1.3UE-05
1.12liE-05
8.6032-06
6.231E-06
4.315E-06
2.88UE-06
1.8772-06
1.31*72-06
1.022E-00
7.352E-07
5.036E-07
1.639E-OY
1.136E-07
6.526E-03
13
2.5982-15
1.118E-12
8.220E-08
8.7673-07
1.8li*E-o6
3.600E-06
5.83}tE-o6
8.909F.-06
1.20l(E-05
1.27UE-05
1.175E-05
9.363E-06
6.928E-06
J».859E-o6
3.271E-06
2.ll*lE-06
1.5'*OE-o6
1.170E-06
5.776E-07
l».o66E-07
2.7632-07
l.G8iiE-07
1.306E-07
7.512E-08
3.279E-08
6.161E-13
3.2'llE-ll
l.l»80E-09
4.666E-08
5.12GE-07
1.082E-06
2.205E-06
3.690K-06
5.930E-06
8.782K-06
1.130E-05
1.180P.-05
9.965E-06
7-593R-06
3.681E-06
1.331E-06
9.597E-07
6.588E-07
2.15UE-07
1.'19515-07
8.606E-08
3.760E-08
15
6.122E-16
2.67J*E-13
l.UlGE-ll
6.550E-10
2.101E-08
l.lO'lE-07
2.39QE-07
5.152E-07
1.080E-06
1.865E-06
3.1^12-06
5.006E-06
1.096E-05
1.0lf7E-05
8.1(21E-06
6.1G2E-06
2.061E-06
1.572E-06
1.13613-06
7.8ll)E-07
3.75^-07
2.5652-07
1.782E-07
1.028E-07
16
2.085E-16
2.280E-10
3.950E-08
8.6l»OE-o8
1.890E-07
1».032S-07
7.105E-07
1.229E-06
2.03'/E-06
3.139E-06
5.328E-06
9.817E-06
7.500E-06
5.1*132-06
3.691B-06
2.703E-06
2.075E-06
1.507E-06
1.0l*lE-06
7.3653-07
5.033E-07
3.JA9E-07
2.1»03E-07
1.391E-07
6.106E-08
17
7.202E-17
1.295E-15
3,l82E-llf
1.705E-12
7.992E-11
2.617E-09
l.l*09E-08
3.101E-08
6.8l»2E-08
2.630E-07
U.620E-07
7.815E-07
1.238E-06
2.2lliE-06
9.108E-06
8.6i}3E-o6
3.56GR-06
2.76612-06
2.0252-06
l.JfOSR-06
6.871E-07
3.307E-07
1.92CE-07
18
2.539E-17
1;.573E-16
2.845E-U
9-368E-10
5.073E-09
1.121E-08
2.1i8GE-08
5.U02E-OS
1.719E-07
4.733B-07
8.700E-07
2.063E-06
8.50-'*E-o6
7.G97E-06
6.051E-06
U.630E-06
2.705E-06
1.900E-06
1.361E-06
6.517E-07
k.^QkE 07
2.693E-07
1.195E-07
19
9.3382-18
1.683E-16
'f.lI*7E-15
2.233E-13
1.053E-n
1.892E-09
4.193E-09
9-337E-09
2.036E-08
3.670E-08
1.L2CS-07
1.832E-07
8.1*592-07
7.892E-06
7.l'f3S-o6
5.760E-06
3.536E-06
2.528E-06
1.6331:-06
1.26QE-06
8.959E-07
6.359S-07
3-7C3E-07
20
3.625E-18
6.537E-17
1.612E-15
8.693E-11*
1.360E-10
7.1413E-10
3.6732-09
7.
1.
3-
1.U52E-08
2.597E-08
4.500E-08
.3!*9E-o8
.3G6E-07
.50GE-07
8.639E-07
2.088E-06
l».3'»OE-o6
7.221E-06
6.569E-06
3.27GE-06
2.1*25E-06
1.725E-06
1.227E-06
8.83l»E-07
5.360E-07
2.U37E-07
00
-------
Transfer Matrix (Continued)
21
22
23
25
28
1
2
3
k
5
6
7
8
9
10
li
12
13
Ut
15
16
17
18
19
20
21
22
23
2k
25
26
27
28
29
30
1.829E-18
3.298E-17
8.137E-16
k.3}OE-lk
2.075E-12
6.883J3-11
3-755E-3.0
8.3'(013-10
1.66l;33-09
'(.07913-09
7.301(13-09
1.3^313-08
2.29713-08
3.76013-08
7.11&3-08
1.C20K-07
1(.660J3-07
1.132E-06
2. '.'7713-06
';.6l;5J3-0($
6.07^-06
5.90C3-06
5.0MI3-06
3. 89923-06
2.979--3-06
2.1YYJ3-06
1.586E-06
l.i.6':E-06
7.25';l3007
3.36CG-07
l.K&S-lS
1.991E-17
U.913E-16
2.65333-!^
1.25'iE-12
lt.l£ni?-Ti.
2.272E-10
5.0'|Q3-1D
1.359E-C9
2.1)73E-O~>
'». '(7913-02
8.032K-0?
1.39G3-C3
2.289E-0;
l».3'i'iE-Cc
1.117E-07
2.888E-07
7.123E-07
1.602E-C-S
3.178E-C-S
'i.699l3-Co
5.5l'(2-C-5
5.119K-C5
lf.l96E-c'
3.306E-C5
2.VC5E-C-S
1.81(132-C5
1.373E-C:
8.735E-"
^.1033-07
6.569S-19
1.1895-17
2.93^E-i6
1.58'iE-l4
7.'(95E-13
2.'(87E-11
1.359E-10
3.021E-10
6.759E-10
l.':8lE-09
2.68513-09
lt.83GE-09
8.302K-09
1.376E-08
2.616E-08
.6.7562-08
1.759K-07
1(.390E-07
1.007E-06
2.07'fl3-0o
3.290I3-06
li.PO': 13-06
H.7'i7E-06
k, 263D-06
3.5^32-06
2.75013-06
2.1052-06
1.605E-06
1.0'iGE-06
5.009E-07
3.87l(E-l9
6.989E-18
1.725E-16
9.316E-15
1».1(09E-13
lJ(6'lE-ll
8.002E-11
1.78013-10
3.982E-10
8.73213-10
1.563E-09
2.8'(3f';-09
k. 95013-09
8.13313-09
l.S^E-OS
'(.015E-08
1.051E-07
2.6ii6E-07
6.1592-07
1.30213-06
2.162E-06
2'.90U3-Oo
3.587K-06
3.99813-06
3.6o6E-o6
2.9603-06
2.363E-06
1.856E-06
1.255E-06
6.128E-07
2.512E-19
l».586E-l8
1.132E-16
6.1l'tE-l5
2.69'iE-l3
9.613E-12
5.256E-11
1.169E-10
2.61713-10
5.73913-10
1.0'(3J3-09
1.67033-09
3.257E-09
5.35^3-09
1.023E-08
2.652R-08
6.966E-08
1.76313-07
1|-.1'>1E-07
8.690E-07
1.515E-06
2.o81iE-06
2.6692-06
3.2J(7E-o6
3.3o7E-o6
3.011iE-06
2.533E-06
2.05913-06
lJi'(5E-o6
7.21213-07
1.6kkE-l9
2.966E-18
7.323E-17
3.955E-15
1.873E-13
6.220E-12
3.1(02E-11
7.567E-11
1.69'iE-lO
3.717E-10
6.7'(3E-io
1.21P.E-09
2.3HK-09
3.1(72E-09
6.623E-09
1.724E-08
!(.5UlE-o8
1.15'(E-07
2.730E-07
5.93.^-07
1.02913-06
l.'(l;2E-06
1.91313-06
2.'(25E-06
2.722E-06
2.666E-06
2.6362-06
2. 259F-06
1.672E-06
8.587E-07
1.CA8E-19
1.892E-18
4.671E-17
2.523E-15
1.195E-13
3.969E-12
2.171E-11
t.830E-n
1.082E-10
2.37^-10
1*.307E-10
7.7'(1E-10
l.3i(9E-09
2.22013-09
1..237E-09
1.10'iE-08
2.916E-OS
7.438E-08
1.76913-07
3.880E-07
6. 833E-07
9.7012-07
1.31^-06
1.720E-06
2.020E-06
2.31GE-06
2.563E-06
2.399E-06
1.917E-06
1.023E-06
6.73^-20
1.215E-18
3.001E-17
1.621E-15
7.676E-li
2.553£-l2
1.395E-31
3.105E-11
6.953E-11
1.526E-10
2.769E-10
1+.979E-10
8.6COP.-10
1.1(29E-09
2.728F.-09
7.320E-09
1.883E-OG
*(.8l7E-o8
1.15.U3-07
2.5U3E-07
k. 528E -07
6.'i9liE-07
8.92'(E-07
1.195E-06
l.i)ii5E-o6
1.7'i'iE-o6
2.106E-06
2.35213-06
2.123E-06
1.193E-06
3.316E-20
5.983E-19
1.U78E-17
7.983E-16
3.783E-1H
1.257E-12
6.875E-12
1.530E-11
3.'»27E-11
7.5235-11
1.336E-10
2.1(5613-10
1(.283E-10
7.05?J3-10
1.31)813-09
3-523E-09
9.3l(2E-09
2.399E-08
5.767H-08
1.287E-07
2.326E-07
3.379^-07
4.731E-07
6.5132-07
8.1'(3E-07
1.036E-06
1.356E-06
1.720E-06
2.21(713-06
I.iil6i:-o6
1.177E-20
2.124E-19
5.2U6E-18
2.835E-16
1.3U3E-1U
1».1*63E-13
2.M2E-12
5.'i36E-l2
1.218E-11
2.67^-11
U.SjUE-ll
8.73333-11
1.523S-10
2.509E-10
k. T9713 -10
1.256E-09
3.336E-09
8.593S-09
2.075E-06
U.665E-08
8.52i(E-o8
1.2502-07
1.772E-07
2.1J85E-07
3.173B-07
4.165E-07
5.6912-07
7.655E-07
1.129E-06
l.i(15E-o6
-------
100,
TABLE 23
Mileages from Polluters to Sections on the Delaware Estuary
Polluter 1
Nuiiiber
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
31
32
33
34
35
36
37
30
39
40
41
42
43
44
0
4
4
5
5
8
28
28
23
35
35
37
37
37
37
40
40
40
40
44
44
48
48
43
48
48
4G
48
52
52
56
53
53
56
56
5G
63
63
63
65
65'
67
73
70
2
4
0
0
4
4
6
25
25
25
31
31
33
33
33
33
36
36
33
36
40
40
44
44
44
44
44
44
44
48
4C
52
52
52
52
52
52
59
59
59
61
61
63
69
74
3
G
4
4
0
0
3
21
21
21
27
27
29
29
29
29
32
32
32
32
36
36
40
40
40
40
40
40
40
44
44
48
48
48
48
48
48
55
55
55
57
57
59
65
70
4
12
8
D
4
4
0
17
17
17
23
23
25
25
25
25
28
28
28
23
32
32
36
36
36
36
36
36
36
40
40
44
44
44
44
44
44
51
51
51
53
53
55
61
66
5
16
12
12
C
8
4
13
13
13
19
19
21
21
21
21
24
24
24
24
23
23
32
32
32
32
32
32
32
36
36
40
40
40
40
40
40
47
47
•
-------
101,
TABLE 23 CON'T.
Mileages from Polluters to Sections on the Delaware Estuary
Polluter 13
Kupber
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
10
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
3G
37
30
39
40
41
42
43
44
42
31
3G
33
33
20
15
15
15
9
9
7
7
7
7
4
4
4
4
0
0
0
4
4
4
4
4
4
8
8
12
12
12
12
12
12
19
19
19
21
2-1
23
29
33
17
46
38
40
37
37
32
19
19
19
13
13
11
11
11
11
7
7
8
8
4
4
2
0
0
0
0
0
0
4
4
8
8
8
8
8
8
15
15
15
17
17
19
25
29
1C
50
42
44
41
41
34
23
23
23
17
17
15
15
15
15
12
12
12
12
8
*-»
o
4
4
4
4
4
4
4
0
0
0
0
4
4
4
4
11
11
11
13
13
15
21
25
19
54
46
47
45
45
38
27
27
27
21
21
19
19
19
19
1C
18
16
1C
12
12
8
8
8
8
8
8
8
4
4
2
2
0
0
0
0
7
7
7
9
9
11
17
21
20
57
50
50
49
49
42
30
30
31
23
23
23
21
21
21
20
20
18
IS
1G
14
10
10
10
12
12
12
12
8
8
4
4
4
4
4
A
3
3
3
5
5
7
13
17
21
GO
52
53
52
52
44
31
31
34
25
25
25
23
23
23
23
23
SO
20
1C
17
13
13
13
14
14
14
14
10
10
6
C
6
6
6
e
0
0
0
?,
2
4
10
14
Estuary
22 23
62
53
55
53
53
45
33
33
36
27
27
20
25
25
25
25
25
22
22
20
19
15
15
15
IS
1G
16
16
12
12
8
8
8
8
8
8
2
2
2
0
0
2
8
12
64
55
57
55
55
47
34
34
3G
29
29
30
27
27
27
27
27
24
24
22
21
17
17
17
IS
18
18
18
14
14
10
10
10
10
10
10
4
4
4
2
2
0
G
10
Section
24 25
6G
57
59
57
57
4S
36
3G
40
31
31
32
29
29
29
29
29
2G
2G
24
23
19
19
19
20
20
20
20
16
16
12
12
12
12
12
12
6
6
6
4
4
2
4
8
GO
59
61
59
59
51
3C
30
42
33
33
34
31
31
31
31
31
28
28
25
25
21
21
21
22
22
22
22
IS
18
14
14
14
14
14
14
8
8
8
6
G
4
2
6
Numbers
26 27
68
59
C2
GO
60
51
39
39
44
34
34
36
32
32
32
33
33
29
29
28
25
21
21
21
24
24
24
24
20
20
16
16
16
16
16
16
9
9
9
7
8
G
0
4
63
51)
G4
62
62
51
39
39
46
34
34
38
32
32
32
35
35
29
23
30
25
21
21
21
26
23
23
23
22
22
18
18
18
16
18
IS
11
10
10
8
10
6
2
2
2G
69
60
GG
G4
54
52
39
39
40
35
35
40
33
33
33
37
37
30
30
32
26
22
22
22
23
23
28
28
24
24
20
20
20
1G
20
20
13
11
11
9
12
8
4
0
29
71
G2
G9
67
67
54
39
39
51
30
38
43
33
36
36
40
40
33
33
35
29
25
25
25
31
31
31
31
27
27
23
23
23
20
23
23
16
15
15
13
15
10
7
3
30
74
65
72
71
71
57
42
42
55
41
41
47
39
39
39
44
44
36
30
39
32
28
28
28
35
35
35
35
31
31
27
27
27
23
27
27
20
18
18
16
19
14
11
7
-------
102.
TABLE 24
Mileages from Polluters to Regional Plants on the Delaware Estuary
Polluter 101
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
9
5
CO
CO
CO
4
21
21
CD
27
27
CO
29
30
31
09
CD
32
32
CO
36
38
40
41
03
CO
CO
00
00
CD
CD
CO
CO
48
CO
CD
CO
55
56
58
CO
61
CD
69
102
CO
CO
6
4
2
03
00
CO
18
oo
00
26
CO
00
00
33
32
00
25
00
oo
CO
CQ
33
32
31
30
33
34
35
36
37
CO
38
39
43
CO
CO
CO
47
oo
55
CO
Regional Plant Number
103 104 105 106
18
14
oo
CO
oo
4
11
11
CO
17
17
CO
19
20
21
CO
09
22
22
oo
26
28
30
31
CO
CO
CO
00
oo
CO
CO
OB
CD
38
CO
CD
CD
45
46
48
09
51
00
59
35
31
00
CD
CO
23
7
7
00
2
2
00
1
2
3
00
oo
3
5
00
7
10
11
12
00
00
00
CO
00
oo
CO
00
CO
19
CO
00
00
26
28
30
00
33
00
24
42
38
CO
CO
CO
28
15
15
00
9
9
CO
7
6
5
CO
00
4
3
CO
3
3
4
5
cc
oo
CO
CO
CO
CO
CD
00
00
10
00
CO
oo
17
19
21
oo
24
00
23
48
44
CO
CO
00
34
21
21
00
15
15
03
13
12
11
CO
CO
10
8
CO
4
3
2
1
oo
oo
oo
oo
oo
CO
CO
00
oo
5
GO
00
CO
11
13
15
CO
19
oo
19
107
CO
00
38
36
34
00
CO
CO
18
co
CO
10
CO
00
CO
6
7
CO
CO
3
cr
CO
CO
CO
2
1
1
2
5
6
7
8
9
00
9
10
14
CO
00
oo
18
oo
26
00
108
00
CO
45
43
41
00
00
00
25
CO
00
17
CO
CO
00
13
14
00
CO
10
00
CO
CO
CO
8
7
6
5
2
1
1
1
2
CD
3
4
8
00
CO
oo
12
oo
20
CO
109
63
59
00
oo
CO
49
34
34
CO
29
29
CO
27
26
25
OD
00
24
21
CO
18
18
17
16
OD
CO
CD
CO
CD
CO
CO
CO
00
11
CO
00
CO
4
2
1
CO
2
CO
7
-------
103.
TABLE 2 5
Mileages from Regional Plants to Sections on the Delaware Estuary
Section 101
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
8
4
0
4
8
12
15
17
18
20
22
23
25
27
29
32
34
38
42
46
50
52
53
55
57
59
59
59
60
62
102
9
8
4
0
4
7
9
11
13
14
16
18
20
22
25
29
33
37
41
44
48
49
51
53
55
56
58
60
63
67
Regional Plant -Number
103 104 105 106
12
10
7
4
0
3
5
6
8
10
11
13
15
17
20
22
28
30
34
38
40
41 "
43
45
47
47
47
48
50
53
37
33
29
25
21
17
14
12
10
8
6
4
2
0
3
7
11
15
19
21
23
25
27
29
31
32
32
33
36
39
44
40
36
32
28
24
21
19
17
15
13
11
9
7
4
0
4
8
12
14
17
19
21
23
25
25
25
26
29
32
48
44
40
36
32
28
25
23
21
19
17
15
13
11
8
4
0
4
8
10
13
15
17
19
21
21
21
22
25
28
107
48
44
40
36
32
28
25
23
21
19
17
15
13
11
8
4
0
4
8
12
14
16
18
20
22
24
26
28
31
35
108
52
48
44
40
36
32
29
27
25
23
21
19
17
15
12
8
4
0
4
8
10
12
14
16
18
20
22
24
27
31
109
67
63
59
55
51
47
44
42
40
38
36
34
32
30
27
23
19
15
11
7
4
2
0
2
4
6
6
8
10
14
-------
TABLE 26
Solution for Treatment at the Polluter for the Delaware Estuary
104,
Polluter
Numbe r
1
2
3
4
5
6
*7
8
*9
*10-
*11
12
*13
14
15
16
17
18
*19
* 20
*21
22
BOD
rag/1
18.6
33.4
1.3
0.8
82
80
104
67
62
45
50
190
126
40.3
43.1
93.4
30.5
69.9
31
40
92
79.9
T
% removal
85
85
85
0
85
55
85
75
85
90
85
35
85
0
70
65
35
35
85
85
80
35
Polluter
Number
*23
*24
25
26
*27
*
28
29
*30
*31
32
*33
*34
*35
36
37
38
39
40
41
42
43
44
BOD
mg/L
648
10.5
32.7
32.3
93
43
62.7
41
2.5
70.1
36
1233
23
10.4
172
171.5
212
128
80.2
187.1
1.1
159.9
% removal
93
87
85
85
80
80
60
80
57
35
85
72
82
0
20
40
0
35
80
40
35
40
* polluters that increased levels of treatment
-------
105
This solution, obtained originally by Thomann and
Marks, was until recently the best known solution. How-
ever, as we are about to show, the general mixed case
gives further cost reductions which are quite significant.
Of the 44 polluters located along the Delaware Estuary,
22 treat domestic wastes and 22 treat industrial wastes
(see Table 10). We have restricted ourselves by not allow-
ing any mixing of the domestic and industrial wastes.
Thus, all the polluters which pipe to regional plants will
be domestic waste polluters. This exclusion of the indus-
trial polluters from the regional plants adversely effects
the optimal solution that we can reach. However, at the
present time we do not have any cost data for treating joint
industrial-municipal wastes. Fortunately, the Delaware River
Basin Commission has recently been awarded a Research and
Development Grant to construct a pilot plant to investigate
a chemical-biological treatment process for joint industrial-
municipal wastes, capable of attaining at least 88 to 93 per-
cent removal of major pollutants. Design, operating and cost
information is to be obtained for a 80 MGD regional treat-
ment complex. When this cost data becomes available, indus-
trial wastes could also be piped to the regional plants with-
out any complications in our model and would undoubtedly im-
prove the final optimal solution. The total annual cost of the
-------
106
optimal solution (for the 9 fixed regional plant locations)
is 2.291 million dollars and is the best known solution
at this time. The details of the solutions are given in
Figure 18 and Table 27.
In Figure 18, we see that polluters 23 and 34, both in-
dustrial waste polluters, mix in a common by-pass pipe. This
unexpected introduction of by-pass piping into the optimal solu-
tion is easily explained. Both polluters 23 and 34 have low
flows with tremendously high BOD concentrations. Since piping
costs are based strictly on flow we get a lot for our money by
piping a low flow (high amount of pollution) a large distance.
Since most of the pipelines proposed in Figure 18 are not in
heavily populated areas, we feel that our piping costs are rea-
sonably accurate.
-------
FIGURE 18
Optimal Solution
107
-------
TABLE 27
Data for Optimal Solution
108,
regional plant % removal
number
1
2
3
4 70
5
6
7 70
8
9 57.9
polluter
number
7
9
12
13
24
27
31
34
Bon
mg/L
105
62
155
126
10.5
93
2.5
1777
% removal
85
85
47
85
87
80
57
60
Note: The remaining 36 polluters don't increase treatment.
Regional Plant 4;
MGD
(22.4 + 7.0)
(2.0 + 0.7)
32.1
Regional Plant 7:
MGD
118.0
(6.6 + 8.6 + 1.4)
134.6
Regional Plant 9;
MGD
(0.902 + 0.467)
1.369
from polluters 10 and 11
from polluters 19 and 21
discharge to section 14
from polluter 20
from polluters 28, 30 and 33
discharge to section 17
from polluters 38 and 39
discharge to section 23
Note: polluter 38 discharges 0.097 MGD to section 21 and
polluter 39 discharges 4.332 MGD to section 21.
By-pass piping:
MGD
(0.3 + 1.1)
0.46514
0.93486
from polluters 23 and 34
discharge to section 23
discharge to section 27
-------
109,
To complete this section, we intend to make a longhand cal-
culation of the optimal solution to verify that it is correct.
Costs at 8 polluters that treat:
Polluter
Number
7
9
12
13
24
27
31
34
$ x 106
.026359
.127177
.477324
.013179
.413668
.013085
.024661
.249282
[170,110 - (107.6-155-8.33)]
199
• [32650 - (1.1-1777-8.33)]
1.344735
Costs at 3 treatment plants:
Plant 4:
r4 = 0.7
V4
i f22-4'322 + 7-217 + 2-124 + 0.7-32ll _
~ [_22.4-456 + 7-333 + 2-207 + 0.7-458J
v4 = .7 + .30689 - (.7 x .30689) = .792067
cost = .39376 (32.1)0'75 (.2920673 + .193113) = .170539
Plant 7:
r? = 0.7
v = 1 PH8-159 + 6.6-141 + 8.6-142 + 1.4-156"] _
7 [_118-265 + 6.6-217 + 8.6-203 + 1.4-240J
v? = .7 + .39247 - (.7 x .39247) = .81774
cost = .39376(134.6)°*75 (.317743 + .107533) = .518496
-------
110,
Plant 9:
rA = .579
- - T F. 902.171 + .467>212~| = .
V9 L.902-286 + .467-212J
vg = .579 + .290579 - (.579 x .290579) = .701333
cost = .39376 (1.369)0*75 (.2013333 + .2094213) = .008644
Consolidation of pipes:
.598
Plant 4: 1865 x 2 x 29.4
,.598
= .028169 "\
1865 x 2 x 0.7
1865 x 5 x 2.7
Plant 7: 1865 x 3 x 118
1865 x 3 x 1.4
,.598
.598 =
.598
.598
1865 x 4 x 10'
1865 x 2 x 16.6
,598
.003013
.016889 y
.097003
.006842
.029562
.020013
Plant 9: 1865 x 2 x .902'598 =
1865 x 2 x 1.369'598=
.003506
.004500
i
,048071
153420
,008006
By-pass piping:
1865 x 7 x .3
.598
1865 x 10 x 1.4
.598
1865 x 6 x .93486
.598
.006354
.022806
.010748
,039908
TOTAL = 2.291819
-------
Ill
Referring back toSchaumburg's [12] cost comparisons in
Table 21, we see that this total of 2.291819 million annual
cost is only about half the previously best known solution
of 4 million annual cost for pure treatment at the polluter.
It is only about a fourth of the 8.153 million annual cost
of "uniform treatment" the currently favored policy. Even
allowing for some crudeness in our cost data and transfer
coefficients, it would seem evident from this result that
Interstate Agencies administering regional treatment com-
plexes could substantially reduce the enormous costs invol-
ved in maintaining and improving water quality.
The present mathematical model also provides an excell-
ent planning tool for evaluating the relative costs of alter-
nate water quality goals. It should help establish more
realistic goals and development plans. In future work, we
shall concentrate on incorporating projected future demand
and larger numbers of alternate plant sites with a view to
establishing long term time dependent development plans. We
shall undoubtedly also have to incorporate into our proce-
dure a periodic updating of the transfer coefficients to
account for any substantial alterations in estuary flow caused
by the use of by-pass piping or regional treatment plants.
-------
VIII. ECONOMIC CONSIDERATIONS
"Preliminary Ideas"
The economic implications of regional water pollution con-
trol systems lead to some very difficult questions. For example,
who will pay for the regional treatment plants? Do new polluters
wishing to locate along the Delaware have the same right to pol-
lute as polluters already established? Can this right to pollute
be taken away or can a polluter be penalized by being made to
pay for some part of a regional treatment system while other pol-
luters pay nothing? We do not intend to resolve any of these
issues at this time. However., we sketch out here an approach
to implementing an efficient regional treatment system to show
the relevance of the results of section VTL The problem is
viewed here in terms of collective or joint decision making
between the public which desires better water quality and ex-
presses this via governmental agencies and the polluters or
users of the water resource. The latter may be roughly broken
down into industrial and municipal users. The issues involved
in determining the "true" goals of various groups competing
for the resource are complex since these groups may overlap.
Decision making in a democratic society is often a diffi-
cult process. It is well-known that there are problems in-
volved in finding a satisfactory voting mechanism to arrive at
agreement and the determination of a concensus is even more
complex. Suppose governmental decision makers must decide on
a suitable system of treatment and charges based on water quali-
ty standards. One possible plan would be to ask each group
involved for their proposal on what water quality should be and
-------
113.
how to achieve it. Each proposal could then be evaluated to
determine its costs (both social and direct). This plan could
become a highly expensive exercise and would perhaps not yield
satisfactory proposals since groups often have very hazy views
on the social implications of their positions. Thus the effec-
tive use of sophisticated mathematical techniques and models
which are reasonable approximations to reality may be a critical
ingredient to efficient decision making in a democratic society.
Without this available, groups may be making choices in al-
most complete ignorance of their implication or the decision
makers may be forced into a quasi-dictatorial situation where
only one or two alternatives can bs examined. By building a
model representing the various decision making units, indus-
trial and municipal waste dischargers and using goals specific
to various regions of the estuary, we can know the rough im-
plications of any proposal in a micro-sense in contrast to
the cost benefit analysis where some mystic aggregated figures
are conjured up to count the total costs and total benefits.
Even ignoring the somewhat crude measures used in cost benefit
studies, the usefulness of an aggregation of that type as a
basis for decision making has been highly overrated. If de-
cisions are to be made on the basis of a reconciliation of the
positions of various groups, each one must be able to judge the
benefits and effects of any proposal on its own situation.
The most popular scheme for the elimination of water pollu-
tion is the effluent charge. The argument usually made is that
water quality is a good that is being used up by polluters and
should be priced as any other good. The price is called an
-------
114.
effluent charge. Setting water quality standards in effect
fixes the amount of water quality for sale. One proposed way
of determining the correct effluent charge is by an iterative
process. Some measure of water quality utilization is first
defined as a basis for pricing. Prices are set at an arbitrary
level for the last section of the waterway. For a given pricing
scheme if the amount of water quality used in some section is
greater than that available, the price is increased; if more
is available than required, the price is decreased. Assuming
convergence of this process we end with a situation that not
only meets the goals set forth but does so in a manner that
minimizes the total waste treatment costs for industry in
that region. The approach has a certain appeal for the follow-
ing reasons.
1. As a starting procedure cost estimates of pollution
control can be made to determine the approximate
range for correct prices. Furthermore by using mathe-
matical programming, dual prices are obtained which
can be used to judge how far from optimal the esti-
mated prices are.
2. The effluent charge system is efficient in the sense
that other approaches such as uniform treatment will
cause total waste treatment costs to be higher.
3. The prices give correct incentives to the polluters.
A polluter may reduce his costs by discovering alter-
nate products or processes.
What are some of the objections to an effluent charge approach?
1. The concept of auctioning off the water quality re-
source has some defects. In a static model of firm
-------
115
behavior an equilibrium price is determined and then
a supply of the good is produced in response. In
reality, behavior is dynamic. There is the possibi-
lity of sudden surges or peaks of demand perhaps re-
sulting from the variability of production processes.
Here there is no limitation on the amount of the
water resource used by a firm; it may be purchased
until total exhaustion of the water quality resource
is achieved. The problem is we are trying to super-
impose a technique that has validity in a static
situation on a highly dynamic phenomenon.
2. It may be to the interests of the polluters to mis-
lead the regional authority as to their actual treat-
ment costs. By indicating that treatment costs are
very high a low effluent charge would be set. There
may be a considerable institutional lag before the
charge can be reset to a higher value.
3. The validity of the static pricing allocation model
is doubtful. The requirement for convexity is proba-
bly not met.
4. There is a substantial equity issue. Assuming a
charge scheme were to work correctly, firms might in-
cur substantial cost differences. For this reason
political objections to the effluent proposal may be
considerable.
With an ordinary good there can be no justification for re-
fusing to pay a price. The case of water pollution is some-
what different since originally water quality was not priced.
-------
116,
The good was available and served as an inducement for indus-
tries to locate on the waterway. Implicitly a right was given
to use whatever water quality was available. With effluent
charges there is no cognizance of these prior rights. A new
factory would be subject to the same type of effluent charge
as earlier polluters.
As indicated earlier we shall approach the problem of a
viable regional treatment system in terms of a collective de-
cision making problem. The solutions to the mathematical mo-
del will provide guidelines and implication of various situa-
tions. We shall use a combination of the following:
1. Effluent standards for each polluter.
2. Marginal effluent charges (payments) to induce appro-
priate reductions in discharge.
3. Taxation on polluters based on the benefits derived
from the pollution services of the regional authority.
In the present study we consider the benefits of re-
gional plants and piping.
4. The use of licenses that are saleable in order to
reallocate effluent standards and rights to utilize
the regional treatment facilities.
For the proposed implementation of the regional treatment
system on the Delaware we indicate how a specific financing
scheme could be developed.
-------
117
(To Be Filled In)
Cost Data For
Financing Regional Treatment System
Polluters 112 3
2
3
1 Costs to each polluter without regional treatment
in order to meet the standard (i.e., costs to pollu-
ters using treatment at polluter only obtained from
optimization - should total 4 million annual cost)
2 Costs to each polluter with regional treatment
3 Net Savings - difference of 1 and 2
The total of column 3 should be positive reflecting the savings
obtainable with regional treatment. Let a be the factor which re-
duces the total of 3 to the cost of the regional treatment system.
This factor is the proportion of net savings that each pollu-
ter must contribute to the cost of the regional treatment sys-
tem. Note that a polluter may contribute to regional treat-
ment even though he does not actually use the regional treat-
ment plant since under the effluent standard he benefits by
being able to discharge a greater amount of effluent than with-
out regional treatment.
In return for payment of the benefit tax polluters are
allocated allowable amounts of effluent discharge. These
amounts are listed below:
(To Be Filled In)
Polluters Discharge Amounts
-------
118,
Some.of the polluters, because of high costs of treatment at
their own plant, use the regional plant. The benefit tax covers
their use of this facility.
The right to discharge wastes and to use the regional
treatment plant will be considered as saleable goods. The
regional authority sets a limit purchase price based on the
value of the dual variables.
(To Be Filled In)
Dual Variable Dual Variable
Polluter Effluent Discharge Regional Plant
The granting of transferable property rights may be an
efficient method of coordinating different standards and use
of the regional treatment system. For a given quality goal the
model determines effluent discharge standards. Each of the
present polluters may be issued a license that conveys a right
to discharge an amount of waste not to exceed the effluent
standard. This right may be given for a fixed number of years.
However during that time the polluter would have the right to
sell off part of all of his rights to pollute. For example, a
new plant finding the Delaware an exceptionally desirable loca-
tion would have the possibility of purchasing pounds of BOD
from existing polluters. The possibility of this sale would
-------
119,
induce present dischargers to find process changes or treatment
economies which would allow for a profitable sale. Over time
the use of the stream for waste disposal is allocated to the
users deriving the most benefit. A useful byproduct of having
a market in waste discharge rights would be the possibility of
valuing the stream among different uses. License fees for fish-
ing and other recreational uses may be measured against the to-
tal market value of pollution rights and the relative movements
of these numbers may at least give some aid in resetting the
quality goals.
License fees may also be used to finance regional treatment
plants. According to the solution of the model polluters are
allocated licenses to use the treatment plant based on the pollu-
ter's cost of waste treatment. After the inital allocation has
been made., any other plant wishing to make use of the regional
treatment facility would have to purchase a license from a cur-
rent user for the BOD amount and pay for piping. Only in the
case where the polluter purchasing the license would have greater
savings than a current user would there be a sale. Thus with
license fees we obtain from among the groups of potential users
the ones with the greatest savings from use of the regional
facility. Note that all users may pay for piping connections
to the plant.
-------
APPENDIX A
(Nomenc la ture )
Section II and the notation in that section are inde-
pendent of the rest of the report. The meaning of symbols
used in Section II apply only to that section while the sym-
bols defined below are used throughout the remainder of the
report.
a. . the change in DO in section i of the estuary due to
*•* a change in input of BOD in section j of the estuary,,
mg/1/lb/day where a. .eA
B. the maximum amount of BOD that can be removed for a
1 particular s. at polluter i, Ib/day
c. the required change in DO in section i, mg/1 where
ciec
d. . the distance from polluter j to the center of section
i, miles
the dist
plant i_, miles
2
d. . the distance from polluter j to regional treatment
d. . the distance from regional treatment plant j to sec-
1^J tion i, miles
]^
D . . the piping cost coefficient associated with a parti-
-1 cular d* . , $/MGD-yr
f . . the flow from polluter j to section i, MGD where
13
G the set of material balance at a polluter constraints
o
G the set of material balance at a regional plant con-
straints
3
G the set of quality constraints
j . the BOD concentration of effluent from polluter i,
1 Ib/MG where j . e J
j . the untreated BOD concentration of influent to pollu-
1 ter i, rag/1
-------
j J2, J3, and J4 the BOD concentrations used in the
piecewise linear cost function, lb/MG-10~3
K a matrix with all elements zero except for a unit
m'n element in the m*-*1 row and the nfch column
K the piecewise linear cost function, $xlO
K the cost of BOD removal at a waste water treatment
c plant $xlQ3/MGD
kl, k2, and k3 the cost coefficients used in the piece-
wise linear cost function $/lb/MG-10~"3
p.. the flow from polluter j to treatment plant i, MGD
1-) where P^6?
Q.. a flow corresponding to a particular f.., p.., or
1-} t- • (not necessarily equal) , MGD ^ 13
Q the design capacity of a waste water treatment plant,
MGD
q. the fractional removal of BOD at polluter i, dimen-
1 sionless
r. the fractional removal of BOD at regional treatment
1 plant i, dimensionless where r.eR
s. the cost of additional removal of BOD at polluter i,
1 $/lb/day
t.. the flow from treatment plant j to section i, MGD
13 where t..eT
U. the BOD discharge from polluter i, Ib/day
v- the translated fractional removal of BOD at regional
1 treatment plant i dimensionless
w the fractional removal of BOD at a waste water treat-
ment plant, dimensionless
UB. upper bound concentration of polluter i
LB. lower bound concentration of polluter i
-------
APPENDIX B
(References)
[1] ABT Associates Inc., Incentives to Industry for Water
Pollution Control: Policy Considerations, FWPCA con-
tract no. 14-12-138, Dec. 1967.
[2] Wright, J.F., "River Basin Management, An Interstate
Approach", J. Waterways and Harbors Div., ASCE Vol. 94,
No. WW3, August, 1968.
[3] Thomann, R.V., "The Use of Systems Analysis to Describe
the Time Variation of Dissolved Oxygen in a Tidal Stream",
A dissertation in the Dept. of Meteorology and Oceano-
graphy submitted in partial fulfillment of the require-
ments for the degree of PhD at New York University,
Feb. 1963.
[4] O'Connor, D.J., St. John, J.P., and DiToro, D.M., "Water
Quality Analysis of the Delaware River Estuary", J. Sanit.
Eng. Div., ASCE, Vol. 94, No. SA6, Dec. 1968.
[5] Saucedo, R., and Schiring, E.E., Introduction to Continu-
ous, and Digital Control Systems^. The Macmillan Company,
1968.
[6] DiStefano, J.J., Stubberud, A.R., and Williams, I.J.,
Feedback and Control Systems, Schaum Publishing Comp, 1967,
[7] Graves, G.W., Hatfield, G.B., and Whinston, A., "Water
Pollution Control Using By-Pass Piping", Water Resources
Research, Vol. 5, No. 1 Feb. 1969.
[8] Thomann, R.V., and Marks, D.H., Results from a Systems
Analysis Approach to the Optimum Control of Estuarine
Water Quality, For presentation at Third Conference of
International Assoc. on Water Pollution Resource, Munich,
Germany, Sept. 1966.
[9] Dwyer, P.S., and Macphail, M.S., "Symbolic Matrix Deri-
vatives", Annals of Mathematical Statistics, 19, 517-534,
1948.
[10] Smith, A.H., and Albrecht, W.A., Fundamental Concepts of
Analysis. Prentice-Hall, Inc. 1966.
[11] Winograd, S., "A New Algorithm for Inner Product", IEEE
Trans, on Computers C-17, 1968, 693-694.
-------
[12] Schaumburg, G.W., Water Pollution Control in the Dela-
ware Estuary, Submitted to the Dept. of Applied Mathe-
matics in partial fulfillment of the requirements for
the degree with honors of BA at Harvard College, May, 1967
[13] Frankel, R.J., "Economic Evaluation of Water Quality,
An Engineering-Economic Model for Water Quality Manage-
ment", Univ. of Calif. Berkeley, SERL Report No. 65-3,
Jan. 1965.
-------
Accession Number
Subject
Field & Group
06 A
SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Organization
Federal Water Quality Administration
Title
Mathematical Programming for Regional Water Quality Management
10
22
Authors)
Glenn W. Graves
Andrew B. Whinston
Gordon B. Hatfield
11
16
Date
8/70
12
Pages
123
Project Number
16110 FPX
21
,r Contract Numb er
Note
Citation
University of California, Los Angeles
23
Descriptors (Starred First)
'"'Regional Analysis, *Waste Treatment, *River Basin Development,
'•Optimum Development Plans, Estuaries, Simulation
25
Identifiers (Starred First)
Nonlinear Programming, Cost Minimization
27
Abstract
This report is an application of a non-linear programming
algorithm to the problem of optimal water quality control in an
estuary. The mathematical model that was developed gives the
solution to the general mixed case of at-source treatment, re-
gional treatment plants, and by-pass piping. The non-linear
algorithm is developed in considerable detail and a sample pro-
blem is worked out. Detailed results are presented for a pilot
problem to illustrate the method of solution. Actual data from
the Delaware Estuary is used to solve a large scale problem and
the solution is given. The results indicate that a regional
treatment system for the Delaware Estuary is superior, in terms,
of total cost, to other proposed schemes.
Abstractor
Dr. Glenn W. Graves
Institution
University of California, Los Angeles
WRjIOZ (REV. OCT. 1866)
WRSIC
SEND TO: WATER RESOURCES SCIENTIFIC INFORMATION CENTER
U a. DEPARTMENT OF THE INTERIOR
WASHINGTON, D.C. 20240
c U. S. GOVERNMENT PRINTING OFFICE : 1970 O • 401-599
------- |