3.7 Reachability of fed-batch fermentation processes
3.7.3 Reachability analysis
=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
deniteness of itstime-derivative:
dV
dt
=
@V
@x _ x=
@V
@x
f(x)
where
f(x)=f(x)in the open loop case (assuming zero input) and
f(x)=f(x)+
g(x)C(x)intheclosedloopcasewhereC(x)isthestaticlinearornonlinearfeedback
law.
The diagonalweighting matrixQinthequadraticLyapunov functionisselected
in a heuristic way: a state variable which does not produce overshoots during the
simulation experiments gets a larger weight than another state variable with
over-shooting behaviour. In the new norm dened by this weighting, a more accurate
estimateofthestabilityregioncanbeobtained. Withthisanalysiswecannot
calcu-late theexact stability regionbut theresults give valuableinformationforselecting
the controller type and tuning its parameters. The nonlinear stability analysis
re-sults inthe time derivativeof the quadratic Lyapunov functionasa function ofthe
Q=I
state variables,whichisa two variate functioninour case seen ing. 3.5. The
sta-bility region of the open-loopsystem is the region on the (x
1
;x
2
) plane over which
thefunctionisnegative. WeremarkthatthisheuristiccandidateLyapunovfunction
selection method serves mainlyfor illustrational purposes and is not considered as
ascientic contribution.
3.10 Summary
The main contributions of this chaptertothe analysis ofnonlinear process systems
are as follows:
Using the special structural properties of process models, the generally
com-putationallycomplexnonlinear reachability analysis can be performedanalytically.
Furthermore, in sections 3.6 and 3.7 it was shown that the singular points of the
reachabilitydistributionhaveclear physicalmeaning. Theseresults wereillustrated
using the models of continuous and fed-batch fermentation processes.
In section3.7rigorousnonlinear analysiswasused foranalyzingthereachability
propertiesofa simplefed-batchfermenter modelandto relatethemtothe
physico-chemical phenomena taking place in the reactor. With a help of this grey-box
approachwehaveshownthatthe knowndicultiesofcontrollingsuchprocessesare
primarilycaused by the factthat the rank ofthe reachabilitydistribution isalways
two which is less than the number of state variables being three. Furthermore, a
istics that shows the nonlinear combination of the state variables that is constant
independently of the input. The coordinates transformation is independent of the
most uncertain part of the state space model: the source () function, too. The
results are extended to the four state variable non-isotherm case, and to nonlinear
fed-batch chemical reactors with generalreaction kinetics. The rank of the
reacha-bilitydistributionremainstwointhis casegivingrise totwoconserved combination
of the state variables independently of the input. The structural properties of the
process models enabling toapply the proposed analytical technique have alsobeen
described.
In section 3.8 it was shown on the example of continuous bioreactors that the
notion of zero dynamics is very useful for output selectionfor control purposes. It
wascalculated thatthe best controlled outputchoice isthe substrate concentration
since the system is minimum-phasewith respect to this output.
Analysis Based Control Structure
Selection
4.1 Motivation
The control of nonlinear process systems is a challenging and emerging
interdis-ciplinary eld of major practical importance. The most common way to control
nonlinear process systems is either use linear techniques on locally linearized
ver-sions of the nonlinear models or use model-based predictivecontrol[76].
At the same time a number of powerful and theoretically well grounded tools
andtechniques areavailablefornonlinear controlinthe eldofsystems and control
theory ([35], [89]) which are applied successfully in other application areas. These
techniques, however, require most often symbolic computation and may be
non-feasible for real process systems. This may be one of the reasons they are not
known and not appliedextensively for process systems.
Fermentationprocessesinparticularexhibitstrongnonlinearcharacteristicsand
areknown tobediculttocontrol. Theinvestigated simplefermentationprocess is
thereforeused asabenchmark problemforadvancednonlinear analysis and control
techniques. Many authors have examined the various approaches of analyzing and
controlling fed-batch ([45], [9], [39], [92]) and continuous fermentation processes
([86], [44], [87]).
Nonlinear controllersof dierent type are designed and compared on the
exam-ple a simple fermenter near an optimal production operating point which is close
toits foldbifurcationpoint. Nonlinearanalysis of stability,controllabilityand zero
dynamics presented in the previous chapter is used to investigate open-loop
sys-tem properties, toexplore the possible controldiculties and to design the system
outputtobeused forcontrol. A widerange ofcontrollersare testedincluding
pole-placement and LQ controllers, feedback and input-output linearization controllers
and anonlinear controller based ondirect passivation. The comparison isbased on
time-domainperformance and oninvestigating the stability region, robustness and
tuningpossibilitiesof the controllers.
ysis
Variouscontrollersofdierenttype: pole-placementcontroller,LQcontrollers,
feed-back linearization based controllers and nonlinear controllers based on direct
feed-back linearization based controllers and nonlinear controllers based on direct