Yield Trajectory Tracking for HyperbolicAge-Structured Population Systems\thanksreffootnoteinfo

# Yield Trajectory Tracking for Hyperbolic Age-Structured Population Systems\thanksreffootnoteinfo

[    [    [
###### Abstract

For population systems modeled by age-structured hyperbolic partial differential equations (PDEs) that are bilinear in the input and evolve with a positive-valued infinite-dimensional state, global stabilization of constant yield set points was achieved in prior work. Seasonal demands in biotechnological production processes give rise to time-varying yield references. For the proposed control objective aiming at a global attractivity of desired yield trajectories, multiple non-standard features have to be considered: a non-local 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 first-order hyperbolic PDEs and integral delay equations leading to a decoupling of the input-dependent dynamics and the infinite-dimensional internal one. Furthermore, the dynamic control structure does not necessitate exact knowledge of the model parameters or online measurements of the age-profile. With a Galerkin-based numerical simulation scheme using the key ideas of the Karhunen-Loève-decomposition, 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 92093-0411, USA

Key words:  first-order hyperbolic PDE, chemostat, tracking control, time-delay systems, input constraints, Galerkin methods.

11footnotetext: Email: kevin.schmidt@isys.uni-stuttgart.de (Kevin Schmidt), iasonkar@central.ntua.gr (Iasson Karafyllis), krstic@ucsd.edu (Miroslav Krstic).

## 1 Introduction

We design an asymptotically tracking control for age-structured 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 Galerkin-methods, which ensures accurate asymptotic properties.
Motivation. In the context of mathematical biology and demography age-structured continuous-time 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 age-structured 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 McKendrick-vonFoerster 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 age-distribution. As a result, it is possible to represent all products which are proportional to the overall population as well as possibly age-dependent 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.
Time-varying 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 first-order hyperbolic PDEs to delay models in the sense of integral delay equations (IDEs, see [11]). More specifically, we decompose the PDE problem to an input-dependent finite-dimensional 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 two-degrees-of-freedom (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 output-feedback controller does not demand online measurements of the population’s entire age-distribution. 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 closed-loop system in Section 6 with a Galerkin-based simulation scheme, which conserves important system properties even at low orders and enables independent age and time discretizations.
Notation.

• The set denotes all positive-valued real numbers, all non-negative 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 Dini-derivative .

• 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 Age-Structured Population Models

Consider the McKendrick-vonFoester PDE (2) valid in the age-time domain

 ∂x(a,t)∂t+∂x(a,t)∂a =−[μ(a)+D(t)]x(a,t) (1) x(0,t) =A∫0k(a)x(a,t)da (2) x(a,0) =x0(a) (3)

which describes the evolution of the population density as a part of an initial-boundary 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 age-dependent mortality rate and the dilution rate which is the control input. In particular, the non-local boundary condition (BC) (2) is valid for and models the production of new-born 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 age-domain  at . In addition, the output is defined by the equation

 y(t)=A∫0p(a)x(a,t)da, (4)

which takes the possibly age-dependent 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 single-input-single-output 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]):

1. The parameters functions are restricted to and , where .

1. The control  takes values in where .

2. The initial condition (IC) (2) is compatible with (2), i.e. .

## 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 2-DOF control structure satisfying assumption (A) a priori

 D(y,yref)=sat(DFF(yref)+DFB(y,yref)) (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

 limt→∞|y(t)−yref(t)|=0     ∀x0∈X. (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

 1=A∫0k(a)e−D∗a−∫a0μ(α)dαda=:A∫0~k(a)da. (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

 x∗(a)=e−D∗a−∫a0μ(α)dα∫A0p(s)e−D∗s−∫s0μ(α)dαds. (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 :

 g(a) =x∗(a)⋅p(a), for all a∈[0,A] (9) xref(a,t) =x∗(a)⋅yref(t), for all a∈[0,A],t≥0. (10)

In the following proposition we refer to the open-loop 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

 Π(f):=⟨π,f⟩⟨π,x∗⟩=∫A0π(a)f(a)\emphda∫A0π(a)x∗(a)\emphda (11)

where is the continuous function

 π(a):=A∫ak(s)eD∗(a−s)+∫asμ(α)\emphdα\emphds. (12)

Given an IC for the PDE system (2)–(2), we define

 η0=lnΠ(x0)y\emphref(0),ψ0(a)=x0(a)x∗(a)Π(x0)−1,a∈[0,A]. (13)

For every piecewise continuous input , consider the solution of the delay model

 ˙η(t) =D∗−˙y\emphref(t)y\emphref(t)−D(t) (14) ψ(t) =A∫0~k(a)ψ(t−a)\emphda. (15)

with the ICs and for . Define the functions for ,

 Φx(a,t) =x\emphref(a,t)eη(t)[1+ψ(t−a)] (16) Φy(t) =y\emphref(t)eη(t)⎡⎢⎣1+A∫0g(a)ψ(t−a)\emphda⎤⎥⎦. (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 input-dependent ODE-state  and an infinite-dimensional 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

 δ(t):=ln⎛⎜⎝1+A∫0g(a)ψ(t−a)da⎞⎟⎠ (18)

we rewrite (17) and get following result for the logarithmic tracking error as the basis for the control design

 ln(y(t)yref(t))=η(t)+δ(t). (19)

## 4 Nonlinear Controller Design

For the design of an appropriate controller, we state a second set of assumptions:

1. The equilibrium dilution rate  defined by (7) satisfies .

2. Assume that belongs to the class of valid reference trajectories , which is defined as all positive-valued and continuously differentiable functions with the additional property

 D∗−Dmax
3. We assume the existence of a constant which fulfills the inequality with defined by (7)

 A∫0∣∣ ∣∣~k(a)−λ∫Aa~k(s)% ds∫A0s~k(s)ds∣∣ ∣∣da<1. (21)
4. 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 age-domain , 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 2-DOF structure in (3) with the control gain  and the observer gains for

 DFF(t) =−˙yref(t)yref(t) (22) DFB(t) =z2(t)+γlny(t)yref(t) (23) ˙z(t) =[−l11−l20]Lz(t)+⎡⎢ ⎢⎣l1lny(t)yref(t)+DFF(t)−D(t)l2lny(t)yref(t)⎤⎥ ⎥⎦ (24) D(t) =sat(−˙yref(t)yref(t)+z2(t)+γlny(t)yref(t)). (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 observer-ODE (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

 e(t)=z(t)−[lnΠ(x[t])yref(t), D∗]⊺=z(t)−[η(t), D∗]⊺ (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 :

 ∥∥∥lnx(a,t)x\emphref(a,t)∥∥∥∞+|e(t)|≤ ×\emphe−L(y\emphref)4t (27)

Moreover, for any constants satisfying

 (2+l1p1−2l2p2)2<8l1p1−4l2p21,    p21<4p2, (28)

there exists a sufficiently small constant , sufficiently large constants , such that for every the functional defined by

 V(z,x,t)=(lnΠ(x)y\emphref(t))2+α1√Q(z,x,t)+α2Q(z,x,t) (29)

with

 Q(z,x,t) =⎡⎣z1−lnΠ(x)y\emphref(t)z2−D∗⎤⎦⊺[1−p12−p12p2]⎡⎣z1−lnΠ(x)y\emphref(t)z2−D∗⎤⎦ +M2⎛⎜ ⎜ ⎜ ⎜⎝∥∥∥e−σa(x(a)x\emphref(a,t)−Π(x)y\emphref(t))∥∥∥∞min{Π(x)y\emphref(t),mina∈[0,A]x(a)x\emphref(a,t)}⎞⎟ ⎟ ⎟ ⎟⎠2 (30)

is a Lyapunov functional for the closed loop system (2)–(2), (22)–(25). More specifically, the differential inequality

 ˙V+(z(t),x[t],t)≤−L(y\emphref)⋅V(z(t),x[t],t)1+√V(z(t),x[t],t) (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 closed-loop 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 closed-loop dynamics in dependence of the delay variables and and . We use the control law (22)–(25), ODE (14) and definition (19) and conclude to

 ˙η(t)= D∗  −˙yref(t)yref(t)−sat(z2(t)−˙yref(t)yref(t)+γη(t)+γδ(t)) (32) ˙z1(t)= z2(t)−˙yref(t)yref(t)−%sat(z2(t)−˙yref(t)yref(t)+γη(t)+γδ(t)) −l1[z1(t)−η(t)−δ(t)] (33) ˙z2(t)= −l2[z1(t)−η(t)−δ(t)]. (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 closed-loop dynamics of the PDE system (2)–(2) completely with the delay model:

###### Lemma 3 (controlled dynamics)

Consider the controlled delay model governed by (15), (32)–(34) with the ICs (13), and . The unique solution , and exists for all , and parameterizes the solution of the PDE model (2)–(2), (22)–(25) by virtue of (16)–(17).

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

 dψ(t)dt=~k(0)ψ(t)−~k(A)ψ(t−A)+A∫0~k′(a)ψ(t−a)da (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 well-defined 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 right-hand 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

 η(t)=lnΠ(Φx[t])yref(t),  ψ(t−a)=Φx(a,t)x∗(a)Π(Φx[t])−1 (36)

hold true from (13)–(16) with . In addition, it is straightforward to show the following identities on using (14)–(16)

 ∂Φx(a,t)∂t+∂Φx(a,t)∂a=−[μ(a)+D(t)]Φx(a,t). (37)

In accordance to (2), (2) we get

 Φx(0,t) =⟨k,Φx[t]⟩ for all t≥0, (38) Φx(a,0) =x0(a) for all a∈[0,A], (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 ODE-state . The proofs are given in Appendix A.

###### Lemma 4 (stability of observer dynamics)

Consider the observer error  defined by (26) in context of the closed loop system (15), (32)–(34) under the assumptions of Theorem 2. Define the positive quadratic form

 J(e)=e21−p1e1e2+p2e22=:e⊺Pe>0 (40)

Then, there exist positive constants such that the following inequalities holds for all

 K1∣∣e(t)∣∣2 ≤J(e(t))≤K2∣∣e(t)∣∣2, (41) \emphdJ(e(t))\emphdt ≤−2β1J(e(t))+β2|δ(t)|2. (42)
###### Lemma 5 (stability of scalar dynamics)

Consider the dynamics of the state  governed by (13), (32)–(34) under the assumptions of Theorem 2. Define the positive constant and the functional

 μ1(y\emphref)= min(2,γ)⋅min{1,D∗−Dmin−supτ≥0˙y\emphref(τ)y\emphref(τ), Dmax−D∗+infτ≥0˙y\emphref(τ)y\emphref(τ)} (43)

depending on the reference trajectory. Then, for every  the following estimate holds with for all

 dη2(t)dt≤−μ1(y\emphref)η2(t)1+√η2(t)+μ2∣∣∣e2(t)+γδ(t)∣∣∣. (44)

At last, define the functionals

 ˜V(e,η,ψ)= η2+α1√˜Q(e,ψ)+α2˜Q(e,ψ) (45) ˜Q(e,ψ)= M2⎛⎜⎝∥∥e−σaψ(−a)∥∥∞1+min{0,mina∈[0,A]ψ(−a)}⎞⎟⎠2 +e21−p1e1e2+p2e22. (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

 A∫0eσa∣∣ ∣∣~k(a)−λ∫Aa~k(s)ds∫A0s~k(s)ds∣∣ ∣∣da<1 (47)

is a direct consequence of assumption (B) with a constant . Next, consider the continuous functional defined by

 P(ψt)=1∫A0s~k(s)dsA∫0A∫aψt(−a)~k(s)dsda (48)

using the notation and notice the fact by virtue of the IC (13). Therefore, it follows from Corollary 4.6 in [12] that

 W(ψt)=∥∥e−σaψt(−a)∥∥∞ (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

 C(ψt)=1+min(0,mina∈[0,A]ψt(−a)) (50)

is non-increasing. As a result, it follows for all and . Using the fact for all , we conclude from (18):

 |δ(t)|≤eσA∥∥e−aσψt(−a)∥∥∞1+min{0,mina∈[0,A]ψt(−a)} for all t≥0. (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):

 ˜V(e(t),η(t),ψt) =V(e(t)+[η(t),D∗]⊺,Φx(a,t),t) (52) ˜Q(e(t),ψt) =Q(e(t)+[η(t),D∗]⊺,Φx(a,t),t) (53)

Subsequently, we obtain the following inequality for the Dini-derivative of (5) for any

 ˙˜Q+(e(t),ψt) ≤−2β1J(e(t))−(σM−β2e2σA) ×⎛⎜⎝∥∥e−σaψt(a)∥∥∞1+min{0,mina∈[0,A]ψt(a)}⎞⎟⎠2. (54)

by identifying the non-increasing functional (50) in the denominator of (5). Choose  such that holds. To this end, we conclude that the estimate

 ˙˜Q+(e(t),ψt)≤−2β˜Q(e(t),ψt) (55)

holds with the constant Beyond that, for all with we combine

 ˜Q\nicefrac12(e,ψ)∣∣t+h−˜Q\nicefrac12(e,ψ)∣∣th≤˜Q(e,ψ)∣∣t+h−˜Q(e,ψ)∣∣th ×(˜Q\nicefrac12(e,ψ)∣∣t+h+˜Q\nicefrac12(e,ψ)∣∣t)−1 (56)

with the fact and get the differential inequality

 ˙√˜Q(e(t),ψt)+≤−β√˜Q(e(t),ψt). (57)

Moreover, (55) implies that the function is non-increasing, 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 Dini-derivative of (45) for all based on (42), (44), (55), (57), viz.

 ˙˜V+(e(t),η(t),ψt)≤ −μ1(yref)η2(t