Yield Trajectory Tracking for Hyperbolic
AgeStructured Population Systems\thanksreffootnoteinfo
Abstract
For population systems modeled by agestructured hyperbolic partial differential equations (PDEs) that are bilinear in the input and evolve with a positivevalued infinitedimensional state, global stabilization of constant yield set points was achieved in prior work. Seasonal demands in biotechnological production processes give rise to timevarying yield references. For the proposed control objective aiming at a global attractivity of desired yield trajectories, multiple nonstandard features have to be considered: a nonlocal boundary condition, a PDE state restricted to the positive orthant of the function space and arbitrary restrictive but physically meaningful input constraints.
Moreover, we provide Control Lyapunov Functionals ensuring an exponentially fast attraction of adequate reference trajectories. To achieve this goal, we make use of the relation between firstorder hyperbolic PDEs and integral delay equations leading to a decoupling of the inputdependent dynamics and the infinitedimensional internal one. Furthermore, the dynamic control structure does not necessitate exact knowledge of the model parameters or online measurements of the ageprofile. With a Galerkinbased numerical simulation scheme using the key ideas of the KarhunenLoèvedecomposition, we demonstrate the controller’s performance.
Stuttgart,SanDiego]Kevin Schmidt Athens] Iasson Karafyllis SanDiego] Miroslav Krstic
Institute for System Dynamics, University of Stuttgart, Waldburgstr. 17/19, 70569 Stuttgart, Germany
Department of Mathematics, National Technical University of Athens, Zografou Campus, 15780, Athens, Greece
Department of Mechanical and Aerospace Engineering, University of California at San Diego, La Jolla, CA 920930411, USA
Key words: firstorder hyperbolic PDE, chemostat, tracking control, timedelay systems, input constraints, Galerkin methods.
1 Introduction
We design an asymptotically tracking control for agestructured chemostats modeled by hyperbolic partial differential equations (PDEs) with a bilinearly acting input.
Based on our prior work on the stabilization of constant yield set points, we guarantee global attractivity of output trajectories with an exponential convergence rate. In addition, we developed an efficient numerical scheme based on Galerkinmethods, which ensures accurate asymptotic properties.
Motivation. In the context of mathematical biology and demography agestructured continuoustime models are a common way of describing the evolution of a certain population with respect to the independent variables of age and time [3, 4]. Continuous bioreactors encountered in bioengineering and pharmaceutical research are usually modeled by ordinary differential equations (ODEs) [15, 16]. Since multiple ecological concepts like resistance and resilience of ecosystems are closely related to the framework of robustness in system theory, these aspects have been studied rigorously [6, 10]. An analysis of the ergodicity problem is given in [8, 9]. Moreover, for agestructured models there is an extensive literature on optimal control problems [3, 7] as well as on the stability of certain PDE models [15].
The model. Throughout this paper we focus on the McKendrickvonFoerster PDE, which is introduced in Section 2. For this setting the dilution rate, which is the ratio of the volumetric flow to the constant volume in the growth chamber, is a natural control variable [16, 17]. On the other hand, the output is chosen as a weighted integral of the population’s agedistribution. As a result, it is possible to represent all products which are proportional to the overall population as well as possibly agedependent synthesized products. Furthermore, the dependence of the microorganisms’ growth rate on the nutrient concentration in the bioreactor is not captured in the model and we hence assume that the nutrient concentration is maintained constant.
Timevarying yields. Having a biomass in mind which is used for the production of antibiotics in pharmaceutical industry, the relevance of the trajectory tracking issue can be elucidated in a demonstrative way. Choosing the yield of an antibiotic as a valid output, its production rate is subjected to external influences like seasonal effects or the current demand. With the prediction of these exogenous factors on an adequate time scale it is possible to determine an optimal production rate governing desired yield trajectories. Due to the fact that a periodic reactor operation may produce a higher yield than the yield achieved by an equilibrium point [2], a special focus is placed on periodical reference trajectories. Beyond this, if the considered bioreactor is a part of a cascaded process, the tracking of predefined trajectories makes it possible to accelerate starting processes and changes of operating points.
Results of the Paper.
The present paper extends our prior work [12], which focused on the global stabilization of desired equilibrium points of the system class under consideration. We now aim at ensuring the global attractivity of desired yield trajectories and therefore generalize the already established concepts, such that constants set point are included as a special case. For this purpose, a definition of the control objective is given in Section 3. The suggested approach exploits the relation of firstorder hyperbolic PDEs to delay models in the sense of integral delay equations (IDEs, see [11]). More specifically, we decompose the PDE problem to an inputdependent finitedimensional subsystem and an autonomous delay subsystem which is correlated to the microorganisms’ reproduction. For special cases of integral kernels we are in the position to construct Control Lyapunov Functionals (CLFs).
In contrast to the prior work, we use a twodegreesoffreedom (DOF) control structure with a separate feedforward control part evoked by the reference trajectory [14], as introduced in Section 4. In this case the feedforward controller does not solely enhance the tracking behavior of the closed loop, but plays an essential role in the overall attractivity concept. The consideration of input constraints is a crucial issue of the present control design assuming a bounded interval for the accessible dilution rate. Moreover, our outputfeedback controller does not demand online measurements of the population’s entire agedistribution. Even the knowledge of exact system parameters is not necessary, since the controller handles uncertainties in a robust way.
In addition, it is important to guarantee that the PDE state, which represents the population density, remains positive at all times and ages. This fact, in conjunction with a control input directly acting on the whole profile (not simply on the boundary), differentiates our work to other control problems of hyperbolic PDEs [1, 5].
Lastly, we provide simulation results of the closedloop system in Section 6 with a Galerkinbased simulation scheme, which conserves important system properties even at low orders and enables independent age and time discretizations.
Notation.

The set denotes all positivevalued real numbers, all nonnegative real numbers.

The inner product of is denoted where .

is the maximum resp. norm for .

is the class of all strictly increasing, unbounded functions with .

The saturation function with respect to is defined ; other intervals are explicitly denoted as an index.

Given the functions , with the metric space , we define the right temporal Diniderivative .

For any and , denotes the class of all functions for which there exists a finite (or empty) set such that: (i) the derivative exists at every and is a continuous function on , (ii) all meaningful right and left limits of when tends to a point in exist and are finite.
2 AgeStructured Population Models
Consider the McKendrickvonFoester PDE (2) valid in the agetime domain
(1)  
(2)  
(3) 
which describes the evolution of the population density as a part of an initialboundary value problem (IBVP) on the same domain with an arbitrary large but finite maximum reproductive age . Strictly speaking, the state describes the density of the overall population which has reached a specific age at a certain time . In addition, the function denotes the agedependent mortality rate and the dilution rate which is the control input. In particular, the nonlocal boundary condition (BC) (2) is valid for and models the production of newborn individual determined by the birth modulus resp. the kernel . Furthermore, (2) is the initial condition, i.e. the initial distribution of the population density in the agedomain at . In addition, the output is defined by the equation
(4) 
which takes the possibly agedependent production rate of a specific (bio)chemical species into account. For instance, the choice includes the overall population as a valid output.
The distributed parameter system given by (2)–(2) with the input and the output is of bilinear singleinputsingleoutput type. Subsequently, we introduce three assumptions in order to guarantee the existence of a meaningful unique solution of the IBVP (2)–(2) aware of the state and input constraints (see also [12]):

The parameters functions are restricted to and , where .

The control takes values in where .
3 Control Objective and Delay Model
The asymptotic tracking of a reference trajectory with respect to the output given by (2) defines the key objective of the contribution. For designing an asymptotic tracking control in the context of system , consider the 2DOF control structure satisfying assumption (A) a priori
(5) 
consisting of a feedforward and a feedback component as sketched in Fig. 1. More specifically, the desired global asymptotic attractivity w.r.t. for all valid initial conditions is defined as follows
(6) 
There are basically two different types of relevant reference trajectories in the context of chemostats. The first type focuses on the planning of a trajectory which reaches a set point in a desired transition time starting from an arbitrary initial profile and is motivated by starting processes or set point changes during an operating chemostat. In addition the stabilization of periodical trajectories is relevant due to an efficient harvesting of synthesized products [13] and is considered as a second case.
Before starting with the control design, we analyze the equilibrium profiles and introduce a delay system to describe the dynamics of system (2)–(2). Therefore, define the equilibrium dilution rate as the unique solution of
(7) 
It follows from a steady state analysis of : (2)–(2) that the necessary condition corresponds to a continuum of equilibrium profiles [12], which are proportional to
(8) 
As a prerequisite for the following proposition, we use this fact to define the kernel and a reference profile of the PDE state with the property :
(9)  
(10) 
In the following proposition we refer to the openloop system, i.e. it is valid for any piecewise continuous input satisfying assumption (A). Its proof is omitted, since it follows the arguments of the proof of Lemma 3, which addresses the controlled system (see Section 5)
Proposition 1 (delay model)
Define the functional
(11) 
where is the continuous function
(12) 
Given an IC for the PDE system (2)–(2), we define
(13) 
For every piecewise continuous input , consider the solution of the delay model
(14)  
(15) 
with the ICs and for . Define the functions for ,
(16)  
(17) 
Then, the unique solution of the PDE model (2)–(2) corresponding to the input and the IC is given by for . Moreover, the output is given by for .
With (13)–(15) we decompose the PDE dynamics to an inputdependent ODEstate and an infinitedimensional internal coordinate which cannot be affected by the input . It can be seen from (17), that a steady state of the delay model , implies an exact output tracking, i.e. . With the definition
(18) 
we rewrite (17) and get following result for the logarithmic tracking error as the basis for the control design
(19) 
4 Nonlinear Controller Design
For the design of an appropriate controller, we state a second set of assumptions:

The equilibrium dilution rate defined by (7) satisfies .

Assume that belongs to the class of valid reference trajectories , which is defined as all positivevalued and continuously differentiable functions with the additional property
(20) 
We assume the existence of a constant which fulfills the inequality with defined by (7)
(21) 
The system parameters as well as the initial profile are assumed to be unknown. As a result of (7), the equilibrium input is also unknown.
While assumption (B) ensures the reachability of the plant’s steady states, the extra condition (20) involved in assumption (B) is needed to guarantee the asymptotic tracking in the presence of input constraints. Note that if the birth modulus has a finite number of zeros in the agedomain , assumption (B) holds true by virtue of Proposition 2.3 in [12], no matter what and are. Moreover, this assumption is not necessary to prove stability but to explicitly construct CLFs.
Consider the nonlinear dynamic output controller specifying the 2DOF structure in (3) with the control gain and the observer gains for
(22)  
(23)  
(24)  
(25) 
The proposed controller (22)–(25) consists of three components: Firstly, the feedforward control is constructed by means of the scalar dynamics (14). On the other hand, the feedback part (23) consists of a proportional control of the logarithmic error (19) and an observerODE (24) which adapts the controller to the unknown equilibrium dilution rate . Hereby, the observer state has the IC and its second component is the initial guess for the equilibrium dilution rate . For and , has a steady state . As a result we define the observer error
(26) 
The control law (22)–(25) guarantees an attractivity with exponential convergence of an admissible reference trajectory as characterized by the subsequent theorem:
Theorem 2 (robust global attractivity)
Consider the system (2)–(2) under assumptions – and –. For every , there exists a constant and a function (independent of ), such that the unique solution of the closed loop (2)–(2) with the controller (22)–(25) satisfies the following estimate for any and and all :
(27) 
Moreover, for any constants satisfying
(28) 
there exists a sufficiently small constant , sufficiently large constants , such that for every the functional
defined by
(29) 
with
(30) 
is a Lyapunov functional for the closed loop system (2)–(2), (22)–(25). More specifically, the differential inequality
(31) 
holds for all and every solution of the closed loop (2)–(2) with the controller (22)–(25).
Theorem 2 does not only guarantee the global asymptotic attractivity of the output , but it determines the whole PDE state of the closedloop system to follow the reference profile by providing a Control Lyapunov Functional. In addition, the proof in the next chapter shows that the overshoot in (2) as well as the constants involved in the CLF (29) are independent of a particular reference trajectory. However, the convergence rate in (2) and (31) is governed by .
5 Proof of the Main Results
Before starting with the proof of the main theorem, we compute the closedloop dynamics in dependence of the delay variables and and . We use the control law (22)–(25), ODE (14) and definition (19) and conclude to
(32)  
(33)  
(34) 
Recall that and the mapping are input independent by virtue of (13) and (18). With this result, we are in the position to describe the closedloop dynamics of the PDE system (2)–(2) completely with the delay model:
Lemma 3 (controlled dynamics)
Proof of Lemma 3.
Assumptions (A)–(A) ensure the existence of a unique solution of (2)–(2) as elucidated in Lemma 3.2 of [12]. The solution of (15) belongs to , since it coincides with the solution of the delay differential equation
(35) 
with the same IC. On this basis, the function defined by (16) is of class for any , and any piecewise continuous input . Define the set
where is the finite (possibly empty) set where the derivative of is not defined or is not continuous.
The function is continuous and welldefined by (18) for all and every kernel defined by (9) under assumption (A).
As a result, the solution of the differential equations (32)–(34) is unique and exists locally for every . Since the righthand sides of (32)–(34) satisfy a linear growth condition, the solutions and exist for all .
The next step is to prove that the solutions of (32)–(34) can be utilized for parameterizing the PDE solution. For all and the relations
(36) 
hold true from (13)–(16) with . In addition, it is straightforward to show the following identities on using (14)–(16)
(37) 
In accordance to (2), (2) we get
(38)  
(39) 
which imply the relation in conjunction with (36), (37). At last, we evaluate the expression with (8), (9) and obtain (17), which completes the proof.
For proving Theorem 2, we need further lemmas which characterize the dynamics of the observer error and the scalar ODEstate . The proofs are given in Appendix A.
Lemma 4 (stability of observer dynamics)
Lemma 5 (stability of scalar dynamics)
At last, define the functionals
(45)  
(46) 
with the constants , , to be selected in the following.
Proof of Theorem 2. At first we analyze the dynamics of the function , see (18), governed by the IDE (15). To this end, we make use of results on the stability of IDEs in [11, 12]. More specifically, we use Corollary 4.6 in [12] as a main tool and preliminarily discuss necessary assumptions in order to apply it properly to the IDE subsystem governed by (15) with IC (13). The existence of a constant , satisfying
(47) 
is a direct consequence of assumption (B) with a constant . Next, consider the continuous functional defined by
(48) 
using the notation and notice the fact by virtue of the IC (13). Therefore, it follows from Corollary 4.6 in [12] that
(49) 
is a Lyapunov functional of the delay subsystem (13), (15) by means of the differential inequality for all with . In addition, Remark 4.7 [12] implies that the functional defined by
(50) 
is nonincreasing. As a result, it follows for all and . Using the fact for all , we conclude from (18):
(51) 
The second part of the current proof aims at deriving bounds for the analysis parameters and independent of a reference trajectory, such that is a CLF of the closed loop (2)–(2), (22)–(25). For this purpose, the following relations between the delay model (13)–(17) and the PDE problem (2)–(2) are a direct consequence of definition (45), (5):
(52)  
(53) 
Subsequently, we obtain the following inequality for the Diniderivative of (5) for any
(54) 
by identifying the nonincreasing functional (50) in the denominator of (5). Choose such that holds. To this end, we conclude that the estimate
(55) 
holds with the constant Beyond that, for all with we combine
(56) 
with the fact and get the differential inequality
(57) 
Moreover, (55) implies that the function is nonincreasing, so if we have we also have for every , because . We conclude that (57) holds for all with . Thus, (57) is valid for all .
With the results up to this point, we are in the position to establish the following inequality for the Diniderivative of (45) for all based on (42), (44), (55), (57), viz.