- Email: [email protected]
Simulation of a nonlinear distributed parameter bioreactor by FEM approach
SIMULATION
PRACTICE
ELSEVIER
Simulation
U’THEORY
Practice and Theory 5 ( 1997) 199-216
Simulation of a nonlinear distributed parameter bioreactor by FEM approach M.T. Nihtila a, J. Tervo a,*, J.P. Kaipio b a Department of Computer Science and Applied Mathematics, University of Kuopio. t?O. Box 1627, FIN-70211 Kuopio, Finland b Department of Applied Physics, University of Kuopio, PO. Box 1627, FIN-7021 1 Kuopio, Finland Received
1 December
1994;
revised 19 March 1996
Abstract A two-dimensional model of a distributed parameter bioreactor is studied. The model describes biomass growth and substrate consumption in a continuous-flow fixed-bed denitrification process aimed to lower the level of harmful nutrients (substrates) in communal waste water. The independent variables of the model are the evolving time and the distance measured from the beginning of the reactor tube. The model is being applied to control the substrate level at the output of the reactor. The incoming flow is the control variable, the incoming substrate a disturbance. The existence of solutions of the corresponding system of partial differential equations is considered. A finite element method is applied in discretising the single space variable to obtain a system of ordinary differential equations. The simulation results that complete the paper are in accordance with the physical insight of the real process. @ 1997 Elsevier Science B.V. Keywords: Bioreactors;
Distributed-parameter
systems;
Simulation;
Finite element method
1. Introduction Control theory for distributed parameter systems has been developed for several classes of partial differential equation (PDE) models including boundary, distributed and moving control actions, see e.g. [5,6,15,20]. Processes described by nonlinear PDEs are natural in many chemical and biological systems. Numerical simulation and synthesis of feedback control algorithms for systems of nonlinear PDEs require most often analyses * Corresponding
author. E-mail: [email protected]
0928-4869/97/$17.00 @ 1997 Elsevier Science B.V. All rights reserved PIISO928-4869(96)00008-O
M.T Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
200
based on individual process and nonlinearity structures. Specific approaches are also needed to prove existence and uniqueness of the solutions of the corresponding initial and/or boundary value problems. Several methods for simulation of initial/boundary value problems of systems of PDEs have been developed and applied. Typically, one class of methods is based on approximate separation of space and time variables. A finite separation series like in the collocation methods, see e.g. [ 29,25,9], approximates then the true solution. Another class of methods utilizes the weak formulation of the problem in the corresponding Sobolev spaces. Galerkin’s finite element method (FEM), see [ 14,27,28,19], that belongs to the latter class, is the method to be used in this paper. Actually, all these so-called semidiscretization methods discretise the spatial domain and give systems of nonlinear ordinary differential equations (ODE) as approximations to the original distributed problem. Application of a suitable semidiscretization method is crucial for the controller design in pointwise controlled distributed systems. A typical control objective for a system of nonlinear ODES is that in closed loop control the overall system is linear, i.e. the idea of linearizing control is applied, cf. [ 131 and the references therein. The validity of the Kalman criterion is in many cases a necessary criterion so that the static (and dynamic) state feedback linearization method applies, see [ 21. The so-called minimum-phase property that roughly speaking means that the transfer function of the linearized system around the steady state has all zeros with negative real parts, is also an essential requirement in the design of the closed loop systems [ 131, especially when problems like asymptotic output tracking are studied. The typical control laws in the problems of asymptotic output tracking contain one or more external reference inputs that must satisfy the Hurwitz condition [ 13,161. Input-output behavior of the open loop system may in some cases exhibit a specific feature which indicates that the associated open loop system may be (nearly) of nonminimum phase type or equally the closed loop system may be (nearly) internally unstable. This feature can be described in the following way. A step change in the control variable drives the output (controlled variable) at the beginning of the step response to the direction opposite to the final steady state value as compared to the initial steady state (cf. Fig. 5 of the simulation results). Carefully planned and performed simulations of the original distributed system can reveal this feature as can be seen in this paper. We have found that this feature depends on the number of nodes used in the FEM scheme. The (as yet unpublished) closed loop simulations seem (from the stability point of view) to support this assertion. In the following we describe the model studied in this paper and present an overview of the control theoretic approaches concerning this model. The model consists of a pair of nonlinear partial differential equations au1 -
at
= (P(W,U2)
ih2 7
=
-
kf)w,
iA42 D-man
au2
lx3
-
hrLL(W>U2)W
for x E 10, 1 [, t > 0 with the following
boundary
and initial conditions
M.T. Nihtiki et al. /Simulation
Practice
L
and Theory 5 (1997) 199-216
201
x=0
I.
x=1
Fig. 1. A schematic diagram of the bioreactor.
$(O, t)
=
9
(U2(O*t)
$(l.r)
-S(f)),
=o
for t > 0 and
Ul(X,O) =U1(x),
u2(x,O)
= U2(x)
(3)
for x E IO, 1 [. Here ~1 = ut (x, t) and u2 = u~(x, f) denote
biomass
and the substrate,
respectively.
the concentrations of the The specific growth rate (due to Contois [ 3 1)
u2 Pu(Ul1 u2) = PI?1 k2ut + ~2
makes the system ( 1) nonlinear. The input flow c (control variable) and the input substrate concentration S, (disturbance variable) are generally functions of t (distributed parameters). In control theoretical setting the variable to be controlled is u2( 1, t), see Fig. 1. The other parameters are positive constants, The initial condition functions Ui , U2 are typically chosen as the steady state functions before the simulated step changes the input. In the steady state the parameters c and S, are independent of time t. The values of c and S, which corresponds to the steady state solutions iii and u:! are denoted by C and 3,. For the model study and for the discussion of the parameters of the system (l)-(3) we refer to [8,9,25]. The existence of solutions for problems such as ( l)-( 3) has been extensively studied. For basic monographs relating to general semigroup theory of evolution equations we refer to [ 11,26,23,1]. In the monographs by Smoller [ 241 and Leung [ 181 one has considered from a quite different point of view the existence, uniqueness, stability results
M.T. Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216
202
and their applications for the equations related to the problem ( I)-( 3). In Section 2 we establish an existence result of solutions for ( 1) -( 3) that is proved as an application of the semigroup theory in [ 121. In practice the key technique for studying the behaviour of solutions of the dynamical system (l)-(3) is to approximate the system (l)-(3) by a nonlinear finite dimensional ordinary differential equation system. A method of orthogonal collocation, see [ 25,9], has been applied to construct the approximations. In Section 3 we deal with the Galerkin finite element method (FEM) to compute numerical solutions for the problem ( 1) -( 3). We construct a FEM scheme to evaluate the solutions. There are some convergence results (error estimates) that guarantee that the semidiscrete Galerkin approximation uh tends to the exact solution u = (~1, ~42) when the mesh size h tends to zero. In our case the convergence of the FEM scheme (in the case where c and S, are constant) can be obtained as an application of the theory for semidiscrete Galerkin type approximations to semilinear evolution equations (see [ 141 and the references therein). The related convergence results for a (single) parabolic equation can be found, for example, in [ 27,28,4,7]. In Section 4 we consider an associated linearized model. The behavior of the solutions of the corresponding semidiscrete model as t + co is studied. Finally, in Section 5 we express some numerical results. We compute the output y(t) = ~2 ( 1, t) as a function of time t for different choices of c and S,, see Fig. 2. The simulation results are consistent with the qualitative nature of the physical process. As indicated in Fig. 3 the qualitative (and quantitative) behavior of the numerical solution for the linearized model is similar to that of the original nonlinear model. This is a confirmation of the correctness of the numerical experiments. The numerical results suggest that also the original (partial differential equation) system ( 1) -( 3) is asymptotically stable, see Fig. 4. In Fig. 5 we see the above mentioned “fingerprint” of the potential vicinity of the non-minimum phase property of the system.
2. Existence
of solutions
2.1. The formulation
of the problem
Let G be an open set in lR and let A be any interval in R. By C’(A), 1 E Nc we denote the space of all I times continuously differentiable functions f : A H W. C’( A, X) is the space of all I times continuously differentiable functions f : A H X, where X is a Banach space. In the above spaces we use the standard topologies. The space L2(G) consists of all Lebesgue square integrable functions f : G H C. Lz(G) is equipped with the usual inner product
(f, do
= /-
fWg(x)dx.
G
Denote by H’(G),
2 E NO the Sobolev space
M.7: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216 djW dxi E
WEL2(G)
equipped
203
exists for 0 < j < 1 in the weak sense
L2(G)
>
with the inner product
(w,@), $$,~)).
Consider
the initial/boundary
au1 = (,4uI,u2) at A.42 -_=D_at
-
]
-
h(/-4W,~?))~l
0, T[ with the following
au7__ c(t)
ax
kd)UI,
c(t)%
ax2 x
of the form
(4)
d2Ll2
for (x, t) E G
value problem
- 3U2(OJ)
Z(1.t)
-&I(t)), u2(x,O)
UI(X,O) =ul(x>,
boundary
and initial conditions =0
x E
=I12(x),
fort
E lO,T[,
G.
(5) (6)
Here T > 0 is the time up to which the evolution of the process is observed and G is the interval IO, 1 [ c IR. The distributed parameters c and S, are assumed to bc positive functions in C2 ( [ 0,oc [ ) such that c(0)
= z,
S,(O) = 3,.
The initial functions (/@I,&) Dz;
(7)
El ,u;! are obtained
with the boundary
(8)
I*?(X) =
=0,
x’E G,
conditions
= 3is?(O)
The problem
equations
- kd)& =O,
_ CU’ 2 - ~,(,u(EI,E~))UI
z;(o)
from the steady-state
z;-:( 1) =o.
- S,),
(9)
(8) -( 9) has the closed form solution qcosh(q(
1 - x)) + p sinh(q(
qcoshq Pm -
kd-
El (xl = ------u2(x) k2kd
+ (a + p) sinhq
1 - x)) S,
>
exp(px),
9
where
c P=z'
a=
kl (b&n- kd) k2C
’
q=&+cc
2 )
kl(j-h - kd) ‘I2 kD 2
).
204
M.I: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
2.2. Existence
of solutions
In this section we consider the existence of a dynamical solution for the problem (4)-(6). The existence result guarantees that the model is mathematically relevant and that the simulations make sense. Existence analysis for equations like (4)-( 6) is often based on the application of semigroup theory. From the mathematical point of view the application of semigroup theory to (4)-(6) is interesting since the domain of the infinitesimal generator, say V( A( t) ) depends on the time t. In such a case there are roughly speaking two possibilities: one can either try to apply the “parabolic case” or the “hyperbolic case” [ 171. One of the the characteristic features of the latter approach is that D(A( t)) is or can be made independent of t with some transformations, This method can be said to have been applied e.g. in [ 121 although in that report the infinitesimal generator generated an analytical semigroup. For the “parabolic case” we refer to [ 11. We express shortly the abstract formulation of (4)-( 6) for later needs in Section 3. Define a linear operator L(t) : L2( G) H L2 (G) such that D(L(t)) L(t)w
{w H*(G) 1(g - F-,(O) =o, S(l) co}, = -(Lg - c(t$). =
E
(10)
Note that the expressions (dw/dx - D-‘c( t) w) (0) and (dw/dx) ( 1) make sense since the Sobolev space H2( G) can be continuously imbedded into the space C’ (G) [ 23, Theorem 7.1.21. Let A(t) : L2 (G) 2 ++ L2 (G) 2 be a linear operator defined by D(A(t)) A(t)u=
= L2(G)
x D(Ut)),
for u = (ut,u2)
(kdul,L(t)uz)
Define the mappings
F and P by
F(t,o)=(~u(u1,U2+Sa(t))u1,--kl~U(u1,~2+S,(t))ul ~~~~~~=~~~I~1I~I~2+~~~~~l~l~1I~--klrU~I~1I,I~2+Sa~t~l~lU1J The
problem
(4)-(
du dt + A(t)u
(11)
E D(A(t)).
(12)
-S:,(t)>, -s:,w.
(13)
6) may be put in the abstract form
= F(t,u), (14)
u(0) = uo := (El,i& where u = (ui,u2)
- S,(O)),
= (ut,u2
- So(t)).
Theorem 1. Suppose that c and S, are C2( [ 0, oo[ ) functions such that (7) holds. Then the problem (4)-( 6) has a unique classical global positive solution u = (~1, ~42) such that
M.T Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216 u1
nc([o,~[,c(G)),
Ec'(10,~:,[,C(m
u~EC’(IO,CO[,C(G)
du dt + A(t)u
(15) nc([0,03[,c(G)).
nc(10,m[,C2(G))
The proof follows by considering = I?t,v),
the associated
3.1. Prelimirtary
( 16)
problem
u(0) = ug
and applying the abstract semigroup theory, positivity equations. This proof is given in [ 121.
3. The semidiscrete
205
and regularity
results of parabolic
FEM-scheme
semidiscretization
In the sequel we give a semidiscrete model based on the Galerkin method. Error estimates for the obtained FEM scheme can be found in the literature (see e.g. the above-mentioned references [ 14,27,28,7] ). The existence of error estimates for the obtained scheme provides a stronger basis for the simulation considerations than many other methods. It could also be useful to compare the the results obtained with the FEM scheme and the results obtained by other numerical approaches (such as the collocation and difference methods) in this field of biotechnology. As discussed above the abstract form of our problem is z
+
A(t)u
= F(t,v), (17)
v(0) = (ii,,&
- S,(O)) A uo.
Suppose that u is the solution of ( 17) (guaranteed by Theorem 1) and that w = (WI, ~2) E LZ (G) x H’ (G). Then we find by integrating in parts (here [ . , .] denotes the L2 (G) 2 inner product) [A(t)u(t)>wl=
[(kdul(t),L(t)v2(t)),(w1,W2)1 Gkd(Ul(t),Wl)O
+
=kd(Ul(t),‘+‘I)O
-D
(L(t)U2(f),W2)0
(
=kd(L’,(t),w,)o+D(~(f),~)o
+C(t)U2(Ovt)W2(0)
+C(t~(~~f),W2)0
(18)
because (dvJdx) ( 1, t) = 0 by the definition of D( L( t) ) and (du2ld.x) (0, t) = D-‘c( t)vz(O, t). Since by the Sobolev Lemma H’(G) c C(G) and since the functions in H’(G) are absolutely continuous, the partial integration in (18) is allowable and w:!(0) makes sense.
206
M.7: Nihiilii
: (L2(G)
Let B(t) B(t>(v,w)
et al./Simulation
x H’(G))2
H R be the bilinear
Then we find that the variational
$w1+Btt)tutt),w)
[u(O), w] = [~a, w] Hence the solution
[
form
=B(t)((ul,U2),(Wl,W2))
+ C(f)U2(O)W2(0)
[
and Theory 5 (1997) 199-216
Practice
$3
(19)
w2)o.
form of the problem
(17) is
= [F(t,utt)),wl, for w E Lz(G)
(20)
x H’(G).
u = (~1, ~2) = (ur , v2 + S,(t) ) satisfies the variational
$w1+B(t)(utf),w)
[u(O),wl
+ c(r)(
=-Btt)t(O,-S,(t)),w)
,~,\
+ [F(u(t)),w],
\L’J
for w E Lz(G)
= [~O~Wl
equation
x H’(G),
where ua A (Ui , ii2 ) and F((UlvU2))
U2)Ul,_kl~U(UI,U2)Ul).
A (pu(w
(22)
Let Sh I [@i,...,@~] be a finite-dimensional subspace of L2 (G) Galerkin FEM-approximation Uh of the exact solution u is defined as
x H’ (G).
The
M
Uh(t) A
c
y(t)@j
(23)
j=l
for which
r$@k 1
=-B(t)(tO,-&tt)),@k)
+Btt)(Uhtt),@k)
[Uh(O),@k]
+ [F(uh(t)),@k], (24)
= [uo,@k]
for k = 1,. . . , M.
3.2. The semidiscretization corresponding to the piecewise linear element basis In the sequel we choose the piecewise
linear element
space (here M = 2N) (25)
Sh=[(~l,0),...,(~N,0),(0,~l),...,(O,~N)l,
where (X-Xj-l)/(xj-Xj-1)~ 4jj
=
(Xj+l
i 0, for 2 < j < N - 1 and
-
X)/(Xj+l
X E [xj-l,xjl, -
Xj>,
X E [XjtXj+ll,
otherwise
(26)
M.T Nihtiki et al. /Simulation (x2 $1 (x)
-x)/(x2
-x1),
x E
=
and Theory 5 (I 997) 199-216
[x1,x21,
(~-xN-l)/(xN--N-l),
XE
0,
otherwise,
=
207
(27)
otherwise,
0,
dN(X)
Practice
[XN-l,XN],
(28)
where {x] = 0,x2,. , x,+1, xN = 1) is a partition of the interval G = IO, I [. Supposing that x, -xj-] =h we have N= l/h+ 1. In the case of the element space (25) we find that 2N u,,
=
y(t)@j
c
Z!Z
j=l
j=l
j=l
where G Vj,
U]j
j=l ,...,N.
u2j A y+Nq
We denote N
N
ii,
k
ti]h
=
c
ii2 A U2,, = C
ulj(t)+j*
j=l
U2j(t)4,
j=l
so that i7] and ii2 are approximations of U] and ~2, respectively. Substituting uh = (C], ii2) in Eq. (24) we obtain for any @k A (@]k, @2k), k = 1,...,2N
This system can be expressed
-
c(t)
(ii2(0,
t)
-
+
((cL(hT~2)
+
(-kl~u(iil,ii2)iil,~2k)o.
-
&J(t)
)@2k(O)
kf)h.@lk)O
in matrix form as
MV’=-DKV-c(t)(RV+f(t,V))+F(V), where the prime denotes differentiation v=
V(t) = (u]]
f(r,
V) = (09 . . .
with respect to time and where
“’ U],,! @] ‘.’ U2N)*,
(vN+l
-So(t))
N F(V)
(29)
w*, N-l
= ((f1(@9a2),+1)0
“’
(fl(h,fi2),4N)O
(f2(h9~2),41)0
“’
(f2(h,i2).4N)O)
M.7: Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216
208
where M(j9 k) =
(djv
a
4k)O9
k) = (cq* &Jo,
k(j,
k) = (4j 9 +L>O.
In the above notations M, K and R are 2N x 2N-matrices and A?, R and a are N x Nmatrices whose elements in the jth row and kth column are fi( j, k), Z?( j, k) and k( j, k), respectively. The two functions fi and f2 in the vector F(V) are fl(h,ii2)
= (p(h,a2)
Direct computation
-
kd)iil,
f2(&,&)
= -hpU(h,~22)&.
yields
The integration of (fi( ii, 4)) &)o for 1 = 1,2 is performed approximations: (A) For 1 < k < N we use the Mean Value Theorem:
(B)
With k = 1 we use Simpson’s
(fi(69~2)2),c731)o
Formula:
= ~(fr~~l*~~~~l~o~~u2l~~~~l~o~~
+ 4fi(~ll(t)+l
(h/2)
+ Q2(0$2(h/2)*
by using the following
M.T. Nihtilii et ai./Simulation Practice and Theory 5 (1997) 199-216
and similarly
Finally,
with k = N
we consider
MV(O) =
the initial condition
((El,~l)O
..*
@1,4N)O
from which V(0) can be solved. For simplicity we approximate Then we have v(o)
(24).
Using the notations
(li29431)o
.”
MV’=-DKV-c(t)
(RV+f(t,V))
Y(f)
(31) for uz is
=Uzj(t),
4. Linearization
=
yields
j=l,...,N.
the approximation =U:!(l,t)
u2(Xj)+j).
+F(V),
whichatthenodesxi=(j-l)h,j=l,...,N+l
Especially
N
form
= (a,P).
fiz(xj,t)
(30)
-t&p)>
where aj =Ut(Xj) and pj =Zz(Xj). As a result we obtain the semidiscrete
The approximation
above we obtain
@2?4N)O)*>
ue = (tit, ii2) K5 (~,N=lEl(xj)4j3~j31
=(al,...,~N,Pl,...,PN)
V(0)
209
for the output y(t)
i&(1,t)
is
= !&N(t).
of the model
We outline some basic facts on a linearization corresponding to the system ( l)(3). The linearized model and the associated semidiscrete model can be used in the formulation and the proofs of various stability results. Let I and S, be as above the values of parameters c and S, corresponding the steady state solutions Ut and i&, respectively. Furthermore, let ut =ut
-iit,
u2 = u2 -ii2.
210
h4.7:Nihtilti et al. /Simulation Practice and Theory 5 (1997) 199-216
The linearization
around
the steady state E = (El, i&) leads to the linear system
(see
1221)
-
dt
= -al VI + a2U2,
au2 at
=
D a2u2 --C(f)% ax2
(32) -
(c(t)
-E)z
for x E IO, 1[, t > 0 with the boundary
$(O, t)
= +
$(l,t)
=o
U2(Ovt)
+kd)Ul
- klazU2
(33)
and initial conditions
-C)&(O)
$((c(t)
+
-kl(-al
+ESa),
-c(t)&(t)
(34)
for t > 0, &(x,0)
UI (XV 0) = 0,
=0
(35)
for x E IO, l[. The FEM scheme for (33)-(3.5) is obtained semidiscrete system is of the form (see [ 221)
as in Section
3. The corresponding
w’=E(t)W+H(t),
(36)
W(0) = 0,
(37)
where
W = (wll(t),
. . . . WIN(t),
wz(t),
. . . . wz,v(t))*
and CE,wlj(t)+j
of the U1 and U2, respectively. CFr w2j(t)+j are the approximations concerning the long time behavior of the solutions W of (36)-( 37).
and
We state a result
Theorem 2. Suppose that pu,, > kd holds and that c and S, are positive functions C’( [O,oo[) such that with some 0 < a < 1 /c(t) - c(t’) 16 qt I&(t)
- tqa,
c(m)
(38)
- Sa(t’>l < Clt - try
for t, t’ E [0, oo[. Furthermore, := tlicc(t)
and
in
(39)
assume that the finite positive
limits
&(cQ> :=/i%&(t)
(40)
exist. Let W : [O,oo[ H R 2N be the solution of the system (36) -( 37). Then there exists W(CXJ> E R2N such that IIW(t) where W(m)
- W(~)II
--t 0,
x(t) IIdW
does not depend on WO.
The proof is given in [ 221.
0
-+O,
ift+00,
(41)
M.7: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
211
Since the nonlinearity F(V) in the problem (31) is a C2( IO, co [ 2N) function, it follows that stability results such as Theorem 2 are also valid for the nonlinear system (3 1) We do not give here a detailed verification for this assertion but we remark that this result is apparent from our simulations in Fig. 4. The stability of the semidiscrete model suggests that also the original exact system (4)-( 6) has a similar kind of stability property (in appropriate spaces). The exact stability analysis for (4) -( 6) is complicated and we have not considered it.
5. Numerical
results
Eq. (3 1) is solved with the the standard Euler method with step size At = 0.1. The number of nodes is N = 8. The fixed parameters are chosen as follows D = 0.005,
kd = 0.05,
kl = k2 = 0.4,
/_Lu,= 0.35.
At first we solve the steady-state equations (8)-(9) ? = 0.1 and $, = 5 for the control parameter and for the input respectively. The corresponding solution is denoted At the moment t = 0 we change the control parameter cases: (I) c(t) =O.l +O.OS(l -e-IO’), (2) c(t) =O.l +O.l(l -e-lo’), (3) c(t) =O.l +O.15(1-e-1o’). (4) c(t) = 0.1 +0.2( 1 -e-lo’). The input substrate concentration is in every case chosen S,(t)=5
choosing the constant values substrate concentration at the by (Ul , U2). c = E to one of the following
to be
forO
So(t) =5 - 2( 1 - e-‘“(f-80)‘)
for 80 < t < 160.
The choices for c(t) and S,(t) were made to enable the evaluation of the system in cases of fast transitions that correspond to realistic situations. Step functions were avoided in order to maintain c(t) and S,(t) in C2 ( [ 0,0~) [ ) . We remark, however, that it is also allowable to use stepwise changes that correspond to the recovery of the system from incorrect initial values and with constant c(t) and S,(t). Infrequent stepwise changes in the parameters c and S, would correspond to a situation in which the error estimates given in [ 221 apply directly for adequately small changes. The above selections mean that at the moment r = 0 the control c is rapidly changed from 0.1 to the values 0.15, 0.2, 0.25 and 0.3 corresponding the cases 1, 2, 3 and 4. The parameter S, is rapidly changed at the time t = 80 from 5 to 3. The evolution is simulated to the time t = T = 160. Note that at the times t = 80 and t = 160 the system has nearly achieved a steady state, see Fig. 2. The solutions of the original and linearized systems are shown in Fig. 3. Here the control c(t) corresponds to case 1 but the step in S, (at time t = 80) is changed from 2 to 0.2.
M.T Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216
212
0
50
100
150
0
50
Time
100
15
I 0
Time
2
2
Case3
1.5 2
5 _a
21
a
Case4
1.5
1 0.5
0.5 i
#i I
0’
0
50
100
O0
150
Time
Fig. 2. The output y(t)
Time
= Q( 1, t) with the controls of the cases 1-4.
,
___________---
0.7. s _a a 0.6 -
0.5 -
0.4
f
0.31 0
Fig. 3. The output y(t)
20
40
60
Tvkfe
100
120
140
1 0
= uz( 1,t) (solid) and the output U2( 1, t) + ii2( 1) of the linearized system (dashed).
In Fig. 4 we see the evolutions when c(t) is chosen according to case 1 with S, = 5 and when the initial function is ua = Gh,E2), uo = G,E2)/2 or uo = GG,i72)/4, respectively. Fig. 4 suggests that the system is asymptotically stable. The “fingerprint” of the potential vicinity of instability of the system for feedback-linearizable control design (cf. [ 161) is clearly visible in Fig. 5 where the output y(t) corresponds to
M.7: NihtilSi et al. /Simulation
1.67
1
213
Practice and Theory 5 (1997) 199-216
I
j
5 go.8 0
t
I
01 0
20
40
60
80
100
120
140
160
Time
Fig. 4. The effect of different initial values on the output ~2 ( 1, I). Solid line corresponds and dotted to Z/4.
0.32: 0
1
2
3
4
5
6
to ii. dashed to ii/2
7
Time
Fig. 5. A demonstration of the potential vicinity of the system with .S, E 5 and D = 0.001. c(t) = 0.07 + 0.05( 1 - exp(-t))
instability
corresponding
to control
214
M.7: Nihtilii et al. /Simulation
Practice and Theory 5 (I 997) 199-216
Table 1 Exact and numerical steady state solutions second transitions) in cases l-4
1
Case Exact Galerkin
0.7588 0.7757
1.1742 1.1890
I.541 1 1.5599
steady state solutions 1
when
4
1.8688 1.8799
t = 160 in cases I-4
2
3
4
0.7045 0.7085
0.9283 0.9302
1.1213 1.1219
control c(t) = 0.07 + 0.05( 1 - exp( -t)) The exact solutions to the steady state compared with the asymptotic solutions of exact and numerical steady state solutions
6. Concluding
3
Case
0.4553 0.4619
t = 80 (just before the
2
Table 2 Exact and numerical
Exact Galerkin
when
with S, 5 5 and D = 0.001. of the system calculated in Section 2.1 are the numerical scheme in Tables 1 and 2. The are easily seen to be consistent.
remarks
A widely accepted issue in modelling is that any model constructed for describing some true phenomenon can be valid only in limited situations. Consequently, a mathematical model is an idealization of the real world phenomenon. Thus in the study of the phenomenon via a model it is of utmost importance that the mathematical prerequisites for the model are correct. This means that existence, uniqueness and solvability issues are on a solid basis. If the model - in our case a PDE system - needs to be approximated further due to simulation requirements, these approximations with respect to the more complicated model must also be justified. The Galerkin finite element method is an ideal tool in this respect. Application of the Galerkin FEM in the problem of the equations (4) - (6) has been carried out to obtain results which can be compared with the results obtained by applying the orthogonal collocation method to the same problem (see e.g. [21] ). An advantage of the Gale&in FEM is that rigorous bounds can be found for the approximation errors. On the other hand some laborious pen-and-paper work was needed to obtain the approximating ODE system from the original PDE system. In the collocation method this semidiscretisation is somewhat more straightforward than in the Galerkin FEM. Proper understanding of the Galerkin FEM also needs knowledge of certain functional analytic results of Sobolev spaces. The value of the linearised PDE system studied lies in the fact that asymptotic stability of the semidiscrete system corresponding to the linearised PDE system also implies asymptotic stability of the semidiscrete system corresponding to the original nonlinear PDE system. These issues were studied via simulations in Section 5. The essential requirement in this implication is the attribute “asymptotic”. Bare stability is not implied in this way.
M.7: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
215
Further theoretical research includes convergence of the Galerkin FEM for the case of continuous and time-varying c and S,. This can most probably be shown by using the continuous dependence of the solution of the PDE system on these functions. Then the functions can be approximated by piecewise constant functions which produce a corresponding sequence of PDE solutions. Then, it is believed, this solution sequence converges to the solution corresponding to the original time-varying functions c and S,,. A possible non-minimum phase property of the approximating ODE system as a function of the number of mesh points adds to the complication of the controller design and requires further consideration.
References for semilinear parabolic evolution equations, Ann. Scuola Norm. Sup. Pisa CL Sci. 11 (4) (1984) 593-676. [ 2 ] B. Charlet, J. Levine and R. Marina, Sufficient conditions for dynamic state feedback linearization. SIAM J. Conk Optim. 29 (1991) 38-57. I3 I D. Contois, Kinetics of bacterial growth relationship between population density and specific growth rate of continuous cultures, J. Gen. Microhiol. 21 ( 1959) 40. 141 M. Crouzeix, V. Thomee and L.B. Wahlbin, Error estimates for spatially discrete approximations of semilinear parabolic equations with initial data of low regularity, Math. Comput. 53 ( 1989) 25-4 I. [ 5 I R.F. Curtain and A.J. Pritchard, Infinite Dimensional Linear Systems Theory, Lecture Notes in Control and Information Sciences, Vol. 8 (Springer, Berlin, 1978). 161 RF. Curtain and H. Zwart, Lecture Notes on Distributed Parameter Systems, Part 1, unpublished manuscript, 1992. [ 7 1 R. Dautray and J.L. Lions, Muthematical Analysis and Numerical Methods for Science and Technology, Vol. 2: Functional and Variational Methods, Vol. 5: Evolution Problems I, Vol. 6: Evolution Problems II, (Springer, Berlin, 1988 and 1989). I 8 1 D. Dochain, J.P. Babary and N. Tali-Maamar, Modelling and adaptive control of nonlinear distributed parameter bioreactors via orthogonal collocation, Aufomafica 28 (5) ( 1992) 273-283. 19 1 D. Dochain, Contribution to the analysis and control of distributed parameter systems with application to (bio)chemical processes and robotics, Thesis, UniversitC Catholique de Louvain, Belgium, 1994. I IO] A. Friedman, Parfial Di$erential Equafions (Krieger, New York, 1976). I I I I J.A. Goldstein. Semigroups of Linear Operators and Applications (Oxford University Press, Oxford, 1985). I 12 1 J.A. Goldstein and J. Tervo, Existence of solutions for a system of nonlinear partial differential equations related to bioreactors, Dynamic Syst. Appl., to appear. [ 131 J.W. Grizzle, M.D. DiBenedetto and F. Lamnabhi-Lagarrigue. Necessary conditions for asymptotic tracking in nonlinear systems, IEEE Trans. Automat. Contr. 39 ( 1994) 1782-1794. I 141 H.P. Helfrich, Error estimates for semidiscrete Galerkin type approximations to semilinear evolution equations with nonsmooth initial data, Numer. Math. 51 (1987) 559-569. I IS I K.S. Hong and J. Bentsman, Direct adaptive control of parabolic systems: Algorithm synthesis and convergence and stability analysis, IEEE Trans. Automat. Contr. 39 ( 1994) 2018-2033. I 16 I A. Isidori, Nonlinear Control Systems (Springer, Berlin, 1989), I 17 I T. Kate, Nonlinear evolution equations in Banach spaces, Proc. Symp. Appl. Math. 17 ( 1964) 50-67. ( 18 I A.W. Lang, Systems of Nonlinear Partial D@erential Equations: Applications to Biology and Engineering (Kluwer, Dordrecht, 1989).
[ I I H. Amann, Existence and regularity
I 19 1 Y. Lin and T. Zhang, Finite element methods for nonlinear Sobolev equations with nonlinear boundary conditions, J. Math. Anal. Appl. 165 (1991) 180-191. [ 20 I J.L. Lions, Optimal Control of Sysfems Described by Partial Di$erenfial Equations (Springer, 1971).
Berlin,
216
M.7: Nihtila et al./Simulation Practice and Theory 5 (1997) 199-216
[ 211 M.T. Nihtill, N. Tali-Maamar and J.P. Babary, Controller design for a nonlinear distributed parameter system via a collocation approximation, in: M. Flies% ed., Proceedings of the IFAC Symposium on Nonlinear Control System Design, NOLCOS’92, Bordeaux, France (1992) 447-452. [22] M.T. Nihtila, J. Tervo and J.P. Kaipio, Simulation of a nonlinear distributed parameter bioreactor by FEM approach, Dept. of Computer Science, Applied Mathematics Report Series A/ 1994/5, University of Kuopio, 1994. [ 231 A. Pazy, Semigroups of Linear Operators and Applications to Partial Dt’erential Equations (Springer, Berlin, 1983). [24] J. Smoller, Shock Waves and Reaction-Dtpsion Equations (Springer, Berlin, 1983). [25] N. Tali-Maamar, T. Damak, J.P Babary and M.T. Nihtill, Application of a collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simul. 35 (1993) 303-319. [26] H. Tanabe, Equations of Evolution (Pitman, London, 1979). [27] V. Thomee and L. Wahlbin, On Galerkin methods in semilinear parabolic problems, SIAM J. Numer. Anal. 12 ( 1975) 378-389. [ 281 V. Thomee, Galerkin Finite Element Methods for Parabolic Problems, Lecture Notes in Mathematics, Vol. 1054 (Springer, Berlin, 1984). [29] J.V. Villadsen and Michelsen, Solution of Dtfirential Equation Models by Polynomial Approximation (Prentice-Hall, Englewood Cliffs, NJ, 1978).
PRACTICE
ELSEVIER
Simulation
U’THEORY
Practice and Theory 5 ( 1997) 199-216
Simulation of a nonlinear distributed parameter bioreactor by FEM approach M.T. Nihtila a, J. Tervo a,*, J.P. Kaipio b a Department of Computer Science and Applied Mathematics, University of Kuopio. t?O. Box 1627, FIN-70211 Kuopio, Finland b Department of Applied Physics, University of Kuopio, PO. Box 1627, FIN-7021 1 Kuopio, Finland Received
1 December
1994;
revised 19 March 1996
Abstract A two-dimensional model of a distributed parameter bioreactor is studied. The model describes biomass growth and substrate consumption in a continuous-flow fixed-bed denitrification process aimed to lower the level of harmful nutrients (substrates) in communal waste water. The independent variables of the model are the evolving time and the distance measured from the beginning of the reactor tube. The model is being applied to control the substrate level at the output of the reactor. The incoming flow is the control variable, the incoming substrate a disturbance. The existence of solutions of the corresponding system of partial differential equations is considered. A finite element method is applied in discretising the single space variable to obtain a system of ordinary differential equations. The simulation results that complete the paper are in accordance with the physical insight of the real process. @ 1997 Elsevier Science B.V. Keywords: Bioreactors;
Distributed-parameter
systems;
Simulation;
Finite element method
1. Introduction Control theory for distributed parameter systems has been developed for several classes of partial differential equation (PDE) models including boundary, distributed and moving control actions, see e.g. [5,6,15,20]. Processes described by nonlinear PDEs are natural in many chemical and biological systems. Numerical simulation and synthesis of feedback control algorithms for systems of nonlinear PDEs require most often analyses * Corresponding
author. E-mail: [email protected]
0928-4869/97/$17.00 @ 1997 Elsevier Science B.V. All rights reserved PIISO928-4869(96)00008-O
M.T Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
200
based on individual process and nonlinearity structures. Specific approaches are also needed to prove existence and uniqueness of the solutions of the corresponding initial and/or boundary value problems. Several methods for simulation of initial/boundary value problems of systems of PDEs have been developed and applied. Typically, one class of methods is based on approximate separation of space and time variables. A finite separation series like in the collocation methods, see e.g. [ 29,25,9], approximates then the true solution. Another class of methods utilizes the weak formulation of the problem in the corresponding Sobolev spaces. Galerkin’s finite element method (FEM), see [ 14,27,28,19], that belongs to the latter class, is the method to be used in this paper. Actually, all these so-called semidiscretization methods discretise the spatial domain and give systems of nonlinear ordinary differential equations (ODE) as approximations to the original distributed problem. Application of a suitable semidiscretization method is crucial for the controller design in pointwise controlled distributed systems. A typical control objective for a system of nonlinear ODES is that in closed loop control the overall system is linear, i.e. the idea of linearizing control is applied, cf. [ 131 and the references therein. The validity of the Kalman criterion is in many cases a necessary criterion so that the static (and dynamic) state feedback linearization method applies, see [ 21. The so-called minimum-phase property that roughly speaking means that the transfer function of the linearized system around the steady state has all zeros with negative real parts, is also an essential requirement in the design of the closed loop systems [ 131, especially when problems like asymptotic output tracking are studied. The typical control laws in the problems of asymptotic output tracking contain one or more external reference inputs that must satisfy the Hurwitz condition [ 13,161. Input-output behavior of the open loop system may in some cases exhibit a specific feature which indicates that the associated open loop system may be (nearly) of nonminimum phase type or equally the closed loop system may be (nearly) internally unstable. This feature can be described in the following way. A step change in the control variable drives the output (controlled variable) at the beginning of the step response to the direction opposite to the final steady state value as compared to the initial steady state (cf. Fig. 5 of the simulation results). Carefully planned and performed simulations of the original distributed system can reveal this feature as can be seen in this paper. We have found that this feature depends on the number of nodes used in the FEM scheme. The (as yet unpublished) closed loop simulations seem (from the stability point of view) to support this assertion. In the following we describe the model studied in this paper and present an overview of the control theoretic approaches concerning this model. The model consists of a pair of nonlinear partial differential equations au1 -
at
= (P(W,U2)
ih2 7
=
-
kf)w,
iA42 D-man
au2
lx3
-
hrLL(W>U2)W
for x E 10, 1 [, t > 0 with the following
boundary
and initial conditions
M.T. Nihtiki et al. /Simulation
Practice
L
and Theory 5 (1997) 199-216
201
x=0
I.
x=1
Fig. 1. A schematic diagram of the bioreactor.
$(O, t)
=
9
(U2(O*t)
$(l.r)
-S(f)),
=o
for t > 0 and
Ul(X,O) =U1(x),
u2(x,O)
= U2(x)
(3)
for x E IO, 1 [. Here ~1 = ut (x, t) and u2 = u~(x, f) denote
biomass
and the substrate,
respectively.
the concentrations of the The specific growth rate (due to Contois [ 3 1)
u2 Pu(Ul1 u2) = PI?1 k2ut + ~2
makes the system ( 1) nonlinear. The input flow c (control variable) and the input substrate concentration S, (disturbance variable) are generally functions of t (distributed parameters). In control theoretical setting the variable to be controlled is u2( 1, t), see Fig. 1. The other parameters are positive constants, The initial condition functions Ui , U2 are typically chosen as the steady state functions before the simulated step changes the input. In the steady state the parameters c and S, are independent of time t. The values of c and S, which corresponds to the steady state solutions iii and u:! are denoted by C and 3,. For the model study and for the discussion of the parameters of the system (l)-(3) we refer to [8,9,25]. The existence of solutions for problems such as ( l)-( 3) has been extensively studied. For basic monographs relating to general semigroup theory of evolution equations we refer to [ 11,26,23,1]. In the monographs by Smoller [ 241 and Leung [ 181 one has considered from a quite different point of view the existence, uniqueness, stability results
M.T. Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216
202
and their applications for the equations related to the problem ( I)-( 3). In Section 2 we establish an existence result of solutions for ( 1) -( 3) that is proved as an application of the semigroup theory in [ 121. In practice the key technique for studying the behaviour of solutions of the dynamical system (l)-(3) is to approximate the system (l)-(3) by a nonlinear finite dimensional ordinary differential equation system. A method of orthogonal collocation, see [ 25,9], has been applied to construct the approximations. In Section 3 we deal with the Galerkin finite element method (FEM) to compute numerical solutions for the problem ( 1) -( 3). We construct a FEM scheme to evaluate the solutions. There are some convergence results (error estimates) that guarantee that the semidiscrete Galerkin approximation uh tends to the exact solution u = (~1, ~42) when the mesh size h tends to zero. In our case the convergence of the FEM scheme (in the case where c and S, are constant) can be obtained as an application of the theory for semidiscrete Galerkin type approximations to semilinear evolution equations (see [ 141 and the references therein). The related convergence results for a (single) parabolic equation can be found, for example, in [ 27,28,4,7]. In Section 4 we consider an associated linearized model. The behavior of the solutions of the corresponding semidiscrete model as t + co is studied. Finally, in Section 5 we express some numerical results. We compute the output y(t) = ~2 ( 1, t) as a function of time t for different choices of c and S,, see Fig. 2. The simulation results are consistent with the qualitative nature of the physical process. As indicated in Fig. 3 the qualitative (and quantitative) behavior of the numerical solution for the linearized model is similar to that of the original nonlinear model. This is a confirmation of the correctness of the numerical experiments. The numerical results suggest that also the original (partial differential equation) system ( 1) -( 3) is asymptotically stable, see Fig. 4. In Fig. 5 we see the above mentioned “fingerprint” of the potential vicinity of the non-minimum phase property of the system.
2. Existence
of solutions
2.1. The formulation
of the problem
Let G be an open set in lR and let A be any interval in R. By C’(A), 1 E Nc we denote the space of all I times continuously differentiable functions f : A H W. C’( A, X) is the space of all I times continuously differentiable functions f : A H X, where X is a Banach space. In the above spaces we use the standard topologies. The space L2(G) consists of all Lebesgue square integrable functions f : G H C. Lz(G) is equipped with the usual inner product
(f, do
= /-
fWg(x)dx.
G
Denote by H’(G),
2 E NO the Sobolev space
M.7: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216 djW dxi E
WEL2(G)
equipped
203
exists for 0 < j < 1 in the weak sense
L2(G)
>
with the inner product
(w,@), $$,~)).
Consider
the initial/boundary
au1 = (,4uI,u2) at A.42 -_=D_at
-
]
-
h(/-4W,~?))~l
0, T[ with the following
au7__ c(t)
ax
kd)UI,
c(t)%
ax2 x
of the form
(4)
d2Ll2
for (x, t) E G
value problem
- 3U2(OJ)
Z(1.t)
-&I(t)), u2(x,O)
UI(X,O) =ul(x>,
boundary
and initial conditions =0
x E
=I12(x),
fort
E lO,T[,
G.
(5) (6)
Here T > 0 is the time up to which the evolution of the process is observed and G is the interval IO, 1 [ c IR. The distributed parameters c and S, are assumed to bc positive functions in C2 ( [ 0,oc [ ) such that c(0)
= z,
S,(O) = 3,.
The initial functions (/@I,&) Dz;
(7)
El ,u;! are obtained
with the boundary
(8)
I*?(X) =
=0,
x’E G,
conditions
= 3is?(O)
The problem
equations
- kd)& =O,
_ CU’ 2 - ~,(,u(EI,E~))UI
z;(o)
from the steady-state
z;-:( 1) =o.
- S,),
(9)
(8) -( 9) has the closed form solution qcosh(q(
1 - x)) + p sinh(q(
qcoshq Pm -
kd-
El (xl = ------u2(x) k2kd
+ (a + p) sinhq
1 - x)) S,
>
exp(px),
9
where
c P=z'
a=
kl (b&n- kd) k2C
’
q=&+cc
2 )
kl(j-h - kd) ‘I2 kD 2
).
204
M.I: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
2.2. Existence
of solutions
In this section we consider the existence of a dynamical solution for the problem (4)-(6). The existence result guarantees that the model is mathematically relevant and that the simulations make sense. Existence analysis for equations like (4)-( 6) is often based on the application of semigroup theory. From the mathematical point of view the application of semigroup theory to (4)-(6) is interesting since the domain of the infinitesimal generator, say V( A( t) ) depends on the time t. In such a case there are roughly speaking two possibilities: one can either try to apply the “parabolic case” or the “hyperbolic case” [ 171. One of the the characteristic features of the latter approach is that D(A( t)) is or can be made independent of t with some transformations, This method can be said to have been applied e.g. in [ 121 although in that report the infinitesimal generator generated an analytical semigroup. For the “parabolic case” we refer to [ 11. We express shortly the abstract formulation of (4)-( 6) for later needs in Section 3. Define a linear operator L(t) : L2( G) H L2 (G) such that D(L(t)) L(t)w
{w H*(G) 1(g - F-,(O) =o, S(l) co}, = -(Lg - c(t$). =
E
(10)
Note that the expressions (dw/dx - D-‘c( t) w) (0) and (dw/dx) ( 1) make sense since the Sobolev space H2( G) can be continuously imbedded into the space C’ (G) [ 23, Theorem 7.1.21. Let A(t) : L2 (G) 2 ++ L2 (G) 2 be a linear operator defined by D(A(t)) A(t)u=
= L2(G)
x D(Ut)),
for u = (ut,u2)
(kdul,L(t)uz)
Define the mappings
F and P by
F(t,o)=(~u(u1,U2+Sa(t))u1,--kl~U(u1,~2+S,(t))ul ~~~~~~=~~~I~1I~I~2+~~~~~l~l~1I~--klrU~I~1I,I~2+Sa~t~l~lU1J The
problem
(4)-(
du dt + A(t)u
(11)
E D(A(t)).
(12)
-S:,(t)>, -s:,w.
(13)
6) may be put in the abstract form
= F(t,u), (14)
u(0) = uo := (El,i& where u = (ui,u2)
- S,(O)),
= (ut,u2
- So(t)).
Theorem 1. Suppose that c and S, are C2( [ 0, oo[ ) functions such that (7) holds. Then the problem (4)-( 6) has a unique classical global positive solution u = (~1, ~42) such that
M.T Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216 u1
nc([o,~[,c(G)),
Ec'(10,~:,[,C(m
u~EC’(IO,CO[,C(G)
du dt + A(t)u
(15) nc([0,03[,c(G)).
nc(10,m[,C2(G))
The proof follows by considering = I?t,v),
the associated
3.1. Prelimirtary
( 16)
problem
u(0) = ug
and applying the abstract semigroup theory, positivity equations. This proof is given in [ 121.
3. The semidiscrete
205
and regularity
results of parabolic
FEM-scheme
semidiscretization
In the sequel we give a semidiscrete model based on the Galerkin method. Error estimates for the obtained FEM scheme can be found in the literature (see e.g. the above-mentioned references [ 14,27,28,7] ). The existence of error estimates for the obtained scheme provides a stronger basis for the simulation considerations than many other methods. It could also be useful to compare the the results obtained with the FEM scheme and the results obtained by other numerical approaches (such as the collocation and difference methods) in this field of biotechnology. As discussed above the abstract form of our problem is z
+
A(t)u
= F(t,v), (17)
v(0) = (ii,,&
- S,(O)) A uo.
Suppose that u is the solution of ( 17) (guaranteed by Theorem 1) and that w = (WI, ~2) E LZ (G) x H’ (G). Then we find by integrating in parts (here [ . , .] denotes the L2 (G) 2 inner product) [A(t)u(t)>wl=
[(kdul(t),L(t)v2(t)),(w1,W2)1 Gkd(Ul(t),Wl)O
+
=kd(Ul(t),‘+‘I)O
-D
(L(t)U2(f),W2)0
(
=kd(L’,(t),w,)o+D(~(f),~)o
+C(t)U2(Ovt)W2(0)
+C(t~(~~f),W2)0
(18)
because (dvJdx) ( 1, t) = 0 by the definition of D( L( t) ) and (du2ld.x) (0, t) = D-‘c( t)vz(O, t). Since by the Sobolev Lemma H’(G) c C(G) and since the functions in H’(G) are absolutely continuous, the partial integration in (18) is allowable and w:!(0) makes sense.
206
M.7: Nihiilii
: (L2(G)
Let B(t) B(t>(v,w)
et al./Simulation
x H’(G))2
H R be the bilinear
Then we find that the variational
$w1+Btt)tutt),w)
[u(O), w] = [~a, w] Hence the solution
[
form
=B(t)((ul,U2),(Wl,W2))
+ C(f)U2(O)W2(0)
[
and Theory 5 (1997) 199-216
Practice
$3
(19)
w2)o.
form of the problem
(17) is
= [F(t,utt)),wl, for w E Lz(G)
(20)
x H’(G).
u = (~1, ~2) = (ur , v2 + S,(t) ) satisfies the variational
$w1+B(t)(utf),w)
[u(O),wl
+ c(r)(
=-Btt)t(O,-S,(t)),w)
,~,\
+ [F(u(t)),w],
\L’J
for w E Lz(G)
= [~O~Wl
equation
x H’(G),
where ua A (Ui , ii2 ) and F((UlvU2))
U2)Ul,_kl~U(UI,U2)Ul).
A (pu(w
(22)
Let Sh I [@i,...,@~] be a finite-dimensional subspace of L2 (G) Galerkin FEM-approximation Uh of the exact solution u is defined as
x H’ (G).
The
M
Uh(t) A
c
y(t)@j
(23)
j=l
for which
r$@k 1
=-B(t)(tO,-&tt)),@k)
+Btt)(Uhtt),@k)
[Uh(O),@k]
+ [F(uh(t)),@k], (24)
= [uo,@k]
for k = 1,. . . , M.
3.2. The semidiscretization corresponding to the piecewise linear element basis In the sequel we choose the piecewise
linear element
space (here M = 2N) (25)
Sh=[(~l,0),...,(~N,0),(0,~l),...,(O,~N)l,
where (X-Xj-l)/(xj-Xj-1)~ 4jj
=
(Xj+l
i 0, for 2 < j < N - 1 and
-
X)/(Xj+l
X E [xj-l,xjl, -
Xj>,
X E [XjtXj+ll,
otherwise
(26)
M.T Nihtiki et al. /Simulation (x2 $1 (x)
-x)/(x2
-x1),
x E
=
and Theory 5 (I 997) 199-216
[x1,x21,
(~-xN-l)/(xN--N-l),
XE
0,
otherwise,
=
207
(27)
otherwise,
0,
dN(X)
Practice
[XN-l,XN],
(28)
where {x] = 0,x2,. , x,+1, xN = 1) is a partition of the interval G = IO, I [. Supposing that x, -xj-] =h we have N= l/h+ 1. In the case of the element space (25) we find that 2N u,,
=
y(t)@j
c
Z!Z
j=l
j=l
j=l
where G Vj,
U]j
j=l ,...,N.
u2j A y+Nq
We denote N
N
ii,
k
ti]h
=
c
ii2 A U2,, = C
ulj(t)+j*
j=l
U2j(t)4,
j=l
so that i7] and ii2 are approximations of U] and ~2, respectively. Substituting uh = (C], ii2) in Eq. (24) we obtain for any @k A (@]k, @2k), k = 1,...,2N
This system can be expressed
-
c(t)
(ii2(0,
t)
-
+
((cL(hT~2)
+
(-kl~u(iil,ii2)iil,~2k)o.
-
&J(t)
)@2k(O)
kf)h.@lk)O
in matrix form as
MV’=-DKV-c(t)(RV+f(t,V))+F(V), where the prime denotes differentiation v=
V(t) = (u]]
f(r,
V) = (09 . . .
with respect to time and where
“’ U],,! @] ‘.’ U2N)*,
(vN+l
-So(t))
N F(V)
(29)
w*, N-l
= ((f1(@9a2),+1)0
“’
(fl(h,fi2),4N)O
(f2(h9~2),41)0
“’
(f2(h,i2).4N)O)
M.7: Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216
208
where M(j9 k) =
(djv
a
4k)O9
k) = (cq* &Jo,
k(j,
k) = (4j 9 +L>O.
In the above notations M, K and R are 2N x 2N-matrices and A?, R and a are N x Nmatrices whose elements in the jth row and kth column are fi( j, k), Z?( j, k) and k( j, k), respectively. The two functions fi and f2 in the vector F(V) are fl(h,ii2)
= (p(h,a2)
Direct computation
-
kd)iil,
f2(&,&)
= -hpU(h,~22)&.
yields
The integration of (fi( ii, 4)) &)o for 1 = 1,2 is performed approximations: (A) For 1 < k < N we use the Mean Value Theorem:
(B)
With k = 1 we use Simpson’s
(fi(69~2)2),c731)o
Formula:
= ~(fr~~l*~~~~l~o~~u2l~~~~l~o~~
+ 4fi(~ll(t)+l
(h/2)
+ Q2(0$2(h/2)*
by using the following
M.T. Nihtilii et ai./Simulation Practice and Theory 5 (1997) 199-216
and similarly
Finally,
with k = N
we consider
MV(O) =
the initial condition
((El,~l)O
..*
@1,4N)O
from which V(0) can be solved. For simplicity we approximate Then we have v(o)
(24).
Using the notations
(li29431)o
.”
MV’=-DKV-c(t)
(RV+f(t,V))
Y(f)
(31) for uz is
=Uzj(t),
4. Linearization
=
yields
j=l,...,N.
the approximation =U:!(l,t)
u2(Xj)+j).
+F(V),
whichatthenodesxi=(j-l)h,j=l,...,N+l
Especially
N
form
= (a,P).
fiz(xj,t)
(30)
-t&p)>
where aj =Ut(Xj) and pj =Zz(Xj). As a result we obtain the semidiscrete
The approximation
above we obtain
@2?4N)O)*>
ue = (tit, ii2) K5 (~,N=lEl(xj)4j3~j31
=(al,...,~N,Pl,...,PN)
V(0)
209
for the output y(t)
i&(1,t)
is
= !&N(t).
of the model
We outline some basic facts on a linearization corresponding to the system ( l)(3). The linearized model and the associated semidiscrete model can be used in the formulation and the proofs of various stability results. Let I and S, be as above the values of parameters c and S, corresponding the steady state solutions Ut and i&, respectively. Furthermore, let ut =ut
-iit,
u2 = u2 -ii2.
210
h4.7:Nihtilti et al. /Simulation Practice and Theory 5 (1997) 199-216
The linearization
around
the steady state E = (El, i&) leads to the linear system
(see
1221)
-
dt
= -al VI + a2U2,
au2 at
=
D a2u2 --C(f)% ax2
(32) -
(c(t)
-E)z
for x E IO, 1[, t > 0 with the boundary
$(O, t)
= +
$(l,t)
=o
U2(Ovt)
+kd)Ul
- klazU2
(33)
and initial conditions
-C)&(O)
$((c(t)
+
-kl(-al
+ESa),
-c(t)&(t)
(34)
for t > 0, &(x,0)
UI (XV 0) = 0,
=0
(35)
for x E IO, l[. The FEM scheme for (33)-(3.5) is obtained semidiscrete system is of the form (see [ 221)
as in Section
3. The corresponding
w’=E(t)W+H(t),
(36)
W(0) = 0,
(37)
where
W = (wll(t),
. . . . WIN(t),
wz(t),
. . . . wz,v(t))*
and CE,wlj(t)+j
of the U1 and U2, respectively. CFr w2j(t)+j are the approximations concerning the long time behavior of the solutions W of (36)-( 37).
and
We state a result
Theorem 2. Suppose that pu,, > kd holds and that c and S, are positive functions C’( [O,oo[) such that with some 0 < a < 1 /c(t) - c(t’) 16 qt I&(t)
- tqa,
c(m)
(38)
- Sa(t’>l < Clt - try
for t, t’ E [0, oo[. Furthermore, := tlicc(t)
and
in
(39)
assume that the finite positive
limits
&(cQ> :=/i%&(t)
(40)
exist. Let W : [O,oo[ H R 2N be the solution of the system (36) -( 37). Then there exists W(CXJ> E R2N such that IIW(t) where W(m)
- W(~)II
--t 0,
x(t) IIdW
does not depend on WO.
The proof is given in [ 221.
0
-+O,
ift+00,
(41)
M.7: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
211
Since the nonlinearity F(V) in the problem (31) is a C2( IO, co [ 2N) function, it follows that stability results such as Theorem 2 are also valid for the nonlinear system (3 1) We do not give here a detailed verification for this assertion but we remark that this result is apparent from our simulations in Fig. 4. The stability of the semidiscrete model suggests that also the original exact system (4)-( 6) has a similar kind of stability property (in appropriate spaces). The exact stability analysis for (4) -( 6) is complicated and we have not considered it.
5. Numerical
results
Eq. (3 1) is solved with the the standard Euler method with step size At = 0.1. The number of nodes is N = 8. The fixed parameters are chosen as follows D = 0.005,
kd = 0.05,
kl = k2 = 0.4,
/_Lu,= 0.35.
At first we solve the steady-state equations (8)-(9) ? = 0.1 and $, = 5 for the control parameter and for the input respectively. The corresponding solution is denoted At the moment t = 0 we change the control parameter cases: (I) c(t) =O.l +O.OS(l -e-IO’), (2) c(t) =O.l +O.l(l -e-lo’), (3) c(t) =O.l +O.15(1-e-1o’). (4) c(t) = 0.1 +0.2( 1 -e-lo’). The input substrate concentration is in every case chosen S,(t)=5
choosing the constant values substrate concentration at the by (Ul , U2). c = E to one of the following
to be
forO
So(t) =5 - 2( 1 - e-‘“(f-80)‘)
for 80 < t < 160.
The choices for c(t) and S,(t) were made to enable the evaluation of the system in cases of fast transitions that correspond to realistic situations. Step functions were avoided in order to maintain c(t) and S,(t) in C2 ( [ 0,0~) [ ) . We remark, however, that it is also allowable to use stepwise changes that correspond to the recovery of the system from incorrect initial values and with constant c(t) and S,(t). Infrequent stepwise changes in the parameters c and S, would correspond to a situation in which the error estimates given in [ 221 apply directly for adequately small changes. The above selections mean that at the moment r = 0 the control c is rapidly changed from 0.1 to the values 0.15, 0.2, 0.25 and 0.3 corresponding the cases 1, 2, 3 and 4. The parameter S, is rapidly changed at the time t = 80 from 5 to 3. The evolution is simulated to the time t = T = 160. Note that at the times t = 80 and t = 160 the system has nearly achieved a steady state, see Fig. 2. The solutions of the original and linearized systems are shown in Fig. 3. Here the control c(t) corresponds to case 1 but the step in S, (at time t = 80) is changed from 2 to 0.2.
M.T Nihtilii et al./Simulation Practice and Theory 5 (1997) 199-216
212
0
50
100
150
0
50
Time
100
15
I 0
Time
2
2
Case3
1.5 2
5 _a
21
a
Case4
1.5
1 0.5
0.5 i
#i I
0’
0
50
100
O0
150
Time
Fig. 2. The output y(t)
Time
= Q( 1, t) with the controls of the cases 1-4.
,
___________---
0.7. s _a a 0.6 -
0.5 -
0.4
f
0.31 0
Fig. 3. The output y(t)
20
40
60
Tvkfe
100
120
140
1 0
= uz( 1,t) (solid) and the output U2( 1, t) + ii2( 1) of the linearized system (dashed).
In Fig. 4 we see the evolutions when c(t) is chosen according to case 1 with S, = 5 and when the initial function is ua = Gh,E2), uo = G,E2)/2 or uo = GG,i72)/4, respectively. Fig. 4 suggests that the system is asymptotically stable. The “fingerprint” of the potential vicinity of instability of the system for feedback-linearizable control design (cf. [ 161) is clearly visible in Fig. 5 where the output y(t) corresponds to
M.7: NihtilSi et al. /Simulation
1.67
1
213
Practice and Theory 5 (1997) 199-216
I
j
5 go.8 0
t
I
01 0
20
40
60
80
100
120
140
160
Time
Fig. 4. The effect of different initial values on the output ~2 ( 1, I). Solid line corresponds and dotted to Z/4.
0.32: 0
1
2
3
4
5
6
to ii. dashed to ii/2
7
Time
Fig. 5. A demonstration of the potential vicinity of the system with .S, E 5 and D = 0.001. c(t) = 0.07 + 0.05( 1 - exp(-t))
instability
corresponding
to control
214
M.7: Nihtilii et al. /Simulation
Practice and Theory 5 (I 997) 199-216
Table 1 Exact and numerical steady state solutions second transitions) in cases l-4
1
Case Exact Galerkin
0.7588 0.7757
1.1742 1.1890
I.541 1 1.5599
steady state solutions 1
when
4
1.8688 1.8799
t = 160 in cases I-4
2
3
4
0.7045 0.7085
0.9283 0.9302
1.1213 1.1219
control c(t) = 0.07 + 0.05( 1 - exp( -t)) The exact solutions to the steady state compared with the asymptotic solutions of exact and numerical steady state solutions
6. Concluding
3
Case
0.4553 0.4619
t = 80 (just before the
2
Table 2 Exact and numerical
Exact Galerkin
when
with S, 5 5 and D = 0.001. of the system calculated in Section 2.1 are the numerical scheme in Tables 1 and 2. The are easily seen to be consistent.
remarks
A widely accepted issue in modelling is that any model constructed for describing some true phenomenon can be valid only in limited situations. Consequently, a mathematical model is an idealization of the real world phenomenon. Thus in the study of the phenomenon via a model it is of utmost importance that the mathematical prerequisites for the model are correct. This means that existence, uniqueness and solvability issues are on a solid basis. If the model - in our case a PDE system - needs to be approximated further due to simulation requirements, these approximations with respect to the more complicated model must also be justified. The Galerkin finite element method is an ideal tool in this respect. Application of the Galerkin FEM in the problem of the equations (4) - (6) has been carried out to obtain results which can be compared with the results obtained by applying the orthogonal collocation method to the same problem (see e.g. [21] ). An advantage of the Gale&in FEM is that rigorous bounds can be found for the approximation errors. On the other hand some laborious pen-and-paper work was needed to obtain the approximating ODE system from the original PDE system. In the collocation method this semidiscretisation is somewhat more straightforward than in the Galerkin FEM. Proper understanding of the Galerkin FEM also needs knowledge of certain functional analytic results of Sobolev spaces. The value of the linearised PDE system studied lies in the fact that asymptotic stability of the semidiscrete system corresponding to the linearised PDE system also implies asymptotic stability of the semidiscrete system corresponding to the original nonlinear PDE system. These issues were studied via simulations in Section 5. The essential requirement in this implication is the attribute “asymptotic”. Bare stability is not implied in this way.
M.7: Nihtilii et al. /Simulation Practice and Theory 5 (1997) 199-216
215
Further theoretical research includes convergence of the Galerkin FEM for the case of continuous and time-varying c and S,. This can most probably be shown by using the continuous dependence of the solution of the PDE system on these functions. Then the functions can be approximated by piecewise constant functions which produce a corresponding sequence of PDE solutions. Then, it is believed, this solution sequence converges to the solution corresponding to the original time-varying functions c and S,,. A possible non-minimum phase property of the approximating ODE system as a function of the number of mesh points adds to the complication of the controller design and requires further consideration.
References for semilinear parabolic evolution equations, Ann. Scuola Norm. Sup. Pisa CL Sci. 11 (4) (1984) 593-676. [ 2 ] B. Charlet, J. Levine and R. Marina, Sufficient conditions for dynamic state feedback linearization. SIAM J. Conk Optim. 29 (1991) 38-57. I3 I D. Contois, Kinetics of bacterial growth relationship between population density and specific growth rate of continuous cultures, J. Gen. Microhiol. 21 ( 1959) 40. 141 M. Crouzeix, V. Thomee and L.B. Wahlbin, Error estimates for spatially discrete approximations of semilinear parabolic equations with initial data of low regularity, Math. Comput. 53 ( 1989) 25-4 I. [ 5 I R.F. Curtain and A.J. Pritchard, Infinite Dimensional Linear Systems Theory, Lecture Notes in Control and Information Sciences, Vol. 8 (Springer, Berlin, 1978). 161 RF. Curtain and H. Zwart, Lecture Notes on Distributed Parameter Systems, Part 1, unpublished manuscript, 1992. [ 7 1 R. Dautray and J.L. Lions, Muthematical Analysis and Numerical Methods for Science and Technology, Vol. 2: Functional and Variational Methods, Vol. 5: Evolution Problems I, Vol. 6: Evolution Problems II, (Springer, Berlin, 1988 and 1989). I 8 1 D. Dochain, J.P. Babary and N. Tali-Maamar, Modelling and adaptive control of nonlinear distributed parameter bioreactors via orthogonal collocation, Aufomafica 28 (5) ( 1992) 273-283. 19 1 D. Dochain, Contribution to the analysis and control of distributed parameter systems with application to (bio)chemical processes and robotics, Thesis, UniversitC Catholique de Louvain, Belgium, 1994. I IO] A. Friedman, Parfial Di$erential Equafions (Krieger, New York, 1976). I I I I J.A. Goldstein. Semigroups of Linear Operators and Applications (Oxford University Press, Oxford, 1985). I 12 1 J.A. Goldstein and J. Tervo, Existence of solutions for a system of nonlinear partial differential equations related to bioreactors, Dynamic Syst. Appl., to appear. [ 131 J.W. Grizzle, M.D. DiBenedetto and F. Lamnabhi-Lagarrigue. Necessary conditions for asymptotic tracking in nonlinear systems, IEEE Trans. Automat. Contr. 39 ( 1994) 1782-1794. I 141 H.P. Helfrich, Error estimates for semidiscrete Galerkin type approximations to semilinear evolution equations with nonsmooth initial data, Numer. Math. 51 (1987) 559-569. I IS I K.S. Hong and J. Bentsman, Direct adaptive control of parabolic systems: Algorithm synthesis and convergence and stability analysis, IEEE Trans. Automat. Contr. 39 ( 1994) 2018-2033. I 16 I A. Isidori, Nonlinear Control Systems (Springer, Berlin, 1989), I 17 I T. Kate, Nonlinear evolution equations in Banach spaces, Proc. Symp. Appl. Math. 17 ( 1964) 50-67. ( 18 I A.W. Lang, Systems of Nonlinear Partial D@erential Equations: Applications to Biology and Engineering (Kluwer, Dordrecht, 1989).
[ I I H. Amann, Existence and regularity
I 19 1 Y. Lin and T. Zhang, Finite element methods for nonlinear Sobolev equations with nonlinear boundary conditions, J. Math. Anal. Appl. 165 (1991) 180-191. [ 20 I J.L. Lions, Optimal Control of Sysfems Described by Partial Di$erenfial Equations (Springer, 1971).
Berlin,
216
M.7: Nihtila et al./Simulation Practice and Theory 5 (1997) 199-216
[ 211 M.T. Nihtill, N. Tali-Maamar and J.P. Babary, Controller design for a nonlinear distributed parameter system via a collocation approximation, in: M. Flies% ed., Proceedings of the IFAC Symposium on Nonlinear Control System Design, NOLCOS’92, Bordeaux, France (1992) 447-452. [22] M.T. Nihtila, J. Tervo and J.P. Kaipio, Simulation of a nonlinear distributed parameter bioreactor by FEM approach, Dept. of Computer Science, Applied Mathematics Report Series A/ 1994/5, University of Kuopio, 1994. [ 231 A. Pazy, Semigroups of Linear Operators and Applications to Partial Dt’erential Equations (Springer, Berlin, 1983). [24] J. Smoller, Shock Waves and Reaction-Dtpsion Equations (Springer, Berlin, 1983). [25] N. Tali-Maamar, T. Damak, J.P Babary and M.T. Nihtill, Application of a collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simul. 35 (1993) 303-319. [26] H. Tanabe, Equations of Evolution (Pitman, London, 1979). [27] V. Thomee and L. Wahlbin, On Galerkin methods in semilinear parabolic problems, SIAM J. Numer. Anal. 12 ( 1975) 378-389. [ 281 V. Thomee, Galerkin Finite Element Methods for Parabolic Problems, Lecture Notes in Mathematics, Vol. 1054 (Springer, Berlin, 1984). [29] J.V. Villadsen and Michelsen, Solution of Dtfirential Equation Models by Polynomial Approximation (Prentice-Hall, Englewood Cliffs, NJ, 1978).