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

-------