A systematic approach for reducing complex biological wastewater treatment models

Wat. Res, Vol. 31, No. 3, pp. 590-606, 1997

Pei'gllmon P I I : S0043-1354(96)00273-4

© 1997 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0043-1354/97 $17.00 + 0.00

A SYSTEMATIC APPROACH FOR REDUCING COMPLEX BIOLOGICAL WASTEWATER TREATMENT MODELS M. A. S T E F F E N S @, P. A. L A N T *@ and R. B. N E W E L L Advanced Wastewater Management Centre, The University of Queensland, Brisbane, Australia

(First received June 1996; accepted in revised form August 1996) Abstract--Biological wastewater treatment systems comprise a variety of processes which occur at vastly different rates. Biological growth, mass transfer, hydraulics and chemical reactions all occur simultaneously and are all inter-dependent. In this paper we address the question "to what extent can we de-couple these processes, and what are the associated issues? We aim to introduce people who work with biological wastewater treatment models to analytical tools which may be used for model reduction. We present a quantitative technique to compartmentalise states into fast, medium and slow. From this we have provided an algorithm for eliminating state variables from a model based on whether they affect the process in the selected "time scale of interest". Through the technique presented we provide a means of quantifying the interaction between state variables, the "speed" of a state and whether it is a candidate for reduction. A simple case study of a biological wastewater treatment process is presented. We were able to reduce four biological and 19 settler differential equations into algebraic equations. This resulted in significant savings in integration time. Application of the technique also highlighted the strong coupling between the slower biomass states and the rest of the model. © 1997 Elsevier Science Ltd. All rights reserved

Key words--activated sludge, ASM No. 1, model reduction, process dynamics NOMENCLATURE

S~ = S~ = X~ = Xs = Xh = Xa = Xp = Sno3= Snh4= So = S, =

soluble inert COD (mg COD L -~) readily biodegradable COD (mg COD L- ~) suspended inert COD (mg COD L -~) slowly biodegradable COD (mg COD L- ~) heterotrophic biomass (mg COD L -~) autotrophic biomass (mg COD L- ~) particulate decay products (mg COD L- ~) nitrate and nitrite nitrogen (mg N L- ~) ammonium (mg N L- ~) dissolved oxygen (mg COD L- ~) soluble organic nitrogen (mg N L -~)

X, = particulate organic nitrogen (mg N L- ~) X~k = alkalinity (mg HCO3 L -t) C, = settler states--solids concentration on layer i (mg SS L -~) E = maximum relative error

Greek letters 2 = eigenvalue (day- ~) z = time constant (days) W = slow mode reduction relative error

available information" (Dochain et al., 1995, p. 98).

INTRODUCTION With the increased complexity of biological wastewater treatment models, there is now a growing need for analytical tools which can examine process dynamics and coupling and which may enable model reduction. This was highlighted in a recent report published by the European commission C O S T 682 program (Dochain et al., 1995). The report identifies "systematic model reduction" as a future research need. "It would be interesting to develop a systematic procedure to reduce mathematical models o f the processes depending on the model objectives and

A "full" model is not always necessary for all users and all applications. As such, the cost associated with parameter identification, wastewater characterisation, model construction, validation and simulation may be substantially reduced if simpler, suitable models can be developed (Holmberg and Ranta, 1982; Jeppsson and Olsson, 1993). In this paper we introduce a systematic approach for reducing complex biological wastewater treatment models. The work presented here aims to answer two questions in relation to biological wastewater treatment processes:

*Author to whom all correspondence should be addressed [Fax: + 39 2 6447 4300].

1. Can the dynamics of the various system components be categorised on a time basis? That is, how do we divide a model of the process 590

Systematic reduction of complex models

591

Table 1. Qualitative estimates of the time scale of various activated sludge variables (Olsson and Jeppsson, 1994) Manipulated variable

States and outputs

lnfluent flow rate Waste sludge flow rate Return sludge flow rate Step feed Nitrate recycle Recycle flows or filter backwash Supernatant recycle (e.g. sludge digestor) Chemical dosing Carbon addition Air flow rate SBR cycle time

into a slow, medium and d e p i c t e d in e q u a t i o n (1)?

Time scale

Hydraulic retention time MLSS; SRT; Filament content Sludge blanket level; Recycle concentration; Effluent SS Flow rates; MLSS distribution Nitrate; DN rate; redox Flow rate; Effluent susp. solids

h days-weeks h~lay

COD; DN rate; P removal


Effluent susp. solids; P removal N removal; P removal DO cone; OUR; Oxygen transfer rate Outlet nitrate, ammonia, susp. solids


f a s t s u b s y s t e m as

dXslow

dt

h h
--fi(XFa~,,XModi.m,XSio., U)

where x represents variables describing the "state" of the process (state variables) and u the input variables.

dxras, = jq (XFast, XMedium,XSlow,U) dt

2. S u b s e q u e n t l y , c a n t h e s u b s y s t e m s b e decoupled a n d t h e c o m p l e t e m o d e l be t r e a t e d as s e v e r a l

dXM~lium--J~(XFast, XMedium,XSIow,U) dt

I'inear,q, Model I Eigenvalueto State 1 Association I ~,.,, = f~(x,~=,,xu~,,.,,Xs~,.,u) ~' | I X slow = A ( X Ftgt , X Medium, X stow , U )

~

(1)

ingularPerturbation-FastMode~ "Thereductionofstates | possessinglastdynamicsthatl quicklydieout" __j

S~gular Perturbation-Slow Mo~ | " T o reducestateswhich are I lappreximatelyconstant duringthe I L ~ timescale of interest" __j

Fig. 1. Systematic approach for model analysis.

592

M.A. Steffens et al.

0.3

'

'

State 1 '

'

'

x l 0.2 0.1

I

I

I

I

0.5

1

1.5

2

2.5

I

'

'

0 0

3

time State 2

0.2

'

'

'

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0.15 0.1 x2 0.05 0

0

i

I

P

L

i

0.5

1

1.5

2

2.5

3

time

Fig. 2. The response of the states to a unit step change in the input variable u. smaller models? In other words, can xeas,, xM~iu, and Xslowbe treated independently? One benefit which could result from "decoupling" the process in this manner is simplification of numerical solution. Models exhibiting a wide range of dynamics are often termed numerically "stiff" and are difficult to solve. Biological wastewater treatment models are inherently stiff due to the wide range of processes occurring. Several authors have attempted to ease the solution of such models by qualitatively "decoupling" the system as shown in equation (1) (Billing and Dold, 1988; Kabouris and Georgakakos, 1992). In this paper we use quantitative techniques to achieve the same end. Another major driving force for reducing biological wastewater treatment models, and to decouple states as shown in equation (1), is to assist with control system design. Olsson and Jeppsson (1994) recently attempted to qualitatively classify manipulated variables and states into time scales for the purpose of control system design (Table 1). It is our aim in this paper to extend this analysis via the use of quantitative tools, The work discussed in this paper relies on two mathematical techniques, namely eigenvalue to state association and the method o f singular perturbations. The latter technique depends on the successful application of an eigenvalue to state association algorithm. Eigenvalues are a mathematical property

of a model which provide modellers with valuable information on the transient response of the whole system. If they can be linked to the various model components (states), then the dynamics o f each state can be quantified. Singular perturbation techniques can be used to reduce model complexity by subdividing the model based on the response speeds of the different states. States with dynamics at either extreme (fast or slow) can be removed from model equations. It is then possible to generate a reduced order model applicable to the "time scale of interest". A summary of the systematic approach performed in this paper is given in Fig. 1. A simple carbon and nitrogen removing process employing the activated sludge model No. 1 (ASM No. 1; Henze etal., 1987)is used here to demonstrate the utility of these tools. Through this analysis we aim to determine benefits for simulation and control of biological wastewater treatment processes. Before we do this, however, it is important to introduce the reader to the concept of model eigenvalues and outline their relationship with process dynamics. EIGENVALUES AND PROCESS DYNAMICS

The objective of this section is to answer the question: "What are eigenvalues and how do they relate to process dynamics?"

Systematic reduction of complex models Biological wastewatertreatmentmodelsareusually presented in the form of a non-linear set of differential equations. Because the technique described here is only applicable to linear systems the model must be linearised (Appendix A), which results in the well known "state-space" model format: dx --

= Ax + Bu +

system is used here to illustrate the utility of eigenvalues (equation (3)). Consider: dx~ d--t-= -4x~ + 0.5x2 + u dx_~2= 1.3xl + 2x2 dt

W~w

dt

(3)

Or, in state space notation:

y=Cx+Du+W2w (2) where x is a vector of variables reflecting the systems state, called state variables (e.g. reactor ammonium concentration); u is a vector of input variables we have control over (manipulated variables, e.g. air flowrate); w is a set of input variables we have no control over (disturbance variables, e.g. inttuent flowrate); and y a vector of outputs or measured variables (e.g. efltuent ammonium concentration). It is important to remember that the linearised model will change depending on the operating point that the model is linearised about. Eigenvalues are a well known mathematical property of a matrix (Appendix B). In this work the eigenvalues of the A matrix in a linear model (equation (2)) are analysed and subsequently used to help describe model dynamics. A simple, two variable

dx dt - A , x + B,u

,

,

,

where

A, = [ - 4 1.3

x =

[x,] xz

0_.52]

[10] and

B~ =

(4)

The eigenvalues of A~ (2~ = - 4 . 2 8 4 5 + 0 i , 22 = -1.7155 + 0i) are both negative with zero imaginary components. The response of the states to a unit step change in the input variable u, is given in Fig. 2. The response of both states is smooth up to a new steady state value. This behaviour is common to systems with eigenvalues possessing zero imaginary components. Another two-variable system can be used to examine the response of a system with non-zero imaginary eigenvalue components:

State 0.8

593

~

1 ,

,

~

,

0.6

0.4

xl

0.2

. . . . . .

0

-0.2 0

I

I

I

I

I

2

4

6

8

10

J

I

I

14

16

18

i

i

i

;

i

i

,

,

,

I

I

I

I

I

10

12

14

16

18

State 1.5

-

1

I

12 Time

20

2

-

x2 0.5

0

0

I

I

I

I

2

4

6

8

time

Fig. 3. The response to a unit step in the input, u.

20

594

M . A . Steffens et al. i

-1.5

i

i

i

i

i

i

i

i

-2

"~ -2.5 C Q 0

m Q.

I~ -3.5

'!4.5

0

I

I

I

I

I

I

I

I

I

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Homotopy Parameter (r)

1

Fig. 4. Eigenvalue traces for a coupled two variable system.

. 0 2°

~ ~,,,,

~

i

i

i

]

r

~

,

,

-0.22

~ ~2

•-0.24 "~ -0.26

._=-o.2s Il

.o.3

I-0.32

//

-

@ -0,~

-0.36

/

-0.38

/ /

-0.4

0

~'1 I

0.1

I

0.2

I

0.3

I

0.4

I

0.5

I

0.6

I

0.7

Homotopy Parameter (r) Fig. 5. Eigenvalue trace.

I

0.8

I

0.9

1

Systematic reduction of complex models

responseOSCi I~ ~/~/~,~ lat°ry" '~=~ ~ 1 ~k

E

H(r) = ( 1 - - r ) * A D + r * A ~

0~
Fastresponse l~ ----1

l-J ~

"~ [

Re(;~) "~ ~¢

~

595

~

Unstable response

(6)

where H is the Homotopy matrix; AD is the matrix representing the decoupled equations (diagonal (A)); A is the linearised model matrix; and r a parameter changed to achieve a linear progression from decoupled to coupled system. For the first example (equation (3))the decoupled system is:

T

response I V

v

I

d~ =

-4x,

Fig. 6. Transient response as a function of eigenvalue (2)

dx2- -2x2 and

magnitude,

A2 =

I -°' ']. 1.3

-0.2

21 = - 0 . 3 + 1.358i

heei.enva,uosofthodocoup,odsystemhavo

2= and

22 = - 0 . 3 - 1.358i

(5)

The response to a unit step in the input, u, is given in Fig. 3. On comparison with Fig. 2 several important pieces of information should be noted: e Non-zero

--02](7)

imaginary eigenvalue components

indicate oscillatory response, • Both variables in Fig. 2 decay to steady state at approximately the same time. In other words, the dynamics of both variables are of similar magnitude, • The dynamics in Fig. 3 are significantly slower, This is reflected in the smaller real eigenvalue components in the second system (equation (5)). At this stage of the analysis, the dynamics of the overall system (both states) can be described by the eigenvalues. If each eigenvalue can be linked to a state, then the dynamics of each state can be quantified. The mathematical homotopy method used to do this was proposed by DeCarlo and Sacks (1979) and Wasynczuk and DeCarlo (1981) and later refined by Robertson (1992). It involves starting with a system with an obvious relationship between eigenvalues and states, such as a diagonalised A matrix, or decoupled system. This is then transformed into the actual system while tracing the eigenvalues. The transformation is described by equation (6).

Sludge recycle Fig. 7. Process configuration for the case study.

obvious state association (x~ is associated to 2~ = - 4 and x2 is associated to 22 = - 2 ) . The homotopy matrix (/-/) can now be formed (equation (6)) as r varies from 0 to 1 and the eigenvalues at each r calculated. If the eigenvalues are traced as coupling is introduced (as r changes from 0 to 1), the state with the strongest link to each eigenvalue can be

determined (Fig. 4). The traces show that 2~ can be associated with x~ and 22 with x2. From this we can conclude that xt will respond to transients more quickly than x2 as the eigenvalue associated to xt is more negative. This behaviour is observed in Fig. 2. The eigenvalue to state association for this system is relatively obvious as the eigenvalues change minimally as r varies from 0 to 1. The second system (equation (5)), however, is not as simple. The eigenvalue trace is presented in Fig. 5. The real parts of the eigenvalues converge quickly to the same value as coupling is introduced. This behaviour is indicative of a system with a large degree of coupling. As was discussed above, similar real parts of a system's eigenvalues implies dynamics of the same magnitude. The transient response of the two variables presented previously (Fig. 3), shows this to be the case. In summary, the eigenvalues of the A matrix provide the modeller with a large amount of information pertaining to its transient behaviour and are often termed the system modes. Several important properties of the system eigenvalues in relation to process dynamics include: • positive real components imply instability; -large imaginary components indicate oscillatory dynamics; and • large negative real components imply a fast decay to steady state (desirable). These are summarised in Fig. 6.

596

M.A. Steffens et al. Table 2. Wastewatercharacteristicsand processdesignparameters Wastewater component Value Design parameter Value S, 120 Aerobicvolume 4000 m3 X, 300 Clarifierarea 1000m2 S~ 60 SRT 12 days X~ 120 Influentflowrate 4000 m3day-' Xh 0 Recycleflow (ratio 1:1) 4000m~day-' .!". 0 Clarifierdepth 4m S.M 40 S.o3 0 x. 5 s~ 5 so 0 xp 0 s.~k 260

CASE STUDY DESCRIPTION

where F is the influent flow rate and R(S,) is the rate

Of the many activated sludge process configurations available, the single aerobic basin and secondary settler in series, achieving carbon removal and nitrification, is the simplest (Fig. 7). Anoxic process dynamics can also be evaluated using this configuration by manipulating the dissolved oxygen setpoint. To simplify subsequent analysis, the study will be limited to this configuration. The widely accepted activated sludge Model No. 1 (Henze et al., 1987) is used here to simulate the biological phenomena in the aerated basin. A total of seven dissolved and eight particulate components are used to describe the wastewater and sludge in the activated sludge Model No. 1 (ASM No. 1; Henze et al., 1987). Each of the 13 components are then formulated as mass balances around each reactor (one in this case). For example, the mass balance for readily biodegradable COD is: (8) dSs = ~F ("Ss , , . - Sso.t)+ R(S,) •

term for the readily biodegradable COD, defined by the ASM No. 1. The rate term is made up of any biological or physical (e.g. mass transfer) processes affecting the component. Thirteen mass balance equations such as this make up the biological model. Typical influent wastewater characteristics used in subsequent simulations are given in Table 2 along with process design parameters. Default kinetic parameters for 20°C were taken from Henze et al. 0987). The settler model utilised in this work is based on the one dimensional flux theory concept (e.g. Vitasovic, 1989) with the settler divided into 20 layers. The model is made up of a solids mass balance on each layer, which results in 20 additional state variables; each state being the solids concentration in the layer (Cj - (720). The NIMBUS simulation package (Newell and Cameron, 1991) is used throughout this work for simulations and to generate linearised models. The package utilises a diagonally implicit Rungc-Kutta

.01 0""0"'0"'0

0 0 ......

•0. t~,..

•1

. 0 -O..O..o..,O...O,.

"')t

""""'" ,

.

0--.@'@"

..,..x,

t~

- - v - - Xs ..,.. x, - - o - - xa

Q

°°~ll Xp

,

.......

- - e - - Si

" i..i..m.-u---m"l'"

:"

::

- - 0 - - Xno3 > '.~) I.IJ

• , ~ . - Xni14

lO &- - - & - . . & . - & - . A . . i - . ~lk...b..&--& -. A- . i . • ~ - - - b . - & . . & - - , L - .A--,&---&...

::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 100 0-..0"--0--0--0-.0.1000

0-- 0--0-.0..0--

0..0-. 0--0-.0.-0--0--

0-- 0--(

, 0 H o m o t o p y p a r a m e t e r (r) Fig. 8. E i g e n v a l u e traces for the A S M N o . 1 M o d e l variables (negative values).

- - ~ - . Sn - - e - - Xn

Systematic reduction of complex models

597

Table 3. Associated eigenvalues (real part) and a p p r o x i m a t e time constants for the single basin model Steady state 1 (feed flow = 4000 m 3 d a y -~) eigenvalue ( d a y - ' ) z(h)

Variable Co C7' C2 C4 C~ C, C,8 C,9 C~ S~ S, C2o Xn X~ S,~

- 1.35E - 1.31E -8.17E -7.03E -5.89E -5.06E -4.29E -4.29E -3.13E - 3.09E - 7.38E -6.85E --6.85E --5.02E -3.51E

+ + + + + + + + + + + + +

03 03 02 02 02 02 02 02 02 02 01 01 01 01 + 01

0.017 0.018 0.029 0.034 0.040 0.047 0.055 0.055 0.076 0.077 0.325 0.350 0.350 0.477 0.684

+

- 1.00E + 00

S~o3 S, X,,k

--3.15E - 2.56E - 7.38E - 6.88E

Xi X,

--

- 1.29E - 1.25E -8.72E - 5.99E -6.98E -5.27E -4.00E -4.00E -2.85E -4.66E - 1.10E -5.24E -6.87E - 5.03E -5.52E

23.915 24 24

- 1.00E + 00 -- 1.00E + 00

Xh X.

Steady state 2 (feed flow = 6000 m 3 day -~) eigenvalue ( d a y - ' ) ~(h)

01 01 02 02

+ + + + +

03 03 02 02 02 02 02 02 02 02 02 01 01 01 01

0.019 0.019 0.028 0.040 0.034 0.046 0.060 0.060 0.084 0.051 0.217 0.458 0.349 0.477 0.435

- 1 . 5 1 E + 00 - 1.50E + 00 1.50E + 00

15.943 16.000 16

+

+ + + + + + +

+ +

-

76.152 93.923 325.023 348.779

-3.13E - 2.54E - 7.43E - 6.56E

-

01 01 02 02

76.784 94.440 322.993 366.070

~Settler states Cr.-C~7 are not shown as they have the same eigenvalues as (?7.

(DIRK) integration algorithm (Cameron, 1983) which is particularly suited to stiff systems. All computations were performed on a 486-100 MHz desktop computer,

RESULTS ANDDISCUSSION Prior to the commencement of the dynamic analysis work, problems with slow integration time and solver stability were experienced. Dissolved

Snh4

Xa

0.00

0.10

5oo I t" •

'

~e

•

0.05 +

~-10.oo '

~ 0.oo T

-0.05t =-zo.oo

~.~-25.00

| -15.00 .oo

"

•

_o101 .

•

"~-0.25 -0.151~ i

•

-~.oo

"

IZ: -0.20 ]_

,.

-40.00

•

-45.00

,

0.00

5.00

i

i

10.00

15.00

Time(days)

•

•

-0.30

20.00

0.00

-1.0010

•

t

2.00

4.00

6.00

8.00

Time(days)

10.00

Xh •

@

•

•

-38.64

--~.ee

~,ol 0030 "

•

~-1.0040 .~-1.0050

•

,~-38.70 "= -38.72 ~-38.74 ~"-38.78

•

.

• • t

0.00

•

I

2.00

•

•

•

~ -38.68

•

-1.0100 -1.0110

•

0

-38.82

~-I.OO~

"8~,-1.0060 "~ -I.0070 ,~ -1.00a0 -1.0090

•

,

Sno3 -1.0000

•

•

I

I

I

4.00 6.00 Time (days)

-38.78 -38.80

I

8.00

10.00

• • • q¥O•

0.00

I

2.00

I

I

4.00 6.00 Time (days)

I

8.00

F i g . 9. E i g e n v a l u e t r a c e t h r o u g h a f e e d a m m o n i u m s t e p a t t = 1 d a y : 4 0 - 1 0 0 m g N H 4 - N L -~ (all i m a g i n a r y p a r t s = 0).

10,00

598

M.A. Steffens et al.

1.50'

1.45

~

1.40

~1.35

~

]

- * ' - Initial

-~-

1.3o

Initial

influentNH4=40 mg N/I influentNH4=100m~ N/I

)1.25 ,~r 1.20

~

1.15

~

1.10

1.00

I 25

24

I 26

I 27

Time (hours)

I 28

t 29

30

Fig. 10. Simulated ammonium concentration (12.5% step in feed NH4-N at t = I day). oxygen (DO) kinetics were included in the model with a PI controller used to regulate the DO at its setpoint, DO dynamics, being one of the quickest in the process (Olsson and Jeppsson, 1994), were thought to be the major cause of the numerical instability, DO was, therefore, approximated as a constant. In other words, the DO control was considered tight enough to warrant removing the oxygen mass balance. This "knowledge based" model reduction step drastically improved solution time and solver stability, highlighting the fact that significant real benefits can be gained from temporal decoupling of this system. All subsequent modelling work presented in this paper is without the DO mass balance, The model has 32 states, 12 of which are from the ASM No. 1 ( S o assumed constant) and 20 from the settler model. The NIMBUS simulation package was used to generate the linearised models (Newell and Cameron, 1991). Figure 8 depicts the eigenvalue

Sno3 1400

°.

~1t2°° ~_1000 ~a00 °

600-

j400. 2oo 0

:

~A

,,

2

I

P

411me (~y~) 6

l

8

10

Fig. 11. Eigenvalue trace when DO is reduced to zero (imaginary part = 0).

traces for the biological states in the aeration basin, from a decoupled (left) to a fully coupled system (right). The non-linear model was linearised using the process data given in Table 2. Although only the real part of the eigenvalue has been plotted in Fig. 8, every variable has a zero imaginary component. This fact infers a smooth (non-oscillatory) response to transients (Fig. 6). This is reflected in the simulations which are discussed later. In order to provide a more physically meaningful parameter, we "transform" the eigenvalues into equivalent "time constants". This is a simple extension of the definition of a time constant (Q for a single input-single output system (Appendix C) w h e r e , = - 1/2 and has units of time. For a linear, first-order system, z is the time taken to get to 63% of the steady state value after an input disturbance. The time constants are, therefore, direct measurements of the speed of response of the individual state variables. Table 3 shows the eigenvalues and time constants (Q for each of the states at two operating points. It summarises the information gained from the eigenvalue to state association technique. The eigenvalues correspond to the right-hand side of the eigenvalue traces in Fig. 8. Wastewater characteristics and process parameters defined in Table 2 w e r e used to generate the linear models. • The state variables can be classified, on a time basis, into three groups: the fast states, with time constants of the order of 1-10 rain. These include all settler states, S,, Xs, S, h4, .t", and S.. Oxygen,

Systematic reduction of complex models

Top settler layer

599

Settler Layer 6

1

1

A

A

O.S

~ I11

aS

t=O days ._:

~,~ ._m

t=-2 days >

o

m ~

t=O days

t=2 days >

o

ID v

m

E

0.5

m

.o.5

E

-1 -900

-o.s

-1 -800

-700

-600

-500

-400

-1360

Real(eig) (days)

-1340

-1320

-1300

Settler Layer 11

Bottom Layer 0.1

I=-2 days

-,,o

~

•. ~

i j..1~-

~,

o.os

>~ m

/

t=2 days

'~

f

"~

-... -0.05.

m

E

-1260

Real(eig) (days)

•.435

--

-1280

-45o

E

t=-Odays

-.

.o.1

.t=O days

-455

-I040

-1020

-1000

-980

-960

-940

Real(eig) (days)

-150

-100

-50

Real(eig) (days)

Fig. 12. Eigenvalues associated to settler layer concentration traced through a feed flow step (4000-4/000 m3day-t).

not included here, would also fit into this category; medium states including nitrate (changes significantly), soluble inerts and alkalinity; and slow states including Xa, X~, X~and Xp. • Due to the non-linearity of the model, we can expect that the eigenvalues will vary depending k~t~

__

~.1=-1/xI

L2=l/'t:2

'~

Imag(X)

"~etlt,slow

~ / S ~¢~

/[~/,~ /~/~/: //;~ ,,/, ~' I~L

~f//;~ ~,/

The points to note from Fig. 9 are:

/

• The real part of the eigenvalue associated to

! :"~~/ J

This ammonium ammonium and and autotrophic biomass (X,) significantly (Sn~4) decreases subsequently increases as the transient propagates through the system. implies dynamics are slowing (time constant gets larger), then speeding up. To confirm these results two transients were simulated at the two different steady states (i.e. steady state at a feed NH4-N of 40 mg L - ' and 1 0 0 m g L -~ ). The feed NH4-N was stepped by 12.5% in both simulations (Fig. 10). The simulation at a feed concentration of 100 mg N L-: is significantly faster, a fact that is reflected in the eigenvalue trace.

~S S~ ~ ;:-,.~

~;,~ j Real(X) • , , .,0/,~

Cdtlealregion I I ~-;/y ~ ~ ~. :" ~: "limescale of interest Fig. 13. Singular perturbation eigenvalue map. C/,

on the operating point. The magnitude and direction of such variations should be examined. Figure 9 shows the real part of the eigenvalues (in days) of several biological model states after a step in the feed ammonium concentration from 40 to 100 mg NH4N L-L A relatively severe disturbance is used to ensure excitation of the systems dynamics.

600

M . A . Steffens et al. Table 4. Fast mode reducible states (Cfftop)-C:0 are settler layer concentrations) Maximum relative error (E) (%)

Scenario 1 2 3 4 5 6 7

Time scale of interest (~] = -1/2,)

5 5 5 5 5 10 10

1h 2h 4h 12 h 1 day 1/2 h 2h

The eigenvalue associated to nitrate stays relatively constant and matches the hydraulic residence time of one day. In other words, the nitrate dynamics are not significantly affected by the feed ammonium concentration, • Several other eigenvalues, such as that associated to X., do not change significantly throughout the transient. This lack of sensitivity is not surprising as ammonium concentration should not affect COD kinetics a great deal (assuming adequate DO). •

In order to investigate anoxic process conditions, the DO setpoint was reduced to zero and associated eigenvalues traced as DO is reduced from 1.5 to

S,, Ci-Ci9 S~, C~-C]9 S,, C~-C2o, S~, X~ S~, C~-C2o, S~, X~,, Snb,,, Xs S~, C~-C2o, S., X., S.h,, X, S,, C]-C,9 S,, C~-C2o, Sn, Xn

0 mg L-] (at t = 1 day). Figure 11 shows the real part of the eigenvalue associated to Sno3. The simulation begins at steady state 1 (Table 3). As the oxygen is depleted, the reactor becomes anoxic and nitrate dynamics are controlled by denitrification. This is reflected in Fig. I l by the large increase in the eigenvalue associated to nitrate (time constant gets smaller). This is to be expected as heterotrophic biomass activity controls denitrification. The eigenvalue then decreases again as nitrate is depleted and the environment changes to anaerobic. Figure 12 displays eigenvalue traces for some of the settler state variables for a feed flow step from 4000 to 6000m3/day -~. Flowrate through the settler is

Ss

4.5

States which can be fast mode reduced

Xs

40

4

35.

O 2.5

10 0

5

10

15

20

0

Time (days)

10

15

20

tS

20

Time (days)

Sno3

6O

5

Snh4

2Q

~,,

15

g, 0

0

'

'

'

5

10

tS

Time (days)

0

20

0

10

Time (days)

Fig. 14. C o m p a r i s o n o f full (. - -) and fast m o d e reduced ( - - ) models (all plots are effluent concentrations).

Systematic reduction of complex models

601

Sno3

Snh4

100

10

E

eo c

o

e

.o_ .

¢~

, 2O

Q 0

4

U

2 0

2

4

6

8

10

0

Time (days)

2

4

6

8

10

Time (days)

Fig. 15. Ammonium step 40--* 100 mg L-' in feed (reduced (--), full (.-.). expected to have a significant effect on dynamics. These plots show this to be the case. The real part of three out of four eigenvalues increases significantly, which implies a slowing of dynamic response through the transient. Conversely, the bottom layer "speeds up" as flow is increased. This is due to the modelling assumption of zero flux due to gravity exiting the bottom layer, thereby making this layer more dependent on hydraulics, Although the dynamics do change, they are still significantly faster than many of the ASM No. 1 state variables throughout the transient. Information such as this will be useful when identifying which states can be reduced,

The method of singular perturbations The method of singular perturbations is a model reduction technique which can be applied to models displaying a wide range of dynamics. In this section we aim to reduce the model using the singular perturbation technique based on the eigenvalue to state association described above. Major contributions to this field were made by Kokotovic et al. (1976), Kevorkian and Cole (1980) and later by Robertson (1992). The general form of non-linear differential-algebraic-equations (DAE) found in process systems can be represented as follows:

/~crit,fast

/

./

//

~,]=-1/x]

Imag(;~)

.///

/

Fast modereducibleregion \

//

/

/

/

)

/

I Real(Z)

~idays'1 [

-30

" - ~

-25

-20

X

-15

~,.~ at t=0

'

-5

X

' -~

5

days

Zn~ at t=2 days Fig. 16. Movement of the eigenvalue associated to ammonium through a feed ammonium step.

M.A. Steffens et al.

602

Layers top to 3 (2 above the feed)

Layers 4 to feed (5)

..., 0.025

~ 20

E 0.02 -.-,.o

.-. "" f~

.~-'~ - - ~

E ~'15 .o

~"

~ "~"• ~"~

t

0.015

I

0.01

o°

I

~ 10

i

o c

__1

o =

0.005

"o

o 03 - -

~

0

i

i

5

o 03

_ 0

~] _~

,

2 4 Time (days)

- - - -

6

i

0

Layers feed(5) to 17

i

2 4 Time (days)

J

6

Layers 18 to bottom(20)

._, 600

,_, 8000

E

~

"="500 0

~'~"

E

-" ~ > "

"~"6000 0

~- ~'~'~'~

~ ~ "~ ~

t

o

400

e '~ ~

c

~ 3 0 0 _ _,"g

'~ 4000 o

)

o 03

°

o 2000 ~ "g o 03

20E0

i

i

i

2

4

6

,,,

•

0

Time (days)

i

i

=

2

4

6

Time (days)

Fig. 17. Reduced order model (.-) compared with full model (--) response to feed flow step. dx --~ = f ( x , u, w)

y = g(x, u, w)

(9)

where x, u, y and w are defined as in equation (2). However, the singular perturbation technique requires the model equations to be in the following explicit form:

set of variables (Xs) are slow enough to be assumed constant. The reduced order model is formed as follows:

0 =fi(XF, XM, XS, U, W) dx =f2(x~,x~,xs, u, w) dt dxs _ 0 dt

dxF d--t-= Jfl"~-F , XM, XS, U, W)

dxM

y

dt --J~(XF, XM, XS, U, W)

_~t=fi(XF,

XM,

XS, U, I¢)

y = g(xF, xu, Xs, u, w)

(10)

where xF, xu and Xs represent groups of states classified according to response time (fast, medium and slow). This classification may be easily extended to any multiple of time groupings, The reduction technique assumes the first (XF) set of variables reach steady state quickly and, therefore, the derivative term is negligible. In other words, the fast subsystem is converted from differential to algebraic equations by setting the derivative to zero. The variables represented by xM are irreducible as their dynamics are in the range of interest. The third

=

g(xF, xM, Xs, u, w)

(11)

While the singular perturbation theory is conceptually simple, problems arise when the user has to convert from the implicit to the explicit format. In other words, how to decide which states are fast, medium and slow. This is where we intend to exploit the eigenvalue to state association results. To do this a critical region between 2C,,.F~,, and J~Crit,$1ow is introduced in the eigenvalue space to define the regions of "reducible" and "irreducible" states. Any states with associated eigenvalues outside the critical regions are considered reducible. Figure 13 depicts the reducible regions and the time scale of interest in the eigenvalue space. Fast mode reduction. The aim of this section is to identify which state variables can be reduced using the following approximation: 0 =J~(xF, xM, xs, u, w). The problem is to decide how large the eigenvalues

Systematic reduction of complex models have to be to qualify for fast mode reduction or, in other words, what is 2Cr~,,Fast.This decision is dictated by model accuracy requirements and the lower bound on the "time scale of interest" (2~ = - 1/~). Before proceeding we must clearly define the error incurred when reducing model order: Maximum relative error = E =

603

The fast mode reduction procedure, then, is as follows: 1. Choose the time scale of interest and maximum relative error for the reduced order model. 2. Calculate the critical eigenvalue, 2Cr~t.Fast

(12)

(equation (13)). 3. Associate eigenvalues to states and determine which states fall in the fast mode reducible region.

where x refers to any state variable and the maximum is evaluated over the full integration time.

4. Convert the fast mode reducible differential equations to algebraic form (equation (1 l)).

Maximum(xF,,Mode=- XRedu~Moa=~) x~.. Mode~

A relationship between fast mode reduced model accuracy requirements and t~,Crit,Fast was defined by Robertson (1992) using a Monte-Carlo simulation, The relative error was calculated for a large number of systems and correlated to the critical eigenvalue (~,Crit,Fast). 2Crit,Fast :

2~*0.5 E

(13)

Table 4 summarises the fast mode reducible states for the biological wastewater treatment model as a function of E and the lower bound on the time scale of interest, z~. It is interesting to note that all but one of the settler states can be reduced in all of the situations encountered in Table 4. This, however, will change with flow rate and plant design. Readily biodegrad-

Snh4

Sno3

15 - -

70

E

E - - 50

5

= 40 0 to 3O

0 0

0

4.s

20

s 5.5 Time (days)

6

4.s

'

s s.s Time (days)

Xh

Xa

844

50

4o

o tO 842.5

o 35 tO

''''''

3O 842 4.5

5

5.5

Time (days)

6

4.5

' 5

' 5.5

Time (days)

Fig. 18. Slow mode reduced model (--) compared with full model (-.-).

6

604

M.A. Steffens et al.

able COD, &, is also considered fast enough for reduction in all cases. This is not surprising as the real part of its eigenvalue is larger than all other biological states by an order of magnitude (Table 3). Scenario four, presented in Table 4, will now be examined in detail here. The outcome of the fast mode reduction technique using a relative error of 5% and lower bound on the time scale of interest of 12 h (zt) is that S~, X~, Xn, Snh4 and CrC2o can be approximated with algebraic equations. Rigorous model comparison requires excitation of the various system dynamics. To do this, several influent step disturbances were introduced in series (at t = 0.5 days feed TCOD 600-1200mgL-~; at t = 5 days feed NH4-N 40-80mgNL-~; at t = 10 days feed flow 4000--6000 m 3 day-~; and at t = 15 days DO setpoint 1-2 mg L-~). Figure 14 compares the outputs of the "reduced" and "full" models in response to the disturbances, From Fig. 14 it can be seen that fast mode reduction does not incur any steady state modelling errors. The reduced order model matches the full model response for all variables except ammonium (S, h4). Ammonium concentration displays the largest relative error which exceeds 100% at 5 days. To further examine this phenomenon, a feed ammonium step (from 40 to 100 mg N L -~) was simulated over l0 days (Fig. 15). The ammonium concentration predicted by the reduced order model again incurs a relative error greater than 100%. This is clearly contradictory to the fast mode reduction recommendation. Reasons for the errors in the ammonium concentration are due to the non-linearity of the model and are reflected in the eigenvalue trace presented in Fig. 9. As feed ammonium is increased, the ammonium state "slows down" significantly, moving from the reducible into the irreducible region (Fig. 16). This suggests that ammonium is not a good state for fast mode reduction. However, all other fast mode reductions result in minimal errors and can be considered successful. Selected settler state responses to a step in feed flow (150%) are presented in Fig. 17. The fourth scenario in Table 4 is again used here. Both the reduced order and complete model predict similar behaviour with a maximum relative error of less than 5% observed towards the end of the integration, With all four biological and 19 settler states reduced, integration time is significantly less for the reduced model (approx. half the time for a 10 day simulation). It is also interesting to note that the number of integration steps used by the reduced order model (41) is significantly less than that used by the full model (212). More savings may be seen when simulating multiple basin configurations, as are common at present. Slow mode reduction. The aim of this section is to identify which state variables can be reduced using the following approximation: dxs/dt = 0.

If a state is identified as being slow mode reducible it may be approximated as constant for the duration of the simulation. The crux of the problem, as with fast mode reductions, is to identify how slow a state must be before it can be reduced. This is determined by approximating the maximum difference between the full and reduced model solutions, W,, for each state. For details of the calculation see Appendix D. The upper bound on the time scale of interest, z2, must be specified to calculate W~. As with the fast mode reduction work, a critical tF is empirically related to the maximum relative error for any state using a Monte-Carlo simulation (Robertson, 1992): Wc~t = Maximum relative error 10

(14)

The slow mode reduction technique is summarised as follows: ' 1. Define values for the maximum relative error and time scale of interest (z2). 2. Calculate udcdt using equation (14). 3. Calculate W~ for each state. 4. If W~
Systematic reduction of complex models treatment systems. The eigenvalue to state association technique enables us to quantify the "speed" of a state. This information is then used for reducing the model via singular perturbation analysis. The singular perturbation technique enables us to " r e d u c e " state equations which are faster than the "time scale of interest" to algebraic equations and those which are slower to constants, Application of the technique to our case study resulted in four biological and 19 settler state variables being identified as fast mode reduction candidates. That is, 23 differential equations could be approximated as algebraic equations without incurring significant modelling errors. Significant reductions in numerical solution time were achieved using the reduced order model. It was interesting that we were unable to "slow m o d e " reduce any state variables until the time scale of interest was very short. This was due to the state interaction prevalent in the biological model. Because o f this, we were not able to eliminate any parameters from the model,

Acknowledgements--The authors would like to acknowledge both the CRC for Waste Management and Pollution Control Ltd and an Australian Postgraduate Award (for Marc Steffens) for funding the work presented here.

REFERENCES Billing A. E. and Dold P. L. (1988) Modelling techniques for biological reaction systems 3. Modelling of the dynamic case. Wat. SA 14(4), 207-218. Cameron I. T. (1983) Solutions of DAE Systems using DIRK methods. IMA. JNA 3, 273-289. DeCarlo R. A. and Saeks R. (1979) A Root Locus Technique for Interconnected Systems. I.E.E.E. Trans. Sys. Man. Cyb. SMC-9(1), 53. Dochain D., Vanrolleghem P. and Henze M. (1995) COST 682--Environment. Optimizing the design and operation of biological wastewater treatment plants through the use of computer programs based on dynamic modelling of the process. European Commission Report, 1992-1995. Luxembourg. Henze M., Grady C. P. L. Jr., Gujer W., Marais G. v . R . and Matuso T. (1987) A general model for single-sludge wastewater treatment system. War. Res. 21(5), 505-515. Holmberg A. and Ranta J. (1982) Procedures for parameter and state estimation in microbial growth process models, Automatica 18, 181-193. Jeppsson U. and Olsson G. (1993) Reduced order models for on-line parameter identification of the activated sludge process. Wat. Sci. Technol. 28(11/12), 173-183. Kabouris J. and Georgakakos A. (1992) Accounting for different time scales in activated sludge process control, Wat. Sci. Technol. 26(5/6), 1381-1390. Kevorkian J. and Cole J. D. (1980) Perturbation Methods in Applied Mathematics. Springer, New York. Kokotovic P. V., O'Malley R. E. and Sannuti P. (1976) Singular perturbations in order reduction in control theory--an overview. Automatica 12, 123-132. Newell R. B. and Cameron I. T. (1991) NIMBUS Users Manual. The University of Queensland, Brisbane, Australia. Olsson G. and Jeppsson U. (1994) Establishing cause-effect relationships in activated sludge plants--what can be

605

controlled? Proceedings of the Forum for Applied Biotechnology (FAB). Brugge, Belgium. Robertson G. A. (1992) Mathematical Modelling o f Star tup and Shutdown Operation of Process Plants. Ph.D. Dissertation, The University of Queensland, Brisbane, QLD, Australia. Ross S. L. (1974) Differential Equations, 2nd edn. Wiley International, New York. Vitasovic Z. Z. (1989) Continuous settler operation: a dynamic model. In Dynamic Modelling and Expert Systems in Wastewater Engineering (Edited by Patry G. G. and Chapman, D.), pp. 59-81. 'Lewis, Chelsea. Wasynczuk O. and DeCarlo R. A. (1981) The component connection model and structure preserving model order reduction. Automatica 17(4), 619.

APPENDICES Appendix A--Model linearisation

Biological wastewater treatment models can be generically represented as follows: dx d-7 = f(x, u, w)

(15)

where x is a vector of variables reflecting the systems state, called state variables (e.g. reactor ammonium concentration): u is a vector of input variables we have control over (manipulated variables, e.g. air flowrate); w is a set of input variables we have no control over (disturbance variables, e.g. influent flowrate); and y a vector of outputs or measured variables (e.g. effluent ammonium concentration). The linear model is formed by numerically evaluating the change in function values ( f and g) resulting from small changes in model variables (x, u and w) at any particular operating point. _ Of . ~ = A - tgx oP' B = ~- oP; W~ O_~OwOP ~ C = O_~ D Og ~x oP; =0u oP; W2 =~w oe

(16)

where all of the above partial derivatives are evaluated at the chosen operating point and will, therefore, change depending on the operating point. This results in the well known linear "state-space" model format: dx dt - Ax + Bu + W~w y = Cx + Du + W2w (17) It is important to remember that the variables (x, u, w, y) are now deviation variables. That is they are deviations from the point the model is linearised about, not their original absolute values. Appendix B--Eigenvalues of a matrix The eigenvalues (2) of a matrix, A, are the solution to the equation (A - 2/)X = 0 where 1 is an identity matrix of the same dimension as A and X is a matrix of vectors called eigenvectors. For this equality to be true the determinant of the term in brackets must be zero. This represents a simple technique for calculating the eigenvalues by hand. Other numerical methods, not discussed here, are used in practice to calculate the eigenvalues of larger matrices. Appendix C--Time constants The time constants of a process can be used to quantify the dynamics of the various system variables. The concept of a time constant is best described with an example. Here we will look at the differential equation for soluble inert COD (S0 in a completely mixed activated sludge tank:

M . A . Steffens et al.

606 dS~ d---t = F'SI,I, -- F.St

(18)

where S~ is the soluble inert COD concentration in the tank, SLy, is the influent soluble inert COD concentration and F is the water flowrate (assume in = out). The analytical solution to equation (18) assuming the initial S~ concentration is zero is: St = SL~.(I - e'~) (19) w h e r e z i s t h e t i m e c o n s t a n t ( z = V/F). For a simple, first-order system such as this the time constant is equivalent to the time it takes for a system to reach 63% of the steady state value after a step in the input variable (Su.). The time constant, therefore, provides the modeller with an estimate of the response time of state variables to changes in input variables (manipulated or

disturbance), In a linear, first-order system the negative, inverse of the

series:

x = I + At + A 2~ + HOT x0

where HOT presents the higher order terms. Suppose the objective is to check whether the state x~ is slow mode reducible. We partition the vector of states and the A matrix as follows: - FA"

AI21x = [X'] x2

I

(25)

A - [_A2, A~J

where x2 is a vector containing the remaining state variables. Substitute the partitioned matrix and vector into equation (24) to get:

x,x2 = I+[AA;:

A2~_l

LA2' A'2_] =

eigenvalues, 2, represents the system time constants, z (equation (20)). T= -~

(24)

+ HOT)(Xl_'.: ")/\~/

(26)

(20)

Even though the system of interest possess neither of these characteristics (i.e. coupled, higher order process), the eigenvalues still provide a good indication of the "time constants",

D--Derivation of the slow mode reduction parameter(9) To quantify the accuracy of a slow mode reduction, the difference between the reference model and reduced model solutions is approximated. We start with the linear model:

Equation (26) represents an approximate solution to the reference model. If we set A. and AI2 to zero we get the solution to the reduced order model for all states. By subtracting the two solutions we arrive at the errors incurred by reducing the state of interest (x0, for x~ and the remaining states, x2.

Appendix

~=

Ax + Bu

(21)

[e,] ([A ez = 0

0AI2] [A~,+AmA2, A.Am+A,2A221t2 t+ A21Au AzlAI2 ] -i + HOT'J(x:I:~ /v" /

(27)

where u is a vector including all inputs (manipulated and disturbance). The full model solution is calculated using the following expression (Ross, 1974):

The slow mode reduction parameter, ~P, is then evaluated for each state i, as follows: ~i = max(y~, y2) where

x = e~'xo + e A' j'tf eA'Budz

~'~ = e~

Xl

(22)

to

where x0 is a vector of the initial values of the states. This solution is made up of a complementary and particular solution. The assumption that much the dynamic response is made up of the complementary solution is made here. Equation (22) then reduces to:

of

x = e'~'Xo

(23)

The matrix exponential term is then expanded as a power

and (28) ~,z = max {.~t) for j = 1, 2 . . . . . n - 1 \x2j] That is, the parameter is the maximum of the complete vector of relative errors (including the state of interest and the remaining states).