Nonlinear Least Squares Analysis Of Oxygen Demand Data For Ultimate Demand And Deoxygenation Rate


Nonlinear Least Squares
Analysis of Oxygen Demand
Data for Ultimate Demand
And Deoxygenation Rate
May 1971
Robert L. Crim
Technical Report No. 4
Environmental Protection Agency
Office of Water Programs, MAR
918 Emmet Street
Charlottesville, Virginia 22901

-------
Introduction:
This paper describes application of the non-linear least
squares technique to laboratory test data for biochemical oxygen
demand. The final result is a method for determining the "best"
values for ultimate demand and the coefficient of deoxygenation.
Statement of problem:
Laboratory test data for BOD consists of values of the
total dissolved oxygen utilized (y^) after varying amount of time
(x^). This oxygen uptake is related to BOD by the expression:
L = a-y	(1)
where L is the BOD at any time
y is the total oxygen utilized by the sample from
time zero to the current time
ot is the ultimate oxygen demand of the sample
The standard equation for BOD is of the form
L = ctebx	(2)
where
x is time in the appropriate units
b is the coefficient of deoxygenation in units
of "'"/time.

-------
2,
Substituting (1) into (2) and rearranging gives the relation-
ship between measured oxygen uptake, time and the two unknowns
a and b:
y = a(l-ebx)	(3)
Clearly the classical least squares approach will be un-
productive because of the non-linearity of eq. (3). A bit of
mathematical juggling is required to eliminate the difficulties.

-------
Derivation of Algorithm
Given the equation
y = f (a ,b) = a(l-ebx)	(4)
where it is desired to determine a and b in the sense of least
squares.
Defining initial estimates of a and b as and bQ the
equation
y = f(a,b) can be written as
y = f(a +Aa,b +Ab)	(5)
J	o o
where
a = a +Aa	and	(6)
o
b = b +Ab.	(7)
o
Expanding in Taylor's series about c^and bQ gives:
y = £(a b) +	+ ....	(8)
0	0
where
(|£) = (l-ibx)	(9)
o
and
(|^-) = +axebx	(10)
o
evaluated at a and b
o	o

-------
4.
Substituting into the Taylor's series:
y = f(aQ,bo)+Aa(l-ebx)+Ab(+axeX)	(11)
The above equation is linear in Aa and Ab and the method of
least squares can be applied to determining Aa and Ab.
The function to be minimized is
N	2
£ {y -f (a ,b )-Aa(l-eboX)-Ab(a x^o*)}	(12)
K = 1
Taking the partials with respect to Aa and Ab and rearranging
gives the following normal equations:
N	N	2 N	u	h
I (y -f (a ,b )) (1-e oXic)=Aa £ (1-e oX<) +Ab j> (a xe oXk)(1-e oXk)
, K lv O O	i	. UK
K = 1	K=1	K=1
(13)
N	N
£ (yK_f !^a0,b0^ (a0XK6 °Xk) =Aa ^ °XK)(Ct0XK® °X|<)'
K=1	K=1
N	, 2
+Ab I (a x„eboXK)	(14)
i 0 K
<=1
let U = (l-eboXK)	(15)
<
= (a x eboX<)	(16)
K	OK
and
D = y -f (a ,b )	(17)
K K K 0 0

-------
5.
then the equations are written as:
N	N 2	N
I DkUk = Act I Uk + Ab [ UkVk	(18)
K=1	K=1	K=1
N	N	N 2
y D V =	Aa y UV +	Ab I V (19)
K K	L, K K	<
K = 1	K=1	<=1
Define:
AaCii + AbCi2 = ^13	(20)
and
AaC2i + AbC22 = ^23	(21)
where:
N 2
Cn = I Uk	(22)
K = 1
N
Cl2 = I Vk	(23)
K=1
N
Cl3 = I Vk	(24)
K = 1
N 2
c2l = C12 = I UKVK	(25)
K = 1
N
C23 = I Vk	(26)
K = 1

-------
Solution to the system may be done directly using Cramer's rule:
Aa =
cl3	+ C12
C23 + C22
C11 + C12
C2i + C22
(27)
and
Ab =
C11 + c13
C21 + C23
C11 + c12
C21 + C22
(28)
or
Aa = (+C13C22 ~ ^12^23)/(+Cl1^22" ^21^12)	(29)
Ab = (C11C23- C21C13(./(+Ci1C22- C21C12)	(30)
Refering back to the definitions of a and b:
a = a +Aa	(31)
0
b = b +Ab	(32)
0
Clearly successive values of a and b can be computed by
iterating through the procedure with improved otQ and b^ values
until the Aa and Ab terms are sufficiently small.

-------
7.
Initial estimates for a and b can be obtained directly from
oo
the observed data. Returning to eq. (13)
fi -bx.
y = a(1-e )
(33)
Take the derivative:
dy ,-bx
= abe
dx
(34)
Solve for a
dz
dx
a =
, -bx
be
(35)
Substituting (35) into (34) and rearranging gives the identity
dy -bx _ d_^ -bx
dx 6	dx 6
(36)
which will be used to estimate b for an initial b . This is done by
o
setting up the ratios:
&
dx
1 / "bx\
_ (e )l
dx
dx
. -bx.
(e )2
(37)

-------
8 „
and taking the natural log of both sides,
d
dx
In
= b[ (x), -(x)^]
dj
dx
(38)
By making the following definitions:
d£
dx
"» y ~y n
yn 'n-l
1 x -x
n n-l
(39)
'dj y2-yi
dx.
V
2 *2-*!
(AO)
(x)2 = x2+xi
(41)
(x) = X +x
I n n-l
(42)
bQ can be expressed in terms of the known data points as:
b = 2 In
o
(Vyn-lHx2"Xl)
(Wi)(y2-yi)
/[(x2+x1)-(xr+xn_1)]
(43)
and aQ can be approximated from (35)
as
a
= 1_
o 2
Wl
x -x
n n-l
-b
be
x +x
n n-l
-T- +
y2-yi
x2"xl
-b
be
x2+xi

-------
Subroutine BFIT
Subroutine BFIT will determine the coefficients A and B in the equation
-Bx
y = A(l-e ) in the least squares sense.
Required information:
Variable name
Description
A (returned)
B (returned)
N
EPS (returned)
IMAX
NERR
Ultimate BOD
Coefficient of deoxygenation
Number of x and y points being used.
The array of N values of time.
The array of N values of oxygen uptake.
Convergence criteria; 1/2% of current
standard deviation. If the difference
between successive values of standard
deviations is less then EPS, the system
has converged.
Maximum number of trials allowed*
If NERR=0 convergence was obtained to EPS,
If NERR=1 convergence was not obtained.
If NERR=2
(1)	A or B has reached 0.0
(2)	Error in input data;
y(N)
-------
10,
Form of call:
Any FORTRAN program can use BFIT with the statement:
CALL BFIT (A,B,N,X,Y,EPS,IMAX,NERR,ESQ): where the terms are as
defined above.
Both X and Y are dimensioned with 10 elements in BFIT. The
dimension statement must be the same size as the corresponding statement
in the calling program.
-------
11.
SUBROUTINE BFIT(A,B,N,X,Y,EPS,IMAX,NERR,ESQ)
DIMENSION X(100),Y(100)
Z=(((Y(N)-Y(N-1))*(X(2)-X(1)))/((X(N)-X(N-1))*(Y(2)-Y(1))))
IF(Z.LE.O.O) GO TO 120
BEST=2.0*ALOG(Z)/((X(2)+X(1))-(X(N)+X(N-1)))
AEST =((Y(N)-Y(N-l))/(X(N)-X(N-l)))/(BEST*EXP(-BEST*((X(N)+X(N-l))
1 /2.0)))
CEST= ( (Y (2) -Y ( 1))/(X(2)-X(1)))/(BEST*EXP(-BEST*((X(2)+X(1))
1 /2.0)))
B=BEST
A=(AEST+CEST)/2.0
DO 10 1=1,IMAX
C11=0.0
C12=0„0
C13=0.0
C21=0.0
C22=0.0
C23=0.0
DO 20 K=1,N
EXBX=EXP(-B*X(K))
UK=1.0-EXBX
VK=A*X(K)*EXBX
DK=Y(K)-A*UK
C11=C11+UK*UK
C12=C12+UK*VK
C13=C13+DK*UK
C21=C12
C22=C22+VK*VK
C23=C23+DK*VK
20 CONTINUE
DIV=C11*C22-C21*C12
IF(DIV) 5,120,5
5 DA=(C13*C22-C12*C23)/DIV
DB=(C11*C23-C21*C13)/DIV
A=A+DA*1.10
B=B+DB*1.10
ESQ=0 o 0
IF(B) 120,120,130
130 DO 30 K=1,N
FK=Y(K)-A*(1.0-EXP(-B*X(K)) )
ESQ=ESQ+FK*FK
30 CONTINUE
ESQ=SQRT(ESQ/FLOAT(N))
EPS=0.005*ESQ
IF(I-l) 50,40,50
50 IF(ABS(SAVE-ESQ)-EPS) 100,100,40
40 SAVE=ESQ
10 CONTINUE
NERR=1
GO TO 200
100 NERR=0
GO TO 200
120 NERR=2
200 RETURN
END
-------
12.
Examples of BFIT Results
The following curves are the output of a program utilizing
BFIT to examine BOD data. The solid line is the graphical
representation of the exponential curve, and the data points
are the points used in its derivation.
The runs were made on the USGS IBM 360/65 in Washington, D.C.
-------
CD-
CD
O
GO
CD
o
LU
Ultimate Demand= 10.5 mg/1
Deoxygenation Rate= 0.06/day
o
o
5.00
25.00
35.00
40.00
20.00
30. 00
15. 00
TIME IN DATS
90 9
6/15/70-
- RAPPAHANNOCK
-------
~
Ultimate Demand = 10.8 mg/1
~
Deoxygenation Rate - 0.041/day
>-
S. 00
10.00
^.00
15.00 20.00 25.00
TIME IN DATS
30.00
35.00
40.00
9314
7/16/70-
- RRPPRHPNNOCK
-------
~
~
~
Ultimate Demand = 3.9 mg/1
Deoxygenation Rate = 0.14/day
T
T
15.00 20.00 25.00
TIME IN DAYS
30.00
^3.00
5.00
10.00
35.00
40.00
100 1
4/23/70-
- RAPPAHANNOCK
-------
o_
o
o
oo
o
o
CO
LU
1 o
Q_ .
ZD=^
Ultimate Demand = 10.6 mg/1
Deoxygenation Rate = 0.056/day
o
o
35.00
30.00
.00
5.00
10.00
15.00
25.00
20.00
TIME IN DAYS
9018
8/ 5/70-
- RAPPAHANNOCK
-------
CO
o
CM
I CD
Q_
Ultimate Demand = 4.1 mg/1
Deoxygenation Rate = 0.064/day
o
o
10.00
.00
5.00
15. 00
20.00
25.00
30.00
35.00
40.00
TIME IN DAYS
100 4
5/12/70-
- RAPPAHANNOCK
-------
~
~
~

~
Ultimate Demand =8.4 mg/1
Deoxygenation Rate = 0.053/day

10.00
¦^b.oo
5.00
15.00 20.00 25.00
TIME IN DATS
30.00
35.00
40.00
9014
7/16/70-
- RAPPAHANNOCK
-------