3.6 Reachability of continuous fermentation processes
3.6.2 Controllability analysis
0 0:4011
1:6045 1:2033
3.6.2 Controllability analysis
Controllability properties of the system play key role in designing controllers not
onlyin the desired operating pointbut alsofor the entire operatingregion one can
foresee. Therefore, nonlinear analysis techniques are recommended to complement
the usual analysis based onlocallylinearized models.
Analysis based on local linearization Aftercalculatingthe
Kalman-control-labilitymatrixofthelinearizedmodelwendthatthesystem iscontrollable(inthe
linear sense) in a close neighborhood of the required operating point, because the
controllability matrix has rank 2. However, even in an important operating point
with a maximum reaction rate very near to it S
controllabilitymatrixhas onlyrank1,thatis,the systemislocallynotcontrollable.
tributionsdescribedinsection3.2.3isthenusedforidentifyingthesingularpointsof
the state space around whichcontrolof the system is problematic oreven
impossi-ble. Thelocalreachabilitydistributionisgenerated incrementallyinanalgorithmic
way [35] intwosteps as follows. The initialdistribution is
0
=spanfgg (3.36)
This is extended by Lie-bracket of f and g in the next step
The second step then gives the following
The Lie-product of f and g isas follows
[f;g](
SS0+K2S 2
Singular points At the point
all the elements of
2
(and, of course, all the elements of
0
and
1
) are equally
zero. Itmeans,thatthereachabilitydistributionhasrank0atthispoint. Moreover,
thissingularpointisasteadystatepointinthestatespace. Fromthisitfollows,that
if the system reaches this (non-desired) point, it's impossible to drive the process
out of it by manipulating the input feed ow rate. If
has rank 1. From a practical point of view it means
that if the biomass concentration decreases to 0
g
l
then it can't be increased by
changing the input ow rate. Both of the above singular points are "wash-out"
states in a sense, as the biomass concentration decreases to zero. Therefore, these
are undesirablestates. It'seasyto calculatefromEq. (3.40) thatinadditiontothe
previous results
1
alsoloses rank atthe points
butifwecalculatetheLie-products[f;[f;g]]and[g;[f;g]]wendthatthesesingular
points "disappear" but the previous ones donot.
Non-singular points At any other point in the state space including the
desired operating point
the reachability distribution has rank 2,
whichmeans,that thesystem isreachableinaneighborhoodofthesepointsand we
can apply state feedback controllersto stabilizethe process.
The reachability of a simple nonlinear fed-batch fermentation process model is
in-vestigated in this section. It is shown that the known diculties of controlling
such processes are primarily caused by the fact that the rank of the reachability
distributionis always less than the number of state variables.
Furthermore, acoordinates transformationis calculatedanalytically that shows
the nonlinear combinationof the state variables which is independent of the input.
The results of the reachability analysis and that of the coordinates transformation
are independent of the source functionin the system model.
The results are extended to the four state variable non-isotherm case, and to
nonlinear fed-batch chemical reactors with generalreaction kinetics.
3.7.1 Problem statement
Bio-processes in general and fermentation processes in particular are dicult to
model and tocontroleven inthe simplestcases. Variousdiculties are reported in
theliterature whichincludeinstabilityandcontrollabilityproblems forboth
contin-uous and fed-batch fermenters ([45], [44]).
The dynamic state space modelof a fermenter is derived from rst engineering
principles whichxes certainstructuralelementsinthe model. Thestate equations
are derived from dynamic conservation balances of the overall mass, component
masses and energy if applicable. The speciality of a fermentation model appears
in the so called source function of these balances which is highly nonlinear and
non-monotonous innature.
Theaimofthissectionistouse rigorousnonlinear analysisofasimplefed-batch
fermenter model for analyzing itsreachability properties and to relate them to the
physico-chemicalphenomena taking place inthe reactor.
3.7.2 Nonlinear state space model
The simplestdynamicmodelof afed-batch fermenter consistsof threeconservation
balances for the mass ofthe cells (e.g. yeast to beproduced), that of the substrate
(e.g. sugarwhichisconsumedbythecells)andfortheoverallmass. Hereweassume
thatthefermenter isoperatingunderisothermconditions,that isnoenergybalance
is needed. The cellgrowth rate is described by anonlinear staticfunction .
Initially,asolutioncontainingbothsubstrateandcellsispresentinthefermenter.
Duringthe operationwefeedasolutionofsubstratewithagiven feedconcentration
S
f
tothe reactor.
Undertheaboveassumptionsthenonlinearstatespacemodelofthefermentation
process can be writtenin the followinginput-ane form[45]
_
x=f(x)+g(x)u (3.41)
where
x= 2
4 x
1
x
2
x
3 3
5
= 2
4 X
S
V 3
5
; u=F (3.42)
f(x)=
The variables of the modeland their units in square brackets are
x
1
=X biomassconcentration (state) [g/l]
x
2
=S substrate concentration (state) [g/l]
x
3
=V volume (state) [l]
u=F feedow rate (input) [l/h].
The constant parametersand their typicalvalues are the following
Y =0:5 yieldcoecient
max
=1 maximum growth rate [h
1
]
K
1
=0:03 kinetic parameter (Monodconstant) [g/l]
K
2
=0:5 kinetic parameter [l/g]
S
f
=10 inuent substrate concentration [g/l]
3.7.3 Reachability analysis
We construct the reachability distribution according to the algorithm described in
section3.2.3.
0
=spanfgg (3.45)
The calculation of the Lie-products in
1
and
2
isas follows.
[f;g](x)=
the Lie-product[f;g]has the form
2
i
It follows from Eqs. (3.48)-(3.51) that the distributions [f;[f;g]] and [ g;[f;g]]
willalsohave the same formas (3.51), i.e.
[f;[f;g]]=
On the basis of the above we can denote the coordinate functions of the vector
elds spanning
2
ata given pointx of the state space as follows
whichmeansthatwecouldn'tincreasethedimensionofthereachabilitydistribution
inthe secondstep and the rankof
2
is atmost 2 inany point of the state space.
Singular points There are, however, pointsin the state space where the rank of
the reachability distribution
2
is of dimension 1. This case means that there is no
biomass in the system and since the inlet ow contains only substrate, the
biomassconcentration cannotbeinuenced by manipulatingthe input.
During the following analysis we will consider the open region of the state space
where
1
isnonsingularandthe valueofstatevectorhasreal physicalmeaning(the
concentrations and the liquidvolume are positive) i.e.
U =fx
Since the generation of the reachability distribution stopped inthe second step
1
=span fg;[f;g]g
is the smallest distribution invariant under f;g and containing the vector eld g.
This distribution is denoted by hf;gjspanfggi. Since hf;gjspanfggi is nonsingular
on U and involutive we may use it to nd a coordinates transformation z = (x).
Thesystem inthenew coordinateswillberepresented byequationsof thefollowing
form(see Theorem 3.2.1)
_
inour case.
To calculate , we have to integrate the distribution
1
rst, that is to nd
a single (dim(x)-dim(
1
)=3-2=1) real valued function such that span fdg =
[hf;gjspanfggi]
?
,where the sign? denotes the annihillatorof adistribution. Since
[f;g](x)=
this amounts tosolvethe partial dierentialequations (PDEs)
Solution by the method of characteristics
The method of characteristics (see e.g. [8], [3]or[66]) isused for solving the above
resulted rst order linearhomogeneous partialdierential equationinthe following
generalform:
n
orbriey
(x) 0
(x)=0; (3.66)
where T R is a domain, x 2 T,
i
;i = 1:::n are known functions and is
the unknown. The characteristic equation system of (3.66) is the following set of
ordinary dierentialequations:
_
=(): (3.67)
We call the : R ! R n
solutions of (3.67) characteristic curves. A 2 C 1
(T)
function is called the rst integral of (3.67) if t ! ((t)) is constant along any
characteristic curve. In order to solve (3.66) we have to nd (n 1) linearly
inde-pendentsolutions(
1
;
2
;:::;
n 1
)ofit. Thenthegeneralsolutionof(3.66)willbe
in the form = (
) is an arbitrary function.
We know that a rst integral of (3.67) satises (3.66), therefore we have to nd
(n 1) linearly independent rst integrals toobtain the general solution. This can
be done withoutsolving (3.67) asit is illustratedbelowin our case.
To solve the rst PDE,namely
@
we start fromthe following set of ordinary dierentialequations:
_
It's easyto observe that
_
=const. Moreover,
_
fromwhichit follows that
x
We can see fromthe above that the solution of (3.68)will be inthe form
(x
with anarbitrary C 1
function.
To solvethe secondPDE werst rememberthatinthe reachability distribution
Æ
@
The characteristic equationsare writtenas
_
It's easyto see that
1
Therefore the solution of (3.69) isin the form
(x
with an arbitrary C 1
function
. To give a commonsolution for both (3.68) and
(3.69) wepropose the function
(x
from which we can see that it indeed satises both PDEs. With the help of we
can dene the local(and luckily global)coordinates transformation :R n
the transformed form of the model(3.41)-(3.43)can be writtenas
_
The aim of this section is to show the reasons present in the original state space
model which led to the reachability and state transformation above. This analysis
enablesto nd other models of similarform with the same properties.
Physical analysis of the model and the solutions
The rst importantthing to observe is that the results of the reachability analysis
and thatofthe coordinates transformationdonot depend onthe actualformofthe
functioninEq. (3.44). Theresultsutilizedthefollowingspecialitiesoftheoriginal
state space model(3.41)-(3.43).
(i) the constant coecients inthe 3rd state equation,i.e.
f
3
=0 ; g
3
=1
where f
i and g
i
are the ithentry of the vector functions f and g in the state
space model. This property always holds for the overall mass balanceof
fed-batch reactors.
(ii) the relationbetween the 1st and the 2nd state equation, namely
f
2
= 1
Y f
1
=C
f f
1
where C
f
is a constant. Such a relationship exists if the two related state
variables, x
1
and x
2
are concentrations of components related by a chemical
reactioninthe form 1
Y
S !X [23].
Further we may notice that the quantity in Eq. (3.71)which is conserved
independently of the input consists of two parts corresponding to the substrate
mass and cellmass of the system asfollows:
(x
1
;x
2
;x
3
)=V(S
f
S)+ 1
Y V(X
f
X) (3.75)
with X
f
= 0 because the feed does not contain any cells. The above two terms
originatefromthe (weighted)convectivetermsinthe componentmassconservation
balancesrespectively,that issuchtermswhichareonlycausedbythefeedasinow.
Generalized state space models
We can generalize the original modelin Eqs. (3.41)- (3.43) in two steps if we want
topreserve the special dynamic properties of the model.
1. Generalreaction rate function
As the results do not depend onthe function in Eq. (3.44),we can replace
the fermentationreaction by a generalchemical reactionof the form
1
Y
S ! X
wherethe reactionrate(source) functionis
(x
2 )x
1
with
isanunspecied
possibly nonlinear function.
If we further release the assumption that the fermenter is operating under
isothermalconditions,thenweshouldincludethe energyconservationbalance
totheoriginalmodel. Thenafourstatemodelisobtained[23]inthefollowing
input-aneform:
_
with T being the temperature inthe fermenter.
f
and e.g. the followingadditionalconstant parameters
c
2
reactionenthalpy coecient [m 3
K/J]
T
f
=293 inuenttemperature [K]
Observe,thatnowthereactionratefunction
dependsalsoonthetemperature
x
3
=T givingrise tothe source function
Furthermore, the required structural properties (i) and (ii) are present in the
generalizedmodel. The property (i)nowholdsfortheentries f
4 andg
4
whichisthe
overallmassbalance. Therearetwoindependentpairs,(f
1
;f
2
)(the twocomponent
mass balances)and (f
1
;f
3
)(amass and anenergybalance) whichpossess property
(ii)with dierentconstants.
Analysis of the generalized models
In the above four state variable case the nal reachability distribution after four
steps would be the following
=span fg
Ifwe calculatethe Lie-products [f
[g
Therefore thecalculationof thereachability distributionstopshere anditturnsout
that the dimension of the distribution is 2in this case.
To nd the decomposed system similarlyto(3.74)we haveto nd two
indepen-dentreal-valued functions
1
and
2
suchthat
It's easyto check that the twoindependent functions
satisfy the PDEs in eq. (3.81). Therefore the new coordinate vector z is given by
the function
and the system (3.76)-(3.78) inthe new coordinates is writtenas
_
with the condition z
3 6=0.
ditionsx
1
(0)=2 g
l , x
2
(0)=0:5 g
l ,x
3
(0)=0:5 g
l
3.7.6 Engineering interpretation
The invariance of in eq. (3.71) expresses the fact that the state variables of the
fed-batch fermenter model can only move on a smooth hypersurface in the state
space. The shapeof thishypersurfaceobviouslydepends onthe choiceof the initial
values of the state variables. It means that the initial concentrations and liquid
volume (that are set by the control engineer) uniquely determine the set of points
inthe state space that are reachable during the process. Figs 3.2 and 3.3illustrate
the eect of the initial liquid volume on the reachability hypersurface when the
concentrations are xed. It is shown that if the initial volume is too small then
thepossibilitiestocontrolthe biomassconcentrationx
1
aredramaticallyworsening.
Similarly,the eect of the initialconcentrations can alsobe easilyexamined, since
in(3.71) is aquite simple functionof the statevariables.
Withthe help of controllerdesign becomes easier. If the desirednal point of
the fermentation is given (in the state space) then the initialconditions can be set
insuch a way that the desired point is reachable.
3.7.7 Comments on observability
Due to space limitations we cannot go into details concerning the observability of
fed-batch fermentation processes, we can only briey describe the most interesting
aspect of the relation between reachability and observability. The rough problem
ditionsx
1
(0)=2 g
l , x
2
(0)=0:5 g
l ,x
3
(0)=0:1 g
l
statement of observability is the following: is itpossible todetermine the values of
the state variables of the system if we measure the inputs and the outputs?
Ob-viously, the observability property of linear and nonlinear systems largely depends
on the selection of the output function y = h(x). Let us suppose that the output
of the system (3.41)-(3.43) is chosen to be ineq. (3.71). It's clear without
com-plicatedcalculations, that the system won't be observable because is constant in
timeindependentlyofthe inputuandthereforeitdoesnotprovideanyinformation
abouttheinternal"movement"ofthesystem. It'svalueonlyidentiesthereachable
hypersurface (manifold).
3.7.8 The minimal realization of fed-batch fermentation
pro-cesses
Using the calculated function, it's not dicult to give a minimal state space
realizationoffed-batchfermentationprocessesinthetemperature-independentcase.
Sincethereachabilityhypersurfacedenedbyandshownings3.2and3.3is
two-dimensional,the minimalrealizationwillcontaintwostatevariables(i.e. the
input-to-state behaviour of the system can be described by two dierential equations).
Since isconstant in time, it'sclear that
(x(t))= 1
Y x
1 (t)x
3
(t) (x
2 (t)x
3
(t) S
f x
3
(t))= (3.88)
=
Thereforewecanexpresse.g. thevolumex
3
fromtheaboveequationinthefollowing
way:
and the minimal state spacemodel reads
_
It'swell-known fromsystemtheory that state-spacerealizationsare not uniqueand
it'seasy to see that insteadof x
3
any one of the other two state variables could be
expressed from eq. (3.88). Therefore one can select those two state variables that
are importantfroma certainpointof view(e.g. acontrolproblem) andexpress the
third one fromeq. (3.88).
It'salsoimportanttoremarkthat themodel(3.91)has aspecialstructure, since
it contains the initialvalues of the original model (3.41)-(3.43) in the input vector
eld g
min
(but luckily not in the vector eld f
min ).
3.8 The zero dynamics of continuous fermentation
processes
Inordertoanalyzezerodynamicsasitisdescribed insection3.4,weneedtoextend
the originalnonlinear state equation(3.26) with anonlinear outputequation
y=h(x) (3.92)
where y is the output variable and h is a given nonlinear function. Then the zero
dynamics of aninput-ane nonlinear system containing two state variables can be
analyzedusing a suitablenonlinear coordinates transformationz =(x):
where (x)is a solutionof the followingpartial dierential equation(PDE):
L
g
(x)=0 (3.94)
where L
g
solve the above equationto obtain:
where F is an arbitrary continuously dierentiable function. Then we can use the
simplestpossiblecoordinates transformationz =(x) inthe following form:
3.8.1 Selecting the substrate concentration as output
Ifa linear functionof the of the substrate concentration is chosen asoutput, i.e.
z
wherek
s
isanarbitrarypositiveconstantthentheinversetransformationx= 1
(z)
isgiven by
z2X0ks S
f
Thusthe zero dynamicsin the transformed coordinates can becomputed as
_
which gives
_
since by construction L
g
(x)=0 (see eq. 3.95). The aboveequation is constrained
by y = k
= 0. Then the zero dynamics of the system is given by the
dierential equation
_
which islinear and globally stable. The equilibriumstate of the zero dynamicsis at
z
2
= V
Y
which together with z
1
= 0 corresponds to the desired equilibrium state
x
1
= 0;x
2
= 0 in the original coordinates. The above analysis shows that if we
manage to stabilize the substrate concentration either by a full state feedback or by
an output feedback (partial state feedback) or even by a dynamic controller (which
does not belongto the scope of this chapter) then the overall system willbe stable.
3.8.2 Selecting the biomass concentration as output
The outputin this case is a linearfunction of the biomass concentration:
z
_
which isonlylocallystable aroundthe desiredequilibriumstate and the right hand
side of Eq. (3.104) has singular points (where the denominator is 0). The stability
region is independent of k
x
and can be determined using the parameters of the
system.
3.8.3 Selecting the linear combination of the biomass and
substrate concentrations as output
In this case the outputis the linear combinationof the biomass and substrate
con-centrations:
In this case the zero dynamics is alsolocallystable around the desiredequilibrium
state and it again has singularpoints. Furthermore, anew non-desiredequilibrium
state appears at z
2
= Vkx
k
s
which can be inside the operating region depending on
the values of k
x and k
s
. The results of the analysis of the zero dynamics show that
thebest choice ofoutput tobe controlled isthesubstrate concentrationand involving
the biomass concentration into the output generally brings singular points into the
zerodynamics andmakes thestabilityregion narrower. Thesetheoreticalissueswill
help inunderstanding the following simulationresults.
3.9 Stability analysis of continuous fermentation
pro-cesses
As localstability analysis shows, in a neighborhood of the desired operating point
thesystem isstable,butbecausethe pointisveryclosetothefoldbifurcationpoint
(X
this stability region is small. This is illustrated in g. 3.4 which shows that the
system moves to the undesired wash-out steady state when it is started fromclose
neighborhoodof the desired operatingpoint(X(0)=4:7907 g
Stability analysis based on local linearization aroundthe operating point
depends uponthe eigenvalues of the linearized state matrix A inEqs. (3.32)-(3.33)
which are complex conjugate values inour case:
12
= 0:60170:5306i (3.106)
We can see that the process is indeed stable around the operating point but the
linear analysis doesnot giveany informationonthe extent of the stability region.
0 20 40 60 80 100 120 0
1 2 3 4 5 6 7 8 9 10
Time [h]
Concentration [g/l]
Biomass concentration Substrate concentration
Figure3.4: Open loopbehaviorof the system
Nonlinear stability analysis is based on Lyapunov technique which aims at
ndingapositivedenitescalar-valuedgeneralizedenergy function V(x)whichhas
negativedenitetimederivativeinthewholeoperatingregion. Mostoftenageneral
quadraticLyapunov functioncandidate is used in the formof
V(x)=x T
Qx
with Q being a positive denite symmetric quadratic matrix, usually of diagonal
form. This function is scalar-valued and positive denite everywhere. The
stabil-ity region of an autonomous nonlinear system is then determined by the negative
stabil-ity region of an autonomous nonlinear system is then determined by the negative