Optimal control of a multi-level dynamic model for biofuel production

Dynamic flux balance analysis of a bioreactor is based on the coupling between a dynamic problem, which models the evolution of biomass, feeding substrates and metabolites, and a linear program, which encodes the metabolic activity inside cells. We cast the problem in the language of optimal control and propose a hybrid formulation to model the full coupling between macroscopic and microscopic level. On a given location of the hybrid system we analyze necessary conditions given by the Pontryagin Maximum Principle and discuss the presence of singular arcs. For the multi-input case, under suitable assumptions, we prove that generically with respect to initial conditions optimal controls are bang-bang. For the single-input case the result is even stronger as we show that optimal controls are bang-bang.


1.
Introduction. Biofuels provide a concrete answer to the pressing need for renewable energy. The problem of increasing the efficiency and reducing the cost of biofuel production has been subject of intensive research [2,11]. A feasible source to obtain biofuels consists in using ethanol produced by cyanobacteria and microalgae [15,16]. This paper deals with the application of optimal control techniques to a general bioreactor for biofuel production.
Since the last decades, optimization of bioprocesses has been a line of research connecting optimal control to system biology, see [28] and references therein. Progress in plant genetic engineering has opened novel opportunities to use plants as bioreactors for safe and cost effective production of vaccine antigens. A review of methods and applications of plant, tissue and cell culture based expression strategies and their use as bioreactors for large scale production of pharmaceutically important proteins can be found in [27]. Exploiting bioreactors was also proved a fruitful method to deal with the problem of wastewater treatment in environmental engineering [22]. In this context a dynamic optimization problem arising frequently is the minimization of the time needed to reach a fixed target configuration for the bioreactor, see [25]. In this case, the trajectory evolves according to a system of nonlinear ODEs and may satisfy some boundary constraints as well. The right-hand side of the evolution equations involves not only control functions, i.e. parameters through which we can modify the dynamics, but also some unknown functions (e.g. biomass growth rate) of the evolving quantities themselves. Therefore, computing optimal controls through standard techniques (necessary conditions) is affected by the unknown functions' behavior. An interesting approach to deal with this issue using observability techniques was proposed in [5], where an application to bioreactors is also provided (see also [14]).
Other applications of optimal control techniques have been studied in the field of biological systems [1]. In the context of medical treatment, chemotherapy was used [19] as a dynamic control to optimize treatment scheduling in early stage HIVinfected cases. In that paper the authors use the effect of chemotherapy on viral production to maximize benefits in terms of T cell count and minimize the systemic cost of the treatement. Time dependent control strategies were also exploited in models to contain the emergence of drug-resistant strains of tubercolosis [18]. Here controls are represented by efforts in finding patients in which virus is only latent and in completing treatment for patients in which virus is already active. The objective function balances the effect of minimizing the cases of latent and infectious drug-resistant tubercolosis and minimizing the cost of implementing the control treatments.
In this paper we are concerned with the optimization of a bioprocess for biofuel production. Namely, we consider a model for the metabolic activity of a microorganism (e.g. E. Coli or Saccharomyces cerevisiae) in a fed-batch culture with different feeding substrates. The optimal control problem is to maximize the final amount of a certain side metabolite (e.g., ethanol) through different feed concentrations.
A steady-state approach to model cellular metabolism is Flux Balance Analysis (FBA, see [23]). The main assumption of FBA is that metabolic activity of cells is performed in such a way that the growth rate of cells is maximized. Since genomescale stoichiometric models for bacteria such as E. Coli are available, this translates in a linear program where the objective is cellular growth rate and constraints are given by metabolic reactions. Optimization is then performed by means of genetic manipulations on the bacteria (gene deletions and insertions). Other approaches based on flux balance analysis were proposed that take account of transcriptional and regulatory effects in [8], that couple the steady-state metabolic activity with a dynamic model [20,21] and that integrate both aspects [9]. Optimizing ethanol productivity in fed-batch cultures modeled through Dynamic Flux Balance Analysis [17] consists in using outputs of FBA (cellular growth rate and metabolite fluxes) to update at each time step the dynamics for evolving extracellular quantities. Within this framework, control can be performed at two levels: intracellular controls (genetic modifications), which are implemented by acting on constraints of the F-BA, and extracelullar controls (of dynamic nature), which are implemented through time-dependent feed concentrations. Genetic strategies were deeply investigated in [17] for in silico evolution of a yeast strain in glucose and xylose media to maximize ethanol productivity with constant feed concentrations. Besides genetic strategies, another parameter which is used as a control in [16,17] is the switching time of oxygen concentration. More precisely, the authors treat dissolved oxygen concentration as an independent variable, assuming it could be regulated by a feedback controller, and switch from aerobic to anaerobic growth to promote ethanol production during later stages of the batch. Then, they analyze sensitivity of ethanol productivity (for different genetic strains) to the aerobic-anaerobic switching time, while feed substrates' concentrations are kept constant.
In this paper, we focus on a unified model dealing with extracellular controls, i.e., glucose, xylose and oxygen feed concentrations, which we allow to be time dependent. Namely, we cast the problem of optimizing ethanol for the bioreactor in the language of optimal control, where concentrations of feeding substrates are modified by an external agent in order to enhance ethanol production. Even though we do not consider genetic manipulations, our model can integrate intracellular controls. Our purpose is to develop a general model featuring a full coupling between extracellular dynamics and intracellular metabolic activity: at each time instant solutions of ODEs provide fluxes to constraint the FBA, whereas FBA gives cellular growth rate and ethanol uptake appearing in the ODEs see Figure 1.  objective function and the constraints in FBA are linear (with respect to fluxes), outputs are piecewise linear. Dependence of fluxes on metabolite concentrations is modeled through Michaelis-Menten behavior, i.e., through rational functions. Therefore, FBA outputs appearing in dynamics are piecewise smooth functions of evolving quantities. The main idea is to implement the coupling between the intracellular and extracellular level through a hybrid control system where in each location the outputs of FBA are smooth. One of the advantages of a hybrid formulation is that it allows to take account of different timescales (for the general theory of hybrid control systems we refer the reader to [3,10]). One can either assume that information translates instantaneously from the micro level to the macro one (and viceversa) or one can implement delays, which are observed in experimental data. This is useful as there are configurations of the extracellular environment (such as saturation of glucose or oxygen) which may cause a major change in the metabolic pathway involved inside the cell with a consequent delay in the variation of outputs of the FBA.
As a first step toward a unified model, we focus on the analysis of optimal trajectories contained in a region where FBA ouptuts are smooth. In other words, we study what happens in each location of the hybrid model. More precisely, the location is characterized by a control-affine system of the typė where x ∈ R 9 is a vector containing substrate concentrations (biomass, glucose, xylose, ethanol, oxygen), u 1 , u 2 , u 3 are glucose, xylose and oxygen feed concentrations, y x are parameters coming from FBA, and F i , i = 0, . . . , 3 are smooth vector fields. The controlled vector fields F 1 , F 2 , F 3 are actually constant with respect to x. The dependence x → y x of FBA outputs (in the sequel, the biomass growth rate µ and the ethanol flux v e ) is smooth with respect to the state variable, as this is the case in each location of the hybrid model, where location switchings occur at nonsmoothness points of FBA outputs. Assumptions we require on the function x → y x , see hypothesis (H1) and (H2), are in accordance with the fact that it comes from a linear optimization problem, which we do not solve here. Under these conditions, the control system (1) is characterized by a sublinear growth. We consider an optimal control problem in Mayer form where we optimize the ratio between the total amount of ethanol (which is one component of the state vector x) and the total amount of available substrates at the final time (For the precise formulation we refer the reader to the next section.) We analyze necessary conditions for optimality given by Pontryagin Maximum Principle (PMP) [24]. Since the system is control-affine, the maximization condition of PMP allows to compute the control along an extremal trajectory as a feedback law, as long as the corresponding switching function (derivative of Hamiltonian with respect to control) has only isolated zeros. Otherwise, when the trajectory shows a singular arc, i.e., a time interval where a switching function vanishes identically, the main tool to find controls is exploiting conditions given by annihilation of higher order time-derivatives. As our interest is driven by applications, one would rather avoid singular arcs, which represent an obstacle to efficient and reliable numerical simulations. First, we study the general three-input case where glucose, xylose and oxygen feed concentrations are used as controls. The situation is rather intricate and some specific singular arcs may indeed be optimal. Under an additional assumption, which amounts to disregard the presence of a preferred substrate, we show that at most initial conditions optimal trajectories do not have singular arcs and are concatenation of arcs where all feed concentrations are either absent or maximal. We then consider the single-input case, where optimization is performed under the action of a single control. In this connection, we prove a stronger result stating that every optimal trajectory is a concatenation of arcs where the feed concentration is either absent or maximal. Let us point out that our analysis in the single-input case provides a theoretical justification of the choice of a piecewise constant profile (e.g. for oxygen concentration as in [17]), which is obviously the easiest to implement in practice.
The structure of the paper is the following. In Section 2 we recall the general model of a bioreactor and formulate the optimal control problem under consideration. In Section 3 we recall a classical first order necessary condition for optimal control problems. Then, we consider in Section 3.2 the single-input case and show that every optimal trajectory is bang-bang. In Section 3.1 we deal with the general case of three feeding substrates and analyze possible singular arcs. Finally, in Section 4 we compare our results with the ideas in [17].
2. Problem formulation. Bioreactors are processes where a living microorganism metabolizes some substrates, consequently grows and produces other metabolites. In [17] the authors consider in silico evolution of a yeast strain, Saccharomyces cerevisiae, which grows in a fed-batch culture with glucose and xylose and produces ethanol. The dynamic model can be applied to a general bioreactor.
We denote by V the total culture volume, which is assumed to grow linearly with respect to time with constant rate R. Biomass concentration is denoted by X. The total biomass in the culture evolves exponentially with a growth rate µ depending on substrates concentrations. Concentrations of feeding substrates, glucose G, xylose Z, and oxygen O, are characterized by an evolution which takes account of a feeding rate and a compensation term due to metabolism of the microorganism. Feed concentrations for substrates represent the control we perform on the system and are denoted by u 1 , u 2 , u 3 . The organisms metabolize glucose, xylose and oxygen with specific rates (or fluxes) v g , v z , v o depending on substrates concentrations. The produced metabolite under study is ethanol, with concentration E. We assume that the organism produces ethanol proportionally to the total biomass through a specific rate v e . Therefore, the control system we analyze is First of all, oxygen uptake kinetics follows Michaelis-Menten law whereas glucose uptake kinetics has an additional regulatory term to capture growth rate suppression due to high ethanol concentration, i.e., Xylose uptake kinetics has a similar form with another regulatory term to account for inhibited xylose metabolism in presence of the preferred substrate (glucose), Parameters v o max , v g max , v z max , k o , k g , k g ie , k z , k z ie are positive and constant. As for µ, v e , the model is based on the principle that the metabolic activity of the microorganism is performed so that biomass growth rate is maximized. This is equivalent to say that µ and v e are outputs of the optimization problem stated as In (3),v ∈ R n is the vector of fluxes (among which glucose, xylose, oxygen and ethanol fluxes) considered in the model; w ∈ [0, 1] n is the vector of weights which determines fluxes producing biomass; S ∈ M r×n (R) is the stoichiometric matrix which encodes the metabolic network inside the cell (involving r reactions and n metabolites); v max is an upper bound associated with the microorganism, (see [23] for an exhaustive treatment of metabolic networks in systems biology.) For our purposes, w, S, v max are considered as given parameters. They depend on the specific strain of microorganism used in the bioreactor and are usually determined through experimental data, see for instance supplementary data of [9] for E. Coli or [16] for Saccharomyces cerevisiae.
Equation (4) . This property is our main assumption below (see (H1) for the multi-input case and (H2) for the single-input case).
We fix the final time t f > 0 and an initial condition (G 0 , Z 0 , O 0 , E 0 , X 0 , V 0 ) and we seek to optimize the ratio between the total amount of ethanol at final time and the total amount of available substrates along a trajectory of (2). More precisely, we optimize the cost where the class of admissible controls is Values 0 and 1 for the controls correspond respectively to no feed or maximal feed concentration. We fix a positive thresholdū, which can be interpreted as the total amount of available substrates, and we restrict the maximization problem (5) to controls (u 1 , u 2 , u 3 ) ∈ U which satisfy the further constraint By construction,ū ∈ (0, 3t f ]. Adding the constraint above means that we exploit the totality of available substrates and oxygen. The problem shares a full coupling between a classical optimal control problem (2), (5), which describes the dynamics outside the cell, and a linear optimization problem (3), which models the metabolic activity inside the cell. In other words, at each time instant t, on the one hand, one needs to solve the linear program (3) to obtain µ, v e at time t to plug into (2); on the other hand, dynamics in (2) must be integrated to in order to get the tuple (G(t), Z(t), O(t), E(t)) that allows to determine constraints in (3).
Since in (3) both the objective function and the constraints are linear with respect tov, the outputs µ and v e are piecewise linear with respect tov. Therefore µ, v e are piecewise smooth as functions of (Z, G, O, E). In the sequel we assume that the outputs of (3) are smooth. This amounts to say that the system evolves in a domain where µ, v e are smooth as functions of (Z, G, O, E). Under this assumption, we study optimal trajectories for problem (2), (5).
Let us give a more compact formulation of the optimization problem above. Let us rename In order to get a cost of Mayer form in (5), i.e., a cost depending only on the final point of the trajectory, it is natural to add three new variables to (2) by setting so that the control system (2) reads where and v e (x 4 , x 5 , x 6 , x 7 , x 9 ), µ(x 4 , x 5 , x 6 , x 7 , x 9 ) are defined in the obvious way by (3), (4) using (6). The equation for x 9 can be trivially integrated. Indeed, the total volume of the culture plays the same role as time along the experiment. Here we prefer to keep the state variable x 9 so that the dynamics is autonomous.
take an initial point and consider the target

ROBERTA GHEZZI AND BENEDETTO PICCOLI
The optimal control problem (5) (with the additional constraint on the controls) is equivalent to where We end this section by stating a result which ensures the existence of optimal solutions to (9) with suitable initial conditions.
Assume v e , µ : R 5 → R are smooth, non negative and bounded from above. Then for every x 0 ∈ D, the optimal control problem (9) admits a solution.
In the proposition above we consider general functions µ and v e of 5 variables (as so is the case in problem (9)). Note that, whenever µ, v e are outputs of the optimization problem (3), they are non negative and bounded from above by construction, thanks to the constraints 0 ≤v j ≤ v max j in (3), so the assumptions of Proposition 1 are satisfied in practical situations.
Proof of Proposition 1. Let x 0 ∈ D. We are going to prove that every trajectory of the control system starting at x 0 is well-defined for all t > 0 and satisfies Indeed we have x 9 (t) = x 0 9 + Rt ≥ x 0 9 > 0 for every t. Using the assumptions on µ, v e , the dynamics is smooth onD and satisfies the sublinear growth condition for a certain constant C > 0. Let T be such that a trajectory starting at x 0 is well-defined on [0, T ]. Then x(t) ∈D for every t ∈ [0, T ]. Indeed, inequalities x i (t) ≥ 0 for i = 1, 2, 3, and x 9 (t) ≥ x 0 9 > 0 are trivially satisfied. Similarly, since the equation for x 8 is linear with respect to x 8 and µ(x 4 , x 5 , x 6 , x 7 , x 9 ) ≥ 0, we have x 8 (t) ≥ 0. By assumption, the ethanol flux satisfies v e (x 4 , x 5 , x 6 , x 7 , x 9 ) ≥ 0, whence x 7 (t) ≥ 0. We are left to prove that x i (t) ≥ 0 for i = 4, 5, 6. By assumption Eventually, trajectories starting at x 0 ∈ D ⊂D satisfy x(t) ∈D and, because of the growth condition (10), they are well-defined for every t > 0. Choosing , the trajectory associated to (u 1 , u 2 , u 3 ) reaches the target T at time t f . Hence for every x 0 ∈ D, a trajectory starting at x 0 and reaching the target always exists. The dynamics (7) is control-affine, whence the set of velocities {f (x, u) | u ∈ [0, 1] 3 } is convex for every x. Moreover, the target T is closed, and the cost function ψ is smooth. Therefore, the existence of an optimal solution is a consequence of classical results, see for instance [4, Theorem 5.1.1].
3. Analysis of extremal trajectories. In this section we consider first order necessary conditions for optimality to identify properties of solutions to (9). We first analyze a simplified case, where only one feeding substrate is consider and then turn to the case of three controls.

Consider a general optimal control problem in Mayer form
where The main tool to look for solutions to (11) is the Pontryagin Maximum Principle [24] which is a first order necessary condition for optimality. The Hamiltonian associated to (11) is H : where λ denotes the adjoint vector. Let u ∈ U be an admissible control whose corresponding trajectory x satisfies the terminal constraint the maximisation condition for almost every t ∈ [0, t f ], and the transversality condition The Pontryagin Maximum Principle, see for instance a version adapted to our problem in [4, Theorem 6.3.1], states that if (u, x) is a solution to (11), then x(·) is an extremal trajectory. Given an extremal trajectory, define the switching functions Since the Hamiltonian is control-affine, whenever t is such that ϕ i (t) = 0, the maximization condition (13) allows to identify u i (t) by u i (t) = 1+signϕi(t) 2 ∈ {0, 1}. We say that an extremal trajectory x(·) has a singular arc on [T 0 , , we say that the control u i is singular.
We say that a measurable control The family of controls (and trajectories) among which one optimizes in (11) is infinite dimensional. Thus, in order to solve (11) a fundamental step is to reduce the problem to a finite dimensional one. For instance, this can be done if, by means of necessary conditions, one deduces that every extremal trajectory is a finite concatenation of arcs where every control is smooth (e.g. bang-bang) and thus has a simple form. To do so, in the control-affine case, it suffices to prove that every extremal trajectory corresponds to a finite concatenation of bang-bang controls. In general this is hard to prove, mainly because either the concatenation of simple arcs may be infinite or because there may be singular arcs, where the control is not uniquely determined.
As concerns the first difficulty, it may indeed happen that an optimal trajectory features an accumulation of arcs where controls are constant. This is the well-known Fuller phenomenon, see [12]. In this work we do not deal with such pathological phenomena. We refer the reader to [6] for new ideas overcoming numerical difficulties caused by Fuller phenomenon.
As concerns the second issue, in the worst case, one can still expect that only a "few" extremal trajectories are "bad", i.e., they feature a singular arc. Roughly speaking, this amounts to prove that extremal trajectories corresponding to a singular control are confined in some subsets of positive codimension. If so, one infers that optimal trajectories starting at most initial conditions do not have singular arcs and correspond to bang-bang controls.
With regard to the presence of optimal singular trajectories for control-affine systems, a genericity result has been shown in [7]. In that paper, the authors prove that, for a quadratic cost, generically with respect to the control system there do not exist nontrivial optimal singular trajectories [7, Corollary 2.9]. Results of this type do not apply to our problem because system (7) is not generic, F i being constant vector fields. However, modeling µ, v e as affine functions of v g , v z , v o , in the threeinput case, we are able to single out only two (out of seven) possible singular arcs and generically exclude all the other ones.

3.1.
Multi-input case. The Hamiltonian associated to the optimal control problem (9) is H : where λ ∈ R 9 denotes the adjoint vector. We model the behavior of µ, v e assuming the following property, see also Remark 1.
One possibility for a singular arc is that u 1 (which corresponds to glucose feed concentration) is singular, u 2 is bang-bang (u 3 may or may not be singular). Otherwise, all other singular arcs are confined in some hypersurface of the type {x ∈ R 9 | P (x) = 0}. Therefore, for most initial conditions, solutions to (9) correspond to controls u such that u 2 is a concatenation of bang-bang arcs and u 1 , u 3 are concatenations of bang-bang or singular arcs.
When we do not account for a preferred substrate and assume v z does not depend on glucose concentration x 4 , then the result states that for most initial conditions optimal trajectories for (9) correspond to concatenations of bang-bang arcs. In other words, generically with respect to initial conditions, singular arcs do not occur.
Proof. The matrix N (x) is given by where 0 l×m denotes the null matrix of dimensions l × m, x 8 ∂vg ∂x7 x 8 ∂vz ∂x7 Let t ∈ [0, t f ] and λ ∈ ker N (x(t)). Then Therefore, a vector λ ∈ R 9 belongs to ker N (x(t)) if and only if it satisfies (20) if µ =v = 0, or it satisfies (21), otherwise.
The first three equations in (24) are trivially integrated and give The proof consists of 8 steps.
Step 1 analyzes equilibria of (24) and ensures that any solution of (24) that is constant on some interval [T 0 , Step 2 is concerned with singular arcs where all the switching functions vanish identically and implies that at least one among the controls is bang-bang. Steps 3,4,5 deal with singular arcs where two switching functions vanish identically and one does not. Steps 6,7,8 deal with singular arcs where one switching function vanishes identically and the other two do not.
Step 1. The adjoint system (24) can be written asλ = N (x)λ, where N (y) is the matrix given in Lemma 3.2. Reasoning as in the proof of Proposition 1, we infer that x(·) is uniformly bounded on [0, t f ] (growth in the control system is sublinear, see (10)). As a consequence, the matrix function t → N (x(t)) is uniformly bounded on [0, t f ]. Hence we have existence and global uniqueness for solutions to the adjoint system (24). Let us now analyze equilibria for (24). By Lemma 3.2, the kernel of N (x(t)) does not depend on t. Hence, by global uniqueness, if a solution λ(·) to (24) is constant on some interval [T 0 , T 1 ] with T 1 > T 0 then λ(·) is constant on the whole interval [0, t f ] and belongs to the subspace defined in (20) or (21). For simplicity, in the sequel we set A = λ 4 − b 1 λ 7 − a 1 λ 8 , B = λ 5 − b 2 λ 7 − a 2 λ 8 , C = λ 6 + b 3 λ 7 − a 3 λ 8 , D =vλ 7 +μλ 8 , see system (45). As it happens, λ is en equilibrium if and only if A = B = C = D = 0, and (24) reads Step 2.
Step 5. Let [T 0 , T 1 ] ⊂ [0, t f ] be a singular arc where ϕ 2 is nonzero and ϕ 1 , ϕ 3 vanish identically. Then u 2 (t) = 1+signϕ2(t) 2 ∈ {0, 1} for almost every t. Moreover, Imposing that λ 4 , λ 6 are constant, by (26) we deduce C = 0 and ∂v g ∂x Differentiating condition C = 0 with respect to time and taking (29), (26) into account we obtain another equation in (A, B, D), namely In order to carry out the same argument as in steps 3 and 4, one should find a third linear and homogeneous condition in A, B, D. The problem is that, if one differentiates (29) or (30) with respect to time, the obtained conditions involvė x 4 , which depends on the control u 1 . Since we are on an arc where ϕ 1 ≡ 0, this control is singular and thus unknown. Therefore, singular arcs of this type cannot be excluded in general.
On the contrary, we can do better under the additional assumption that v z is given as in (19). In this case, since v z does not depend in x 4 , condition (29) is equivalent to A = 0. Differentiating conditions A = 0, C = 0 with respect to time we obtain a linear homogeneous system in B, D, namely The determinant of the system above is ( , only depends on the value of (x 4 , x 7 , x 8 , x 9 ) and does not depend explicitly on t. Step 6. Let [T 0 , T 1 ] ⊂ [0, t f ] be a singular arc where ϕ 1 , ϕ 3 are nonzero and ϕ 2 vanishes identically. Then u i (t) = 1+signϕi(t) 2 ∈ {0, 1}, i = 1, 3 for almost every t. Moreover, From (26) we infer that B = 0. Differentiating this condition with respect to time we obtain Differentiating twice (32) with respect to time gives rise to other two linear homogeneous equations in (A, C, D) whose coefficients depend only on (x 4 , x 6 , x 7 , x 8 , x 9 ) and (ẋ 4 ,ẋ 6 ,ẋ 7 ,ẋ 8 ,ẋ 9 ). Since the controls u 1 , u 3 are constant on [T 0 , T 1 ], these coefficients do not depend explicitly on t. Considering the obtained linear homogeneous system in (A, C, D) there are two possibilities. Either there exists t ∈ [T 0 , T 1 ] such that the determinant of the system does not vanish or the determinant vanishes identically. In the first case we deduce that λ is an equilibrium of (24) and we get a contradiction with the transversality condition. In the second case, there exists a polynomial function P : R 9 → R depending only on (x 4 , x 6 , x 7 , x 8 , x 9 ) such that the trajectory is contained in {x ∈ R 9 | P (x) = 0}.
Since the control u 1 is singular and thus unknown, we are not able to exclude singular arcs of this type in general.
Under the additional assumption that v z is as in (19), the situation simplifies. More precisely, λ 4 constant implies A = 0. The system is linear and homogeneous in (B, C, D) and its coefficients depend only on the value of (x 5 , x 6 , x 7 , x 8 , x 9 ) and of (ẋ 5 ,ẋ 6 ,ẋ 7 ,ẋ 8 ,ẋ 9 ). Recalling that the controls u 2 , u 3 are constant on [T 0 , T 1 ] we can conclude as in step 7 and prove that either such a singular arc does not occur or the trajectory is contained in the hypersurface defined by {x ∈ R 9 | P (x) = 0}, where P is a polynomial function depending only on (x 5 , x 6 , x 7 , x 8 , x 9 ). Remark 2. Concerning singular extremals having two singular controls and a bangbang control there are three possibilities, that we consider in steps 3, 4, 5 of the proof above. In the first and second cases (steps 3, 4 respectively), the only possibility for a singular extremal to exist is that glucose, respectively oxygen, concentration is identically zero since the beginning of the experiment. In view of applications, these situations are rather meaningless, being that glucose is the preferred substrate for biomass growth and aerobic conditions enhance biomass growth. In the last case (step 5) under the simplified assumption (19), the only possibility is that xylose concentration is identically zero since the beginning of the experiment. This represents the only reasonable possibility from the point of view of applications.

3.2.
Single-input case. In this section we study a simpler version of the optimal control problem (9) by considering anaerobic growth and removing one feeding substrate, e.g., xylose. Within this simplified framework, it is possible to go further in the analysis carried out as in Theorem 3.1 and show a stronger result, see Proposition 2 below.
Under this hypothesis, µ and v e are smooth, non negative and bounded from above on the set {y ∈ R 5 | y 5 > 0, y i ≥ 0}, so that existence of solutions to (38) can be proved as in Proposition 1. See Remark 1 for a justification of the affine behavior of µ, v e with respect to v g .
Proposition 2. Let µ, v e satisfy assumption (H2). Then any optimal control for (38) is bang-bang on [0, t f ] and has no singular arcs.
Proof of Proposition 2. By the Pontryagin Maximum Principle, any optimal control for (38) gives rise to an extremal trajectory. Therefore, it suffices to show that given an extremal trajectory y(·) with adjoint vector η(·) satisfying (40), (41), (42), then the switching function does not vanish identically on an interval of positive length.
The proof is split into two steps. First, we show that given an extremal trajectory y(·) then any solution η(·) of (40) which is constant on an interval [T 0 , Second, we prove that given an extremal trajectory y(·) and a corresponding covector satisfying (40), (41), (42), then the presence of a singular arc implies that the covector is constant on [0, t f ] and we get a contradiction imposing the transversality condition.
Equilibria of (40). Let u ∈ U be an admissible control whose corresponding trajectory y satisfies the terminal constraint y(t f ) ∈ T and is extremal. Let η(·) be a solution of (40). Once y(·) is given, system (40) can be written asη = M (y)η where M (y) is the matrix given in Lemma 3.3. Reasoning as in the proof of Proposition 1, we infer that y(·) is uniformly bounded on [0, t f ] (growth in the control system (35) is sublinear). As a consequence, the matrix function t → M (y(t)) is uniformly bounded on [0, t f ]. Hence we have existence and global uniqueness for solutions to the adjoint system (40).
Let us now analyze equilibria for (40). By Lemma 3.3, the kernel of M (y(t)) does not depend on t. Hence, by global uniqueness, if a solution η(·) to (40) is constant on some interval [T 0 , T 1 ] with T 1 > T 0 then η(·) is constant on the whole interval [0, t f ] and belongs to the subspace defined in (45) or (45).
Remark 3. When we consider the single-input system obtained by removing oxygen and glucose, the conclusion of Proposition 2 holds. This follows directly by the fact that, without oxygen and glucose, the equation for v z given in (8) becomes precisely (36) and the analysis is the same as the one above. Furthermore, when we remove glucose and xylose and keep oxygen as the only control, the analysis is even simpler, since v o in (8) does not depend on ethanol concentration. As a consequence, with the same technique, the conclusion of Proposition 2 holds true for the single-input case where the only control is oxygen feed concentration.

Conclusions.
In [17], the authors perform two series of in silico experiments for Saccharomyces cerevisiae. The first one considers glucose media, whereas the second one deals with glucose and xylose media. In both series, feed concentrations (which in our framework are represented by u 1 , u 2 ) are assumed constant and ethanol productivity is maximized among trajectories characterized by a dissolved oxygen concentration of 50% on [0, t s ] and of 0% from [t s , t f ] and the switching time t s is the only control parameter. In Section 3.1 we consider a the more general case where glucose, xylose, oxygen feed concentrations are treated as controls and µ, v e are affine functions of v g , v z , v o . In this unified framework, we prove that, at least when there is no preferred substrate (among glucose and xylose), then optimal trajectories exiting from most initial conditions are characterized bang-bang feed concentrations. In particular, in this simplified context, the only case of interest (in view of applications) for a singular extremal having two singular controls is when glucose and oxygen are singular and xylose concentration is zero since the beginning of the experiment. It would be of interest to check optimality of these particular singular extremals via in silico experiments.
Furthermore, in Section 3.2 we consider the single-input case where the feed concentration of a single substrate (glucose, xylose or oxygen) is treated as a control. In this framework our analysis provides a stronger result: optimal trajectories are bang-bang (see Remark 3 for the discussion concerning different susbstrates). This provides the theoretical background needed to legitimate the choice in [17] of piecewise constant dissolved oxygen concentration: optimizing ethanol productivity among all possible input profiles for the oxygen feed concentration is equivalent to optimizing among profiles that are bang-bang.
The analysis of singular extremals provided here constitutes the building block for the development of the hybrid model where parameters µ and v e are only piecewise smooth. To go further, the starting point is to study necessary conditions for optimality, which, for this generalized context, can be found in [13,26].
Finally, due to the high dimension of our problem, implementing efficient numerical simulations is rather heavy and therefore a numerical study is out of the scope of the present paper.