Cette thèse est divisée en deux parties principales. Premièrement, l'objectif est d'étudier méthodes appropriées au contrôle et l'estimation des bio-procédés complexes. Deuxièmement, le problème du contrôle des systèmes incertains avec contraintes (sur les états et sur la commande) est résolu en fusionnant des estimateurs à intervalles avec la commande prédictive (MPC).

Pour la première partie, les bio-procédés concernés décrivent la croissance de populations microbiennes en culture continue. Deux scénarios différents sont envisagés : la compétition entre deux espèces pour un seul nutrient, et le cas d'un consortium microbien producteurnettoyeur, où deux souches de microbes effectuent des tâches en mutualisme.

Ces systèmes biologiques sont très complexes en raison de leurs non-linéarités, de leur nature contrainte (en raison de leur signification physique, les états sont toujours non-négatifs et caractérisent ces systèmes comme positifs) et sur la commande (les actionneurs sont physiquement limités à des valeurs non négatives), et en raison de l'incertitude inhérente aux processus biologiques.

Tout d'abord, profitant de leurs similitudes structurelles, l'observabilité de tels systèmes est étudiée, suivie par la conception d'estimateurs d'état appropriés utilisant des mesures réalistes. Ensuite, les architectures de contrôle proposées visent à garantir la coexistence de deux espèces (ou souches) et deux scénarios distincts sont considérés : dans le premier, une architecture de contrôle robuste est conçue pour le modèle de compétition, en utilisant un ensemble de différentes lois de commande discontinues et qui prennent explicitement en compte les incertitudes. Puis, dans un second temps, la conception des lois de commande pour le consortium microbien producteur-nettoyeur est présentée en considérant le cas nominal.

Pour la deuxième partie de la thèse, l'idée principale est d'aborder le problème de la commande robuste pour les systèmes contraints en incorporant des estimateurs à intervalles (un observateur et un prédicteur) dans l'algorithme classique du MPC. Ces estimateurs, sous certaines conditions de non-négativité des erreurs d'estimation, fournissent des enveloppes qui contiennent toutes les trajectoires possibles du système incertain. En couplant cette caractéristique avec un problème d'optimisation, les algorithmes proposés offrent faisabilité récursive, ainsi que la satisfaction robuste des contraintes, à une faible complexité de calcul et avec une mise en oeuvre très simple. Deux classes générales de systèmes sont concernées dans ce cadre : les systèmes linéaires invariants dans le temps et les systèmes linéaires à paramètres variables. Pour chaque cas, les conditions de conception pour les estimateurs d'intervalle, ainsi que pour leur commande à retour d'état, sont données sous la forme d'inégalités matricielles linéaires (résolues hors ligne).

These biological systems are highly complex due to their nonlinearities, their constrained nature in both states (due to their physical meaning, which characterizes such systems as positive) and control inputs (the actuators are physically restrained to non-negative values), and due to the inherent uncertainty of biological processes.

First, profiting of their structural similarities, the observability of such systems is studied, proceeded by the design of proper state estimators using realistic measurements. Then, the control solutions aim at guaranteeing coexistence of all species and is then divided in two scenarios: in the first one, a robust control architecture is designed for the competition model, using a set of different discontinuous control laws that take the uncertainties explicitly into account. Then, in a second moment, the design of control laws for the microbial consortium is presented considering the nominal case.

For the second part of the thesis, the main idea is to tackle the robust constrained control problem by incorporating interval estimators (an observer and a predictor) into the classical MPC algorithm. These estimators, under some conditions on non-negativity of the estimation errors, provide envelopes that contain all possible trajectories of the uncertain system.

Casting this feature into an optimization problem, the proposed algorithms offer robust constraint satisfaction and recursive feasibility at a low computational complexity and ease of implementation.

Two general classes of systems are concerned in the framework: Linear Time Invariant and Linear Parameter Varying systems. For each case, the conditions of design for the interval estimators, as well as for their feedback control, are given in the form of (offline) Linear Matrix Inequalities. • The absolute value of an element is denoted | • |, while the Euclidean norm is denoted ∥ • ∥.

List of Tables

• The vector V = vec(v 1 , . . . , v n ) ∈ R ν represents the concatenation of vectors v 1 , . . . , v n and ν = n i=1 ν i , where ν i is the dimension of each vector v i . The matrix M = diag(m 1 , . . . , m n ) is a diagonal matrix whose elements are M ii = m i , where m i are vectors or matrices of appropriate dimension. The transpose of the matrix M is denoted M ⊤ . For a symmetric matrix M , the symmetric entry (i.e., M i,j = M j,i , for i ̸ = j) is denoted by ⋆. The determinant of a matrix M is denoted det(M ).

• A matrix M is said to be non-negative if all of its elements are non-negative. A matrix M is said to be Schur stable if all of its eigenvalues have absolute value less than one.

The identity matrix is denoted by I n ∈ R n×n .

• Let x 1 , x 2 ∈ R n be two vectors and A 1 , A 2 ∈ R n×n be two matrices, then the relations x 1 ≤ x 2 and A 1 ≤ A 2 are to be understood component-wise. For a matrix A we define A + = max{0, A}, A -= A + -A (also understood in a component-wise sense, similarly for vectors), and also denote the matrix of absolute values of all elements by |A| = A + + A -. For a symmetric matrix A ∈ R n×n the relation A ≺ 0 (resp. A ⪰ 0) means that A ∈ R n×n is negative (resp. positive semi-) definite. If A is diagonal, then A > 0 is equivalent to A ≻ 0.

Continuous-time systems

• A continuous function α : R + → R + belongs to the class K if α(0) = 0 and the function is strictly increasing. A function β : R + ×R + → R + belongs to the class KL if β(•, t) ∈ K for each fixed t ∈ R + and β(s, •) is decreasing and lim t→+∞ β(s, t) = 0 for each fixed s ∈ R + . A function β : R + × R + → R + belongs to the class GKL if β(s, 0) ∈ K, β(s, •) is decreasing and for each s ∈ R + there is T s ∈ R + such that β(s, t) = 0 for all t ≥ T s ;

• For a Lebesgue measurable and essentially bounded function x : R → R n , denote

Background and motivation

Biotechnology has applications in several societal domains, such as medicine, agriculture, environmental sciences and industry. Indeed, many bioprocesses -processes in which a living organism (such as bacteria, yeast and microalgae), or their components, participate in the reactions -have been used by humans to obtain special products: ranging from brewing and fermentation to the development of new medicines.

Eventually, the interest in having biological systems with different capabilities or that do not occur in nature gave rise to synthetic biology [START_REF] Fu | Systems Biology and Synthetic Biology[END_REF]. This multidisciplinary and rather new field applies engineering approaches to design (or re-design existing) biological systems [Schmidt, 2012] and leverage new advanced technologies.

A recent trend in modern biotechnology is the concept of microbial consortia (or cocultures), i.e., the association of different species in a community (or the same species, but different strains 1 ) [Rosero- Chasoy et al., 2021]. The objective of such an association is that it might allow more complex tasks to accomplished, or even have a better performance than possible by a single organism [START_REF] Said | Synthetic microbial ecology: Engineering habitats for modular consortia[END_REF]. Indeed, many applications employing microbial consortia, such as bioprocessing [START_REF] Peng | Microbial communities for bioprocessing: lessons learned from nature[END_REF], biorefinery and the production of high-valuable compounds [START_REF] Rosero-Chasoy | Microbial co-culturing strategies for the production high value compounds, a reliable framework towards sustainable biorefinery implementation -an overview[END_REF], biosynthesis and biodegradation [Che and Men, 2019], have shown to be promising in industrial, medical, and environmental sciences.

In this light, synthetic biology and systems biology offer interesting tools regarding bioengineered synthetic microbial communities. For instance, it is possible to establish communication between organisms, "program" reactions to exogenous molecules or stimuli (such as light) to induce specific behaviours or establish syntrophic interactions in a network of microorganisms [START_REF] Mccarty | Synthetic biology tools to engineer microbial communities forbiotechnology[END_REF].

Over the last years, many researchers have devoted their attention to potential applications of control theory in synthetic biology [START_REF] Hsiao | Control theory for synthetic biology: Recent advances in system characterization, control design, and controller implementation for synthetic biology[END_REF]. Indeed, control theory offers several mature tools to deal with modelling, analysis, identification, control and observer design for general dynamical systems. However, several challenges arise when it comes to biological 1 A strain is a genetic variant or subtype of a microorganism. 1.1. Background and motivation processes:

1. Uncertainty: biological processes are inherently uncertain thanks to variations due to internal and external factors. Also, the models available for control and estimation often fail to provide an accurate description of the real process, since many complex functions are still not elucidated by biologists.

2. Monitoring: online measurements of some key variables are often not available, due to the lack of sensors. Depending on the type of measurement, this information might be subject to several sources of noise and also suffer from slow sampling (which might take several minutes, depending on the type of sensor).

3. Control: designing a controller for such systems is not a simple task thanks to its inherent complexity, uncertainty and non-linear nature.

4. Constraints: thanks to their physical meaning, these systems are highly constrained.

For instance, both state and input take only non-negative values.

5. Unusual supervision problems: frequently, the plant model in some operational mode becomes uncontrollable/unobservable (or close to it), which prevents a direct application of many existing approaches, since these properties often serve as obligatory prerequisites in theoretical studies.

There are two different ways of addressing control design for biological processes [START_REF] Hsiao | Control theory for synthetic biology: Recent advances in system characterization, control design, and controller implementation for synthetic biology[END_REF]: in vivo, which implements the control directly at a cell level (i.e., each cell has a copy of the controller), by using a set of chemical species and reactions, or in silico, which implements it at a population level by using a computer on the loop (see the illustration in Figure 1.1). While the first scenario highly constrains the control design, the second one flaws by offering control only up to a community-averaged sense.

Concerning the second approach, it has been shown that it is possible to control communities by using an external stimulus. For instance, the coexistence of co-cultures has been shown to be possible in an early work [START_REF] Davison | Effect of ph oscillations on a competing mixed culture[END_REF]] by oscillating pH levels and, similarly, in [START_REF] Krieger | Temperature regulation as a tool to program synthetic microbial community composition[END_REF] by cycling temperature levels. Microfluidics has also been used for control purposes [START_REF] Postiglione | Regulation of gene expression and signaling pathway activity in mammalian cells by automated microfluidics feedback control[END_REF]. Also, the use of optogenetics is popular, i.e., bio-engineering strains so it expresses certain genes when stimulated by light [START_REF] Toettcher | Lightbased feedback for controlling intracellular signaling dynamics[END_REF], [START_REF] Milias-Argeitis | Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth[END_REF], [START_REF] Scott | Biological signal generators: integrating synthetic biology tools and in silico control[END_REF].

However, despite the huge variety of existing control and estimation algorithms designed for biological processes, there is a persistent need for their development and adaptation to new tasks and tools in the domain.

The défi COSY

This thesis was hosted in the framework of the Inria project lab real-time COntrol of SYnthetic communities (COSY) 2 . The objective of this project was to promote the collaboration between different research teams and to foment interdisciplinary approaches to address synthetic microbial communities and their possible applications in biotechnology.

The central interest of this project was the development of a microbial consortium and its related mathematical model, as well as its subsequent analysis, and design of control and estimation algorithms. This model, presented in [START_REF] Mauri | Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs[END_REF], describes the interaction of two different strains a bacteria E. coli, where one of them is called the producer (since it produces some compound of interest, e.g., a heterologous protein), and a second one called the cleaner. This association is motivated by the fact that, while growing, the producer excretes a toxifying by-product (acetate) that inhibits its growth. The cleaner, having more affinity with the acetate, "cleans" the environment by consuming it and therefore alleviates inhibitory effects on the producer. The interest here is to observe if this association could lead to enhanced production of the compound of interest.

The design of control and estimation algorithms for this consortium based on that model is the primary objective of the present thesis.

Robust output feedback MPC

The control of constrained systems is known to be a hard (or even impossible) task to be tackled using classical feedback tools, such as the ones based on Lyapunov methods [START_REF] Mayne | Constrained model predictive control: stability and optimality[END_REF]. However, in practice, there exist several limitations on control (such as physical limitations on actuators) or states (such as safety margins or natural constraints, such as non-negativity), which can be usually met in the regulation of bioprocesses, as discussed in the previous subsection, for instance. In this light, Model Predictive Control (MPC) has been an active research topic over the last three decades since it offers a way to handle constrained, possibly non-linear and multivariate optimal control problems.

The MPC algorithm is intuitive: at each decision instant, the states are measured and the future behaviour of the system is predicted using a dynamical model. Using this prediction, an optimization problem is solved and the optimal input sequence that solves the control problem is obtained. Finally, the first control move of this sequence is applied to the system, and the procedure is iteratively repeated. However, although simple, two challenges arise with the implementation of this algorithm:

1. The model used for prediction is often uncertain due to unmodelled dynamics, parametric uncertainties, disturbances, and unaccounted delays. These sources of errors lead to discrepancies between the predicted and real behaviour of the system. These discrepancies might lead to the transgression of constraints and even to instability; 2. To achieve constraint satisfaction, MPC requires full state feedback -which is often unavailable. This fact leads to the need for estimation techniques which, consequently, increases the level of uncertainty due to inherent estimation errors;

An MPC scheme is said to be robust if it achieves the control task while robustly respecting all constraints for a given range of all these uncertainties. Furthermore, item (2) characterizes the robust output-feedback MPC (OF-MPC) problem and is the subject studied in this thesis.

Before formulating the thesis' problematics with more details, let us briefly survey the existing results and approaches for monitoring/control of bioprocesses and MPC design. Such a preliminary discussion will help us with the positioning of the thesis, indicating some of the remaining gaps that have to be filled.

State of the art

Control and observation of microbial communities

If only a single species is considered, the existing literature disposes of a myriad of results concerning control and observation. Indeed, the bioreactor model has been extensively used as a benchmark due to its high nonlinearities, yielding results on time-optimal [START_REF] Kokotovic | Time-optimal control for a model of bacterial growth[END_REF][D'ans et al., ] [Bayen et al., 2017]], sliding-mode [START_REF] Kravaris | Tracking the singular arc of a continuous bioreactor using sliding-mode control[END_REF], robust [Robledo, 2006], nonlinear feedback [START_REF] Karafyllis | A new small-gain theorem with an application to the stabilization of the chemostat[END_REF] and adaptive control [START_REF] Mailleret | Nonlinear adaptive control for bioreactors with unknown kinetics[END_REF]. In [START_REF] Mazenc | Stabilization in a chemostat with sampled and delayed measurements and uncertain growth functions[END_REF], the authors propose a nonlinear controller taking into account sampled and delayed measurements, as well as robustness against model uncertainty. The estimation problem (also called the monitoring problem, in which the objective is to design soft sensors) also disposes of a rich literature, such as presented in [START_REF] Bastin | On-line Estimation and Adaptive Control of Bioreactors[END_REF] and reviewed in [Dochain, 2003] [Mohd Ali et al., 2015].

On the other hand, if microbial consortia are considered, the results are more scarce. For instance, in vivo controllers have been proposed by [START_REF] Fiore | In-silico analysis and implementation of a multicellular feedback control strategy in a synthetic bacterial consortium[END_REF] and [START_REF] Kerner | A programmable escherichia coli consortium via tunable symbiosis[END_REF]. The use of in silico control (as well as experiments using such controllers) are yet rare [START_REF] Schlembach | Measurement techniques to resolve and control population dynamics of mixedculture processes[END_REF]. In [START_REF] Hoo | Global linearization and control of a mixed-culture bioreactor with competition and external inhibition[END_REF] the authors propose a control law of a community of yeasts (competing and externally inhibited) by means of global linearization.

In [START_REF] Treloar | Deep reinforcement learning for the control of microbial co-cultures in bioreactors[END_REF], a reinforcement learning technique has been used.

In this sense, the control and the estimation problems applied to the microbial consortia are less studied in the literature. Furthermore, the robust versions of these problems are still relatively open, especially in the context presented by the project COSY.

Robust output feedback MPC

The second part of this thesis is devoted to the robust output feedback (ROF) MPC problem, i.e., using a predictive controller to control a system that is subject to constraints (on control input and states) and does not dispose of full state measurement. The cases of linear time-invariant (LTI) and linear parameter-varying (LPV) systems are addressed.

For LTI systems, the ROF-MPC problem is very mature. Early results ranges from minmax optimization (i.e., the minimization of a cost function while considering worst-case of the disturbances) [START_REF] Bemporad | Min-max control of constrained uncertain discrete-time linear systems[END_REF], Raimondo et al., 2009], to the repeated solution of LMIs [Kothare et al., 1996, Rodrigues and[START_REF] Rodrigues | [END_REF]. Set-membership techniques, such as moving-horizon estimation (MHE) [START_REF] Bemporad | Output-feedback predictive control of constrained linear systems via set-membership state estimation[END_REF], [START_REF] Brunner | Enhancing outputfeedback MPC with set-valued moving horizon estimation[END_REF], [START_REF] Chisci | Feasibility in predictive control of constrained linear systems: the output feedback case[END_REF], have also been applied. These techniques use the available measurements to compute sets of admissible states (based on the available dynamics, bounds on the disturbances and initial uncertainty) and use them to guarantee constraint satisfaction.

Finally, Tube-based approaches [START_REF] Mayne | Robust output feedback model predictive control of constrained linear systems[END_REF]] [Mayne et al., 2009] [Langson et al., 2004] have also grown popular thanks to their reduced computational complexity: by predicting the trajectory of the nominal system, the controller is designed to enforce any deviation of the perturbed one to remain inside a tube.

For the LPV case, the problem becomes more complex since the scheduling parameter is unknown in future steps, making it harder to obtain reliable predictions. Nevertheless, several dynamic output feedback controllers have been proposed over the last years [Ding, 2010] [Ding et al., 2013]. In a different approach, [START_REF] Yang | An optimal approach to output-feedback robust model predictive control of LPV systems with disturbances[END_REF] proposes an observerbased technique relying on input-to-state stability and robust positively invariant sets of the estimation error. In [START_REF] Ping | An observer-based output feedback robust MPC approach for constrained LPV systems with bounded disturbance and noise[END_REF], the authors develop an approach that optimizes, simultaneously, both controller and observer. A tube-based method is presented by [START_REF] Yang | Tube-based output feedback model predictive control of polytopic uncertain system with bounded disturbances[END_REF]. Parameter-dependent Lyapunov functions have also been used [START_REF] Lee | Output feedback model predictive control for LPV systems using parameter-dependent lyapunov function[END_REF]. Min-max optimization has been utilized in [START_REF] Huang | Quasi-min-max dynamic output feedback MPC for LPV systems[END_REF] and [START_REF] Kim | Output-feedback model predictive control for LPV systems with input saturation based on quasi-min-max algorithm[END_REF], although no state constraints are imposed. Most of these works recursively update the estimation error sets and thus require the common assumption that the scheduling parameter is measured, (there are exceptions, e.g., [START_REF] Ding | Dynamic output feedback robust MPC with one free control move for LPV model with bounded disturbance[END_REF], [START_REF] Ding | Output feedback robust MPC for LPV system with polytopic model parametric uncertainty and bounded disturbance[END_REF], but with a constrained prediction horizon).

In the literature mentioned above, sets bounding the uncertain states or the estimation error are used. However, these techniques face some drawbacks such as numerical complexity and an intricate implementation, which can be relaxed by the use of interval tools [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF]. Especially concerning the LPV case, many of these techniques fix the prediction horizon to a one-step-ahead fashion (which can potentially shrink the feasible region of the optimization problem) and assume that the scheduling parameter is available, which is not always true in practice.

Thesis outline

In the light of the discussion above, the following problems were selected to be tackled in this thesis:

1. Estimate, in real-time, the concentration of the species in the consortium by using realistic measurements. These state estimators are then applied to the COSY problem; Control, pp. 2529-2534, 2020. 4 These works are not reported in this manuscript.

Part 2:

Control and Observation of Microbial Communities

Introduction to the overall framework

Introduction

In this chapter, the models used for the observer and control design, as well as their features, are introduced. Two models are considered: (i) the competition model, which describes the competition of two different species for a single limiting substrate, and (ii) the model developed by the COSY project, which describes the syntrophic relationship of two different strains. The main difference between these models is that, while in the former there is only one nutrient and there is competition, in the latter two nutrients are considered but the species do not compete directly for survival.

The competition model

The competition model is described by the following nonlinear differential equations [START_REF] Smith | The theory of the Chemostat[END_REF]:

dS(t) dt = D(t)(S in (t) -S(t)) -µ 1 (S(t))x 1 (t) -µ 2 (S(t))x 2 (t) dx 1 (t) dt = (k 1 µ 1 (S(t)) -D(t))x 1 (t), dx 2 (t) dt = (k 2 µ 2 (S(t)) -D(t))x 2 (t), (2.1)
where the states S and x i are, respectively, the concentrations of nutrient and the concentration of the i-th species. D is the dilution rate, which describes the rate of inflow and outflow of media of the bioreactor. S in is the concentration of nutrient diluted in the inflow of medium.

The functions µ i (S) are the reaction rates (also called specific growth rates) and describe the kinetics of nutrient uptake by each species. Finally, k i are yield coefficients, describing the portion of nutrient uptake that is used for growth by each strain.

This model (and also the one described in the next subsection) is a mass balance model.

From (2.1), one can see that the concentration of nutrient is increased by the inflow D(t)S in (t) and decreased by the consumption of the species (i.e., µ i (S)x i (t)). The concentrations of the species are positively influenced by this consumption. Finally, since a bioreactor usually operates with a constant volume, the outflow is equal to the inflow, meaning that all concentrations are negatively affected by the dilution rate D(t).

The main source of uncertainty in (2.1) comes from the description of nutrient uptake.

Indeed, the functions µ i might depend on a series of factors (such as substrate and population concentration, pH, temperature, inhibitors) [START_REF] Bastin | On-line estimation of microbial specific growth rates[END_REF], not always evident in experimentation. Furthermore, it is also hard to quantify how much of this uptake is indeed used for growth, leading the parameters k i to be also an approximation.

The functions µ i (S) may have many different forms. The most common one, known as Monod's law (or Michaelis-Menten kinetics), describes the growth of a species by consuming a single nutrient, being described by the following equation:

µ(S) = µ max S K m + S , ( 2.2)
where µ max > 0 and K m > 0 are the maximum reaction rate and the half-saturation constant.

Inhibitory effects (i.e., the decrease of the growth rate caused by high concentration of a nutrient) are usually modelled by the Haldane equation:

µ(S) = µ max S K m + S + S 2 K i . (2.3)
where K i > 0 is the inhibition constant. very important concept concerning competition models is the competitive exclusion principle:

the competition between two or more species for a single limiting substrate will result in the survival of, at most, one species.

This principle has been proven mathematically in [START_REF] Hsu | A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms[END_REF] for monotone response function, and in [START_REF] Wolkowicz | Global dynamics of a mathematical model of competition in the chemostat: General response functions and differential death rates[END_REF] for general functions.

The COSY model

The COSY model, presented by [START_REF] Mauri | Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs[END_REF], describes the dynamics of a producer strain that grows and produces a special compound by consuming glucose and, by doing so, excretes acetate (which is toxifying, in the sense that it inhibits its growth). The syntrophic relationship is established by adding a second cleaning strain, that has more affinity in consuming acetate and therefore "detoxifies" the environment and alleviates inhibitory effects on the producer growth (see the pathways illustrated in Figure 2.3 below). The producer and the cleaner strains are henceforth denoted B p and B c , respectively. The glucose and acetate are denoted by G and A, respectively, and the nutrient state vector is denoted vec(G, A). The equations of the model have the following form:

Ġ = D g G in -DG -r 1 (G, A)B p -r 2 (G, A)B c Ȧ = D a A in -DA -r 3 (G, A)B p -r 4 (G, A)B c Ḃp = (1 -k h )µ p (G, A)B p -(k + D)B p Ḃc = µ c (G, A)B c -(k + D)B c (2.4)
In this model, G in and A in are the concentration of each nutrient at the bioreactor inflow.

Also, D g and D a are the dilution rates of glucose and acetate, respectively, while k is a positive scalar representing degradation (or death rate).

In (2.4), the functions µ p and µ c are given by

µ p (G, A) = Y g r 1 (G, A) + Y a r 3 (G, A), µ c (G, A) = Y g r 2 (G, A) + Y a r 4 (G, A), (2.5)
where Y g and the Y a being the yield coefficients, denoting how much of each nutrient uptake is used for growth. Functions r i , i ∈ {1, . . . , 4}, similarly as for the competitive model, are the kinetic rates and denote the rates of exchange between the environment and the microbial cells and defined as

r 1 (G, A) = k p g G G + K g Θ a A + Θ a , r 2 (G, A) = k c g G G + K g Θ a A + Θ a , r 3 (G, A) = k a A A + K a Θ g r 1 (G, A) + Θ g -k over max(0, r 1 (G, A) -l), r 4 (G, A) = k a A A + K a Θ g r 2 (G, A) + Θ g + k ac A A + K ac -k over max(0, r 2 (G, A) -l)), (2.6)
and determine the nutrient uptake and acetate over-expression. The constant, positive, scalar parameters k * g , k a , k ac , k over are yield coefficients, while K a , K g , K ac , Θ g , Θ a are half saturation constants. The constant scalar l represents a threshold in which acetate starts to be secreted.

It is worth noticing that the term r 1 (G, A) indicates that the producer grows by consuming glucose (the first term), while it is inhibited by the concentration of acetate (the second term).

An illustration of this effect is depicted in Figure 2.2. Also, from r 3 (G, A) and r 4 (G, A), it is readily seen that some terms enter positively in the equation for A in (2.4), showing that acetate might be increased by excretion of the microbes. Along with (2.4), additional dynamics can be considered to express production by the producer:

Ḣ(t) = k h µ p (G(t), A(t))B p (t) -DH(t), (2.7) Ḟ (t) = mH(t) -D(t)F (t). (2.8)
where H is the product (e.g., a heterologous protein) and F is a by-product (for instance, a fluorescent reporter resulting of the maturation of H) and k h , m are proportionality constants. Remark 2.1. The dynamics (2.7) basically means that the uptake of substrate G is used for another purposes other than growth. Therefore, similar additional dynamics can be also considered for the competition model (2.1).

Interesting features of the model (2.4) have been pointed out in [START_REF] Mauri | Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs[END_REF]. The coexistence of producer and cleaner strains depends on the growth parameters (such as the product yield k h ), as well as on the operating conditions (i.e., the dilution rate D and nutrient inflow concentration G in ). The cleaner strain grows, primarily, by consuming acetate. Although model (2.4) foresees an input of this nutrient (through the term D a A in ), its concentration in the bioreactor is mainly due to secretion by the producer strain, which might only occur if r 1 (G, A) > l or r 2 (G, A) > l. Therefore, the persistence of the cleaner strain in the bioreactor is guaranteed only if the producer grows sufficiently fast. Furthermore, since this growth is affected by k h , this coefficient also needs to be taken into account when coexistence is analyzed.

Figure 2.4 illustrates the simulation of (2.4) in two scenarios: one with

D = 0.2 [h -1 ] and with D = 0.3 [h -1 ], both considering G in = 20 [gL -1
] and A in = 0. In the left plot, the trajectories of B p and B c are shown and it is clear that, with a lower dilution rate, the cleaner is washed out (i.e., B c → 0). This is due to the insufficient acetate excretion by the producer species: the glucose inflow is not enough to ensure acetate overflow (i.e., if r 1 (G, A) > l).

This acetate production (as illustrated by the solid line in the right plot) allows the dynamics of B c to converge to a positive steady-state.

In the same scenarios, Figure 2.5 illustrates the trajectories of the product H. As it can be seen, although B p reached a higher concentration with D = 0.2 [h -1 ], the scenario in which coexistence happens has yielded a higher production of H.

Following [START_REF] Mauri | Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs[END_REF], Figure 2.6 illustrates the relation between the inputs G in and D in the steady-state values of each species. As it can be seen, there is a narrow region of D that permits coexistence. Obviously, for too high values of D, no species can survive and the biomass is washed-out of the bioreactor.

Clearly, the investigation above is carried out in an open-loop setup. Therefore, inspired by these facts, it is of interest to study the use of closed-loop control algorithms to control this consortium. In addition, the objective is to evaluate the conditions in which the consortium will outperform the producer-only scenario. The control algorithm development also aims to enhance the operation of the consortium in terms of performance and robustness.

In the two next chapters, the estimation and the control algorithms are discussed. First, the conditions of observability and the subsequent state observers are proposed (these results can be applied to both systems). Then, the control algorithms are discussed: first a robust approach for the competition model, then the (nominal) control of the COSY model. For the latter, a nominal scenario (i.e., all parameters and variables are perfectly known) is considered. In this chapter, the design of estimation algorithms is presented for the systems described in Chapter 2. First, the available measurements are introduced, then the questions addressed in this chapter are posed.

Available measurements

In this chapter, the following measurements will be considered:

1. The total biomass y 1 (t) = n i=1 x i (t): obtained by methods such as flow-cytometry or mass-spectrometry. These kinds of measurement are usual in practice, although demanding a higher sampling time;

2. The fluorescent reporter y 2 (t) = F (t): obtained by optical density and obeying dynamics (2.7). Since it uses light to promote an indirect measurement, it can be sampled faster.

However, the strain that produces the protein and expresses this reporter needs to be bio-engineered.

In this light, this chapter aims to address the following questions: Q1-Are systems (2.1)- (2.4) completely or partly observable? Q2-What are the necessary conditions for observability? Q3-How do the available measurements y 1 (t) or y 2 (t) influence observability and observer design? Q4-How to design the appropriate observers?

Preliminaries

Before discussing the questions posed above, some preliminary tools are introduced.

Observability of nonlinear systems

Observability is a property of a dynamical system indicating how well its internal states can be inferred by means of the available measurements. For nonlinear systems, the problem of observability has extensively discussed in the seminal paper by [START_REF] Hermann | Nonlinear controllability and observability[END_REF].

Following that work, let us consider a general nonlinear system endowed with a measurement function as

ẋ = f (x, u) y = h(x) (3.1)
where x ∈ M , a smooth connected manifold of dimension n, y ∈ R p , u ∈ Ω ⊆ R m and f, g being smooth functions. Let Σ(x 0 , u, [t 0 , t 1 ]) denote the solution of (3.1), in the time span [t 0 , t 1 ], with initial condition x 0 and input u. Then, a pair of points x 0 , x 1 are said indistinguishable if they realize the same input-output map, i.e., for every admissible input

u(t) in [t 0 , t 1 ], h(Σ(x 0 , u, [t 0 , t 1 ])) = h(Σ(x 1 , u, [t 0 , t 1 ])).
Let I(x 0 ) define the set of all indistinguishable points at x 0 . Then, observability is defined as follows Definition 3.1. [START_REF] Hermann | Nonlinear controllability and observability[END_REF]

System (3.1) is said observable at x 0 ∈ M if I(x 0 ) = {x 0 }, and observable if I(x) = {x} for every x ∈ M .
The definition above is global, however, it can be also stated locally (i.e., in a open neighborhood U of x 0 ). Furthermore, it might be only necessary to distinguish x 0 from its closest neighbors, leading to the weaker concept of observability called weak observability. Definition 3.2. [START_REF] Hermann | Nonlinear controllability and observability[END_REF] System (3.1) is said to be (i) locally weakly observable at x 0 ∈ M if there exists an open neighborhood U ⊂ M of x 0 such that I(x 0 ) = {x 0 } for every open neighborhood V ⊂ M of x 0 contained in U , and (ii) locally weakly observable if it is so for all x ∈ M .

To introduce a criteria for verification of observability for (3.1), consider the following definition of the Lie derivatives: Definition 3.3. (Lie Derivative): Let h(x) be a continuously differentiable function and f (x, u) be a vector field. The Lie derivative of h(x) along f (x, u) is defined as

L f h(x) = L 1 f h(x) = ∂h(x) ∂x f (x, u), L k f h(x) = L f L k-1 f h(x), k ≥ 2.
Then, the observation space O be the set formed by all linear combinations induced by the Lie derivatives of h i (the components of h) with respect to f , i.e.,

O = [h 1 , . . . , h p , L f h 1 , . . . , L f h p , L 2 f h 1 , . . . , L 2 f h p , . . . ].
Let dO be the gradient of the elements of O. The well-known observability rank condition is stated as follows:

Theorem 3.1. [START_REF] Hermann | Nonlinear controllability and observability[END_REF] The system (3.1) is (locally) weakly observable if dO contains n linearly independent vectors, i.e., rank(dO) = n.

If the system under consideration is time-varying and autonomous (as it will be the case in the following), i.e.,

ẋ = A(t)x(t) y = C(t)x(t)
and considering that A(t), C(t) are analytic functions, then the following algorithm [Sontag, 1998, Ch. 6, p. 279] can be used to verify the rank condition above:

rank                   C 0 (t) C 1 (t) . . . C n-1 (t)                   = n,
where

C 0 = C(t) and C i+1 (t) = d dt C i (t) + C i (t)A(t) for i = 1, . . . , n -1.

Estimating microbial co-cultures

Observability analysis

Since the structure of systems (2.1) and (2.4) are similar, the observability analysis presented in this section will be derived only for the latter [Souza et al., 2021a]. Indeed, the observability conditions for the species concentration are identical (by substituting 1 -k h by k 1 , and µ p , µ c by µ 1 , µ 2 , respectively). The observability condition for the nutrients are different but, nonetheless, follow the same rationale in both cases.

Consider (2.4) decomposed in two sub-systems: the first composed of the strains and products and the second composed of the nutrients. Denoting x 1 = vec(B p , B c , H, F ) and

x 2 = vec(G, A), the first sub-system can be obtained by rewriting (2.4)

as ẋ1 =         (1 -k h )µ p (x 2 ) -(k + D) 0 0 0 0 µ c (x 2 ) -(k + D) 0 0 k h µ p (x 2 ) 0 -(k + D + m) 0 0 0 m -(k + D)         x 1 y =   0 0 0 1 1 1 0 0   x 1 (3.2)
where µ p (x 2 ), µ c (x 2 ) are given by (2.5). In (3.2), the variables G, A (i.e., the elements of x 2 ) are interpreted as external inputs.

Proposition 3.1. Assume D > 0 and that k, m, k h and µ p (x 2 ), µ c (x 2 ) are known. Then,

system (3.2) is observable in an interval [t 0 , t 1 ] if there exists a t ∈ [t 0 , t 1 ] such that (1 - k h )µ p (x 2 (t)) ̸ = µ c (x 2 (t)).
Proof. Since system (3.2) can be interpreted as an autonomous (without inputs) time-varying dynamics, and endowed with a constant measurement matrix, the observability condition to be satisfied in the interval [t

0 , t 1 ] is ∃t ∈ [t 0 , t 1 ] such that rank                 0 0 0 1 1 1 0 0 0 0 m -(k + D) (1 -k h )µ p (x 2 ) -(k + D) µ c (x 2 ) -(k + D) 0 0                 = 4.
Since m, k, k h are positive constants and D > 0 by assumption, this condition is verified

if (1 -k h )µ p (x 2 ) ̸ = µ c (x 2 ).
If the conditions given in Proposition 3.1 are verified, then B p and B c are observable and, therefore, can be considered as available signals. As readily seen in ( 2.4), the state vector x 2 and the measurement y 2 (t) are coupled by the functions µ p (x 2 ) and µ c (x 2 ) through the first time derivative of x 1 .

The fact above simplifies the observability analysis of the pair x 2 since all terms in (3.2), except the arguments of functions µ p (x 2 ), µ c (x 2 ), are known. Denoting the vector-valued

function µ(x 2 ) = vec(µ p (x 2 ), µ c (x 2 )
), this brings us to the following proposition:

Proposition 3.2. Let the conditions given in Proposition 3.1 be satisfied. Then, the state

vector x 2 is observable in a neighborhood of a point x 0 2 if, for x 1 ̸ = 0, the following condition holds det J(µ(x 0 2 )) ̸ = 0
where J(µ(x 0 2 )) is the Jacobian of µ(x 2 ) evaluated at x 0 2 . Proof. As discussed above, if x 1 and its time derivative are available, the observability problem for x 2 reduces to the invertibility of the vector-valued function µ(x 2 ). Evoking the inverse function theorem [Nijenhuis, 1974], the condition given in this proposition assures that µ(x 2 ) is invertible in, at least, a vicinity of the point x 0 2 . Remark 3.1. If only one substrate is considered, i.e., x 2 = S (such as in the case of model (2.1)), then the condition given in Proposition 3.2 reduces to the invertibility of function µ i (S).

A question that might arise concerns the observability conditions if solely y 2 is available.

To tackle this question, a weaker concept of observability, called detectability (or asymptotic observability) can be evoked. Recalling [Sontag, 1998], a system is called detactable if all of its unobservable states are asymptotically stable. This brings us to the following proposition: Proposition 3.3. If only the measurement y 2 is available, then system (2.4) is detectable.

Proof. As it can be directly seen from (2.4), the states H and F are not observable through y 2 , since no relation with this measurement or its derivatives can be directly constructed.

For the remaining states, however, the observability conditions given in propositions 3.1-3.2 remain valid. Indeed, if a subsystem composed of B p and B c is considered, the rank condition simply reduces to rank

    1 1 (1 -k h )µ p (x 2 ) -(k + D) µ c (x 2 ) -(k + D)     = 2
which does not alter the requirement stated in Proposition 3.1.

Noticing that the dynamics of the unobservable states F and H are basically first-order stable filters having B p and x 2 as inputs, this fact implies that the unobservable states have asymptotically stable dynamics. Therefore, system (2.4) is detectable with respect to y 2 .

Using the conditions given in propositions 3.1-3.3, the observable (or detectable) states of (2.4) might be numerically estimated.

Estimation of H, B p and B c

In this subsection, the objective is to develop state observers for system (3.2). The estimation of the nutrients G, A will be addressed in the subsequent section. To construct these observers, the differentiation of the available measurements will be needed, hence it is assumed that the input D is continuously differentiable with a known and bounded derivative.

As it will be shown in the following, the first and second time derivatives of y 1 , as well as the first time derivative of y 2 , will be needed. Therefore, the following differentiators (see Appendix B) will be implemented:

             ż0 = λ 12 L 1 3 1 |z 0 -y 1 |sign(z 0 -y 1 ) + z 1 ż1 = λ 11 L 1 2 1 |z 0 -y 1 |sign(z 0 -y 1 ) + z 2 ż2 = λ 10 L 1 sign(z 0 -y 1 ) (3.3)      ξ0 = λ 21 L 1 2 2 |ξ 0 -y 2 |sign(z 0 -y 2 ) + ξ 1 ξ1 = λ 20 L 2 sign(ξ 0 -y 2 ) (3.4)
and thus, according to [Levant, 2003], if the gains L 1 and L 2 are selected such as

L 1 ≥ m k h µ 1 B p -mH -2(k + D)H -(F + w) Ḋ + (k + D) 2 + (k + D)m , L 2 ≥ (1 -k h )µ 1 B p + µ 2 B c -(k + D)(B p + B c + w) , (3.5)
then, it is guaranteed that z 0 → y 1 , z 1 → ẏ1 , z 2 → ÿ1 , and ξ 0 → y 2 , ξ 1 → ẏ2 in a finite-time in the noise-free case, and the related discrepancies stay bounded for bounded noises. Having these estimates, an observer for H is readily obtained from the dynamics of F in (2.4) as

follows: Ĥ = z 0 (k + D + m) + z 1 m . (3.6) Remark 3.2.
The quantities in (3.5) can be estimated by determining the range of operation of the bioreactor.

Since the dynamics of both B p and B c depends on G and A (which are, a priori, unknown),

denote µ p = µ 1 (G, A)B p and µ c = µ 2 (G, A)B c . Computing the second derivative of F , one gets F = m Ḣ -(k + D) Ḟ = m [k h µ p -(k + D + m)H] -(k + D) Ḟ -ḊF,
and thus, the measurement y 1 (t) can be used to estimate µ p (G, A) by means of the output of the differentiator (3.3) and the estimate (3.6):

μp = z 2 + (k + D)z 1 + Ḋz 0 mk h + (k + D + m) Ĥ k h = z 2 + [2(k + D) + m]z 1 + Ḋ + m + (k + D) 2 z 0 mk h . (3.7)
Now, let us investigate the dynamics of the total biomass. By summing the third and forth equations in (2.4), one gets

Ḃp + Ḃc = (1 -k h )µ p (G, A)B p + µ c (G, A)B c -(k + D)(B p + B c ),
which allows us to use the measurement y 2 , along with the output of the differentiator (3.4)

and the estimate (3.7), to compute an estimate of µ c by

μc = ξ 1 + (k + D)ξ 0 -(1 -k h )μ p .
Therefore, profiting on the structure of Ḃp and Ḃc , the following observers can be proposed: Using the Young's inequality [START_REF] Hardy | Inequalities. Cambridge Mathematical Library[END_REF], the term e p (µ p -μp ) can be upper bounded as follows:

Bp = (1 -k h )μ p -(k + D) Bp , ( 3
e p (µ p -μp ) ≤ ϵ 2 2 e 2 p + (µ p -μp ) 2 2ϵ 2
for any constant ϵ > 0. Applying the same idea to the term e c (µ c -μc ), the derivative V can be upper bounded by

V ≤ -2(k + D -ϵ 2 )V + (µ p -μp ) 2 2ϵ 2 + (µ c -μc ) 2 2ϵ 2 . ( 3

Estimation of G and A

An interesting feature of the observers proposed in the previous section is that they do not depend on any of the reaction rates r i and, therefore, do not require the joint estimation of the nutrients G and A.

Obviously, this is not the case when G and A need to be estimated -what would be of interest when designing control algorithms for system (2.4) (or (2.1)). Indeed, real-time information of these quantities is hard to obtain in practice, making this estimation an interesting option. In this light, following hypothesis is imposed:

Assumption 3.1. The estimates Bp and Bc have converged to the true values of B p and B c .

Assumption 3.1 is technical and is needed only to simplify the analysis of the observers proposed in the following. In practice, all observers will be launched at once. Under assumption 3.1, the following observer can be proposed

Ĝ = -r 1 ( Ĝ, Â)B p -r 2 ( Ĝ, Â)B c -D Ĝ + D g G in  = -r 3 ( Ĝ, Â)B p -r 4 ( Ĝ, Â)B c -D  + D a A in . (3.11)
For the next theorem, let the estimation errors be given by e g = G -Ĝ and e a = A -Â.

Theorem 3.3. Let Assumption 3.1 be satisfied and consider system (2.4) and observer (3.11).

If D > (B c + B p ) max k p g + k c g Θ a , k a k g Θ g K 2 g
, then the dynamics of the estimation errors e g , e a are globally asymptotically stable.

Proof. The dynamics of the estimation errors e g , e a are computed as

ėg = (-r 1 (G, A) + r 1 ( Ĝ, Â))B p + (-r 2 (G, A) + r 2 ( Ĝ, Â))B c -De g ėa = (-r 3 (G, A) + r 3 ( Ĝ, Â))B p + (-r 4 (G, A) + r 4 ( Ĝ, Â))B c -De a .
(3.12)

The objective here is to show that e g is a perturbation on the dynamics of e a (and viceversa), and the overall interconnection is ultimately bounded. Considering the term proportional to B p , summing and subtracting r 1 ( Ĝ, A) leads to

-r 1 (G, A) + r 1 ( Ĝ, Â) = -r 1 (G, A) + r 1 ( Ĝ, Â) + r 1 ( Ĝ, A) -r 1 ( Ĝ, A) = k p g Θ a A + Θ a Ĝ K g + Ĝ - G K g + G + k p g Ĝ K g + ĜΘ a 1 Θ a + Â - 1 Θ a + A = k p g Θ a A + Θ a K g -e g (K g + Ĝ)(K g + G) + k p g Ĝ K g + ĜΘ a e a (A + Θ a )( Â + Θ a )
.

Since the function r 2 is similar to r 1 , an analogous result can be derived for the term proportional to B c . Considering a Lyapunov function candidate given by V g = |e g |, a bound on its time derivative can be computed as

Vg ≤ -D|e g | + B p k p g + B c k c g Θ a |e a |. (3.13)
Now, let us focus on the dynamics of e a in (3.12). Note that, due to the over-expression terms (i.e., r over, * a (G, A)), r 3 and r 4 are composed of two functions each:

-r 3 (G, A) + r 3 ( Ĝ, Â) = -r 3 (G, A) + r 3 ( Ĝ, Â) + r 3 ( Ĝ, A) -r 3 ( Ĝ, A) = r up,p a ( Ĝ, Â) -r up,p a (G, A) + r over,p a (G, A) -r over,p a ( Ĝ, Â).
Consider the difference r up,p a ( Ĝ, Â) -r up,p a (G, A) as above. Summing and subtracting

r up,p a ( Ĝ, A), one can write r up,p a ( Ĝ, Â) -r up,p a (G, A) = k a r up,p g (G, A) + Θ g -e a (A + K a )(  + K a ) + k a   + K a Θ g   r up,p g (G, Â) -r up,p g ( Ĝ, Â) (r up,p g ( Ĝ, Â) + Θ g )((r up,p g (G, Â) + Θ g )   = k a r up,p g (G, A) + Θ g -e a (A + K a )(  + K a ) + k a   + K a Θ g Θ a  + Θ a   k p g eg ( Ĝ+Kg)(G+Kg) (r up,p g ( Ĝ, Â) + Θ g )((r up,p g (G, Â) + Θ g )   .
Proceeding analogously for the remaining difference r over,p a (G, A) -r over,p a ( Ĝ, Â) one has

r over,p a (G, A) -r over,p a ( Ĝ, Â) = k over r up,p g (G, A) -r up,p g ( Ĝ, Â) = k over k p g G G + K g Θ a -e a (A + Θ a )(  + Θ a ) + k p g Θ a  + Θ a K g e g (K g + G)(K g + Ĝ)
.

The last term to be analyzed is the one proportional to B c , i.e., r 4 (G, A) -r 4 ( Ĝ, Â).

However, due to the similarity of the concerned functions, one can profit of the developments made above to write

-r up,c a (G, A) + r up,c a ( Ĝ, Â) = k a r up,p g (G, A) + Θ g -e a (A + K a )(  + K a ) + k a   + K a Θ g Θ a  + Θ a   k c g eg ( Ĝ+Kg)(G+Kg) (r up,c g ( Ĝ, Â) + Θ g )((r up,c g (G, Â) + Θ g )   , and r over,c a (G, A) -r over,c a ( Ĝ, Â) = k over k c g G G + K g Θ a -e a (A + Θ a )(  + Θ a ) + k c g Θ a  + Θ a K g e g (K g + G)(K g + Ĝ)
.

Finally, considering a Lyapunov function candidate V a = |e a |, whose time derivative can be upper-bounded using the developments made above as

Va ≤ -D|e a | + B p k a k p g Θ g K 2 g + k p g K g + B c k a k c g Θ g K 2 g + k c g K g |e g | ≤ -D|e a | + k a K g 1 + Θ g K g Θ g K g (k p g B p + k c g B c g ) |e g |.
(3.14)

In order to prove boundedness of the overall estimation error, let V ag = V g + V a . The time derivative of this function, according to (3.13) and (3.14), is upper bounded by

Vag = Vg + Va ≤ -D(|e g | + |e a |) + B p k p g + B c k c g Θ a |e a | + k a K g 1 + Θ g K g Θ g K g (k p g B p + k c g B c g ) |e g | ≤ -DV ag + B p k p g + B c k c g Θ a |e a | + k a K g 1 + Θ g K g Θ g K g (k p g B p + k c g B c g ) |e g |, leading us to conclude that, if D > (k p g B c + k c g B p ) max 1 Θa , ka(1+ΘgKg) ΘgK 2 g
, then V a,g is a Lyapunov function, proving the claim.

The advantage of measuring y 1

As discussed in Proposition 3.3, the system (2.4) is detectable through y 2 . Then, a question is raised: why not design an observer for the observable states (in this case, the concentration of substrates and microbes) using solely this measurement? By substituting B p and B c by Bp and Bc in (3.11), additional observer equations can be proposed:

d dt B =   (1 -k h )µ p (x 2 ) -(k + D) 0 0 µ c (x 2 ) -(k + D)   B + L(y 2 -C 2 B), d dt Ĥ = k h µ p (x 2 ) -(k + D + m) Ĥ, (3.15) for B = vec( Bp , Bc ), C 2 = 1 1 , x2 = vec( Ĝ, Â)
, and some gain L. In this case, the estimates Bp and Bc would be computed jointly with the estimates Ĝ and Â. Obviously, a more intricate analysis would have to be carried out to identify the stability conditions.

It is interesting to notice that, in addition to allowing the observation of H, the measurement y 1 is helpful when observing B p and B c in the sense that it eliminates the need of knowing the rates r i perfectly.

Therefore, one can conclude that estimation of H, B p and B c using both y 1 and y 2 requires only knowledge on parameters m, k, and k h , while the estimation of G and A requires knowledge on all parameters except Y g and Y a . In counterpart, if solely y 2 is used, then the knowledge of all parameters in (2.4) are required for the estimation of the (observable) state variables. Clearly, in practice, these parameters (especially those related to the functions r i , i ∈ {1, . . . , 4}) may be uncertain and, thus, the first scenario would introduce less error in the resulting estimation.

Numerical example

In this section, a numerical experiment is presented to illustrate the proposed methodology.

For the simulation of system (2.4), the parameters are such as given given in [START_REF] Mauri | Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs[END_REF] and, arbitrarily (but realistically) k h = 0.2 and m = 1. The initial conditions are selected as

G(0) = 5, A(0) = 0, B p (0) = 3, B c (0) = 2, H(0) = 0 and F (0) = 0 [gL -1 ].
Furthermore, the inlet glucose concentration is taken as G in = 20 [gL -1 ] and the dilution rate as D = 0.3 + 0.05 sin(0.1t) [h -1 ]. Aiming at a simulation that is closer to a real-world scenario, it is considered that all measurements are sampled at every 10 minutes by a zeroorder holder. Figure 3.1 shows the sampled measurements y 1 and y 2 , both corrupted by measurement noise (in contrast with the unperturbed signals, depicted by the black curves).

In the following, the simulation results will be presented for two scenarios: using both measurements and applying observers (3.6), (3.8) and (3.9), and using only the measurement of the total biomass (i.e., y 2 ), applying observers (3.15). A discussion on the estimation error introduced by parametric uncertainty is then provided. Finally, the initial conditions for the observers are taken as Ĝ(0) = 0, Â(0) = 0, Bp (0) = 0 and Bc (0) = 0. which are obtained through integration, it does not require any knowledge of the rates r i .

Using only the total biomass

Now, the same scenario described previously will be considered, but using observers given by (3.15). These observers will be initialized with Ĝ(0) = 0, Â(0) = 0, Bp (0) = 0, Bc (0) = 0 and Ĥ(0) = 0. The observer gain will be selected as L = 10. measurement noise is attenuated. The same conclusion is drawn for the estimate of H, as shown in Figure 3.7.

Part conclusion

In this chapter, the observation problem for systems (2.4) (and also for (2.1), thanks to their similarities) has been discussed. First, an analysis has been carried out to identify their conditions of observability considering the total biomass and the fluorescent reporter as available measurements. It has been shown that the system is observable with both measurements, but merely detectable if only the first is available.

Second, considering the two available measurements, state observers were proposed for each variable, and the stability properties are discussed for the case in which the two measurements are available. Finally, the advantage of measuring the fluorescent reporter is highlighted, in the sense that it not only enables observability of the system, but also alleviates the dependence on the parameters of the system. In this chapter, the generic bioreactor model (2.1) and the COSY model (2.4) are reconsidered to address the control design for a co-culture of microbes. Mainly, the interest of designing a controller for a heterogeneous community is to ensure their co-existence, i.e., avoiding that one goes extinct. It has also been shown that an oscillating chemostat [Smith, 1981], or the use of feedback control [START_REF] Leenheer | Feedback control for chemostat models[END_REF], might induce the coexistence of two different species.

Also, another control objective relates to the stabilization (or the optimization) of the production on the bioreactor. Indeed, controlling production requires, in some sense, avoiding that the producing species are outcompeted and that the required nutrient is depleted from the bioreactor (if there is no nutrient, there is no metabolism).

Therefore, the main questions to be tackled in this chapter are:

Q1-Are these systems controllable (or stabilizable) taking into account the control constraints?

Q2-How can the co-population be (robustly) stabilized at desired concentration levels?

Q3-How can the production of the bioreactor be (robustly) stabilized? Is it possible to cast an optimization problem?

These questions are going to be discussed in the sequel, by considering both the generic model (2.1) and also the COSY model (2.4). Before starting to tackle these objectives, a discussion on the available control inputs are presented in the following.

The available control inputs

In both models (2.1) and (2.4), there are two controlled variables, both physically altering the composition of the environment in the bioreactor:

1. D(t): describing the dilution rates, i.e., the rates of transfer (inflow or outflow) of medium in the bioreactor (in unit h -1 ). This variable is often actuated by an electromechanical pump. This is interesting since it allows the design of a control considering variables of similar nature. Furthermore, the time constant of a electromechanical device is much smaller than the one of a biological process.

S in (t) = S max in D s (t) D s (t) + D 0 (t)

A discussion on controllability

Thanks to its physical meaning, systems (2.1) and (2.4) are time-varying, positive systems. Furthermore, they have constrained controls, since D(t) and S in (t) (or G in (t) and A in (t) for the COSY model (2.4)) are also non-negative and possibly upper bounded. In such a case, deriving a rank condition for controllability is not applicable since it cannot be investigated by utilizing Lie algebra, as discussed for general nonlinear systems in [START_REF] Hermann | Nonlinear controllability and observability[END_REF].

Nevertheless, an analysis can still be performed by profiting the structure of (2.1) and (2.4) and it is possible to evaluate the reachability of such a system, according to the following definition: Definition 4.1. A state x 1 is called reachable if there exist a finite time t 1 and an admissible control u(t) on [t 0 , t 1 ] such that x(t 0 ) can be steered to x 1 .

Consider the competition model (2.1). If the nutrient inflow concentration S in can be made sufficiently high, it is possible to prevent the nutrients of being depleted from the bioreactor (i.e., S → 0). Then, the sign of the derivative of each x i can be assigned by a proper selection of D.

Consider the sub-system for x 1 and x 2 in (2.1). Clearly, it can be rewritten as autonomous time-varying linear system, whose solution is

x 1 (t) = e f 1 (t-t 0 ) B p (t 0 ), x 2 (t) = e f 2 (t-t 0 ) B c (t 0 ), (4.1)
where

f i = t t 0 k h µ i (S(s)) -(k + D(s))ds, for i = 1, 2.
From these expressions and information on the functions µ i , it is possible to estimate the reachable set for the sub-system above.

Consider now the COSY model (2.4). Although three inputs are considered available (i.e., the dilution rate D and the inflow substrate concentrations G in and A in ), a similar analysis is more intricate due to the presence of two substrates. Also, considering D = D a = D g and recalling the dynamics for the substrates:

Ġ = DG in -DG -r 1 (G, A)B p -r 2 (G, A)B c Ȧ = DA in -DA -r 3 (G, A)B p -r 4 (G, A)B c
the two right-most terms on the dynamics of G are always negative. This is important since the input DG in , being always non-negative, can be used to compensate these terms and assign the sign of Ġ.

However, the same does not apply to the two right-most terms on the dynamics of A: thanks to the metabolic overflow, acetate is excreted by the species (mostly by the producer strain) after a certain threshold depending on G, A (see (2.6)). Therefore, it can be seen as an extra input of acetate, what might eventually change the sign of these terms.

Clearly, it would be possible to compensate this by a proper selection of D. Nevertheless, if D is fixed (or used for another purposes, as it will be the case in this chapter), it would be impossible to freely assign the sign of Ȧ and drive the concentration of A by using solely the input DA in , since it is always non-negative.

Controlling competing species 4.3.1 A gentle motivation for a robust approach

Before presenting the proposed control architecture, the nuances of the problems are stressed. Consider (2.1) with a single limiting substrate S. As early pointed out by [START_REF] Leenheer | Feedback control for chemostat models[END_REF], coexistence is possible if µ 1 (S) and µ 2 (S) intersect at a point S i . Selecting

D(t) = µ(S i ) implies that ẋi = (µ i (S(t)) -D(t))x i = 0.
However, if each of the competing species are to be stabilized at given concentration levels

x d i , i = {1, 2}
, the control problem becomes more complex. Indeed, this problem can be seen as the stabilization of the origin of the auxiliary dynamics σ(

x i ) = x i (t) -x d i . Computing the time derivative of σ(x i ), one has σ(x i ) = (µ i (S(t)) -D(t))x i .
Evidently, it is possible to steer only one of the species to a certain point in the statespace by setting D(t) = µ i (S). The other species, depending on the functions µ i (S(t)), will be driven to another point and might even be extinct. However, consider that D(t) is set such that one of the species (henceforth, arbitrarily, x 1 ) is stabilized at a given point. Denoting ξ(S) = µ 1 (S) -µ 2 (S) and computing the first and the second time derivative of σ(x 2 ), one

has that σ(x 2 ) = (µ 2 (S) -µ 1 (S))x 2 = -ξ(S)x 2 (4.2) σ(x 2 ) = - ∂ξ ∂S Ṡ + ξ 2 x 2 = x 2 - ∂ξ ∂S µ 1 (S)(S in -S) -k 1 µ 1 (S)x 1 -k 2 µ 2 (S)x 2 + ξ 2 (S) , (4.3)
which indicates that S in can be utilized to stabilize the surface σ(x 2 ) = 0. In this light, the following challenges are met: 3. From (4.3) one observes a heavy dependence on the parameters of µ i (S), which, as discussed previously, are uncertain. This fact motivates the need for a robust approach.

In the remainder of this section, the design of a control architecture aiming to solve this problem is proposed. This control architecture is composed of two phases: in the first moment, the control D(t) is used to stabilize x 1 at the desired position x d 1 > 0 and, simultaneously, the substrate concentration will be shifted to a positive value ensuring the system controllability using the control S in (t). Once these goals are reached, the control goal for S in (t) will commute to the stabilization of x 2 at a constant level.

In this section we will assume that

µ i (S) = a i b i + S , i = 1, 2,
and the sets of admissible values A i and B i are given for a i and b i , respectively, then the following hypotheses are imposed: Hypothesis 4.1 states that the kinetic rates intersect, which is an obligatory requirement to have co-existence of the species (see [START_REF] Leenheer | Feedback control for chemostat models[END_REF] for an extensive discussion), while the latter property can be guaranteed by a proper numbering of each species. Hypothesis 4.2 implies that the function ξ(S) and its derivative with respect to S are not zero at the same point. As it will be shown in the design of the control architecture, the distinction of these points is needed to achieve the stabilization of x 2 (t). For further reference, define

µ i (S) = min a i ∈A i , b i ∈B i µ i (S)
, and µ i (S) = max

a i ∈A i , b i ∈B i µ i (S).

Design of the control architecture

Stabilization of x 1

This subsection will report the results of [START_REF] Souza | Robust stabilization of competing species in the chemostat[END_REF]. Consider the following decision variable:

σ 1 (x 1 ) = x 1 -x d 1 (4.4)
where x d 1 > 0 is an arbitrary point for x 1 to be stabilized. As it can be seen from equation (2.1), if the concentration x i is in a steady-state value, i.e., ẋi (t) = 0 for t ≥ 0, then there are two possibilities: either x i = 0 or D(t) = µ i (S(t)). With the latter case in mind, the following proposition is stated: Proposition 4.1. Consider model (2.1) and the surface (4.4). If the dilution rate is selected such as

D(t) =      µ 1 (S) + χ, if σ 1 ≥ 0 µ 1 (S) -χ, if σ 1 < 0 (4.5)
where χ ∈ (0, µ 1 (S(0))) is a tuning parameter and, then the closed-loop system is globally finite-time stable with respect to the output σ 1 , provided that S(t) ≥ S(0) ∀t ≥ 0.

Proof. Consider the following Lyapunov function candidate for (4.4), which is proper and positive definite with respect to σ 1 :

V 1 = 1 2 σ 2 1 (4.6)
whose time derivative is given by V1 = σ 1 σ1 . By noticing that σ 1 = √ 2V 1 , this relation can rewritten such that D appears explicitly:

V1 = σ 1 (µ 1 (S) -D(t))x 1 (t)
and hence, selecting D(t) as given by (4.5) and stressing that x 1 (0) can be either lower or

greater than x d 1 , then V1 ≤ -χ √ 2V 1 min{x 1 (0), x d 1 }
holds, proving the claim. An immediate consequence of this choice on D is the fact that, when x 1 is stabilized at level x d 1 , the equivalent control is given by D eq = µ 1 (S). Furthermore, control law (4.5) uses information on x 1 and S, and the latter has to be bigger than zero (bigger than S(0) > 0 to have D(t) > 0 for the selected value of χ). Hence, it is necessary to design a control law for S, which can be done using S in (t).

First Control Law -S [1] in

The first control strategy is to be active whenever x 1 (t) is not stabilized in the surface (4.4) or its vicinity. The objective here is to keep S away from zero, otherwise, it would cause D(t) = 0 (due to µ 1 (0) = µ 2 (0) = 0) and no control action on x 1 (t) will be possible. Moreover, it is also necessary to realize the condition imposed in Proposition 4.1 that S(t) ≥ S(0) for all t ≥ 0.

In order to design a control law for S in (t), it is necessary to note that the point S o , defined in Assumption 4.2, corresponds to the value of S at which the system loses its controllability for the variable x 2 with respect to the control input S in (t) (see below the analysis of control for the second phase: at S o the dynamics of ẍ2 is independent in S in , hence, the system is not controllable). In addition, the system is not controllable at all if S = 0. Then the goal is to shift S out of this dangerous region.

To this end, let us introduce three special points: S a,b , S o and S m , which are given by

S o = max a i ∈A i b i ∈B i S o , S a,b = min a i ∈A i b i ∈B i S a,b
hence, S m is selected such as S m > S o . Thus stabilization of S at S m will ensure that, at least at the instant of commutation to regulate x 2 (t), the latter variable is controllable, avoiding the issue remarked above. Hence, in this light, the following result is stated:

S [1] in = 1 D      S max , if σ(S) < 0 0, if σ(S) ≥ 0 (4.7)
where

S max = DS + 2 i=1 ρ i x i µ i (S) + ϵ, ϵ > 0,
the output σ(S) is globally finite-time stabilized provided that D(t) > 0, for all t ≥ 0.

Proof. : Consider a Lyapunov function candidate given by V S (S) = 1 2 (S -S m ) 2 , whose time derivative is given by VS = (S -S m ) Ṡ = (S -S m ) D(S [1] in -S) -

2 i=1 ρ i x i µ i (S)
Then, taking (4.7) into account, one has VS ≤ -ε √ 2V s for some ε ∈ (0, ϵ] dependent on initial conditions, implying the stated result.

T 1 ≤ V 0.5 1 (x 1 (0)) min{x 1 (0), x d 1 }χ √ 2
Consequently, control law (4.7) is well-posed (since D(t) > 0 guarantees that no division by zero will occur). If ϵ > 0, as proven in Proposition 4.2 by means of the Lyapunov function V S , control law (4.7) renders the system with output S -S m globally finite-time stable after a time T 2 > 0 bounded by

T 2 ≤ V 0.5 S (S(0)) ε √ 2
Finally, control laws (4.5) and (4.7) render x 1 (t) = x d 1 and S(t) = S m after a time

T = T 1 + T 2 ≤ 1 √ 2 (min{x 1 (0), x d 1 }χ -1 V 0.5 1 (x 1 (0)) + ε -1 V 0.5 S (S(0))), as claimed.
Second Control Law -S [2] in

The second control has to be active once x 1 is stabilized on the surface (4.4). The objective here is then to stabilize the remaining species x 2 at an arbitrary point x d 2 > 0, which can be done by means of S in , as discussed previously.

Remark 4.2. As x 1 is in sliding motion due to the discontinuous control laws developed previously, an immediate consequence is that after the establishment of the sliding motion for the variable σ 1 , the control D can be considered in its equivalent form [Utkin, 1992] D eq = µ 1 (S), which is assumed in this subsection.

As it can be seen from the model (2.1), if concentration x 2 is in steady-state, i.e., ẋ2 = 0, the equivalent control on D implies that stabilization can be achieved if S is steered to a certain level in which ξ(S) = 0 holds, i.e., S = S a,b . Indeed, if all parameters of the kinetic rates were perfectly known, then the intersection point S a,b would be readily available and the stabilization is easily solved. However, as aforementioned, these parameters are uncertain and therefore the control law designed in this section must provide a robust stabilization of

x 2 .
Before stating the main results of this subsection, let us define the bounds on ξ as follows:

ξ(S) = µ 1 (S) -µ 2 (S), ξ(S) = µ 1 (S) -µ 2 (S),
and also, the bounds on the derivatives of ξ with respect to S are given by

∂ξ ∂S = a 1 b 1 (b 1 + S) 2 - a 2 b 2 (b 2 + S) 2 , ∂ξ ∂S = a 1 b 1 (b 1 + S) 2 - a 2 b 2 (b 2 + S) 2
In this sense, recalling (4.2)-( 4.3) and following the idea of the sub-optimal control (see the Appendix B for the definition of the functional ∆), the control input S [2] in can be selected such as

S [2] in = S + 1 D 2 i=1 ρ i x i μi (S) -k 1 sign[∆(y(t))] (4.8)
where μi (S) =

µ i (S)+µ i (S) 2
is a middle point estimate of the uncertain function µ i (S) that appears in (4.3). Then, by plugging control law (4.8) in the dynamics (4.3), one has that

σ(x 2 ) = x 2 ξ 2 (S) - ∂ξ ∂S 2 i=1 ρ i x i (μ i (S) -µ i (S)) -k 1 sign[∆(y(t))] , = a(t) + b(t)k 1 sign[∆(y(t))] (4.9)
where the terms a(t) and b(t) are given by

a(t) = x 2 ξ 2 (S) - ∂ξ ∂S 2 i=1 ρ i x i (μ i (S) -µ i (S)) b(t) = x 2 ∂ξ ∂S (4.10)
and hence this last expression clearly has the same form as (B.4). Also, in order to eliminate the uncertain term µ i (S), the bounds of functions a(t) and b(t) in (4.10) can be written as

|a(t)| ≤ x 2 ξ 2 (S) + ∂ξ ∂S 2 i=1 ρ i x i µ i (S) -µ i (S) 2 = a max (x 1 , x 2 , S) b min (x 2 , S) = x 2 ∂ξ ∂S (S) ≤ b(t) ≤ x 2 ∂ξ ∂S (S) = b max (x 2 , S) (4.11)
and it is worth noticing that b min (x 2 , S) > 0 while S > S o .

Proposition 4.3. The selection of control law (4.8) renders S in non-negative, for all t ≥ 0,

provided that k 1 ≤ DS + 2 i=1 ρ i x i μi (S).
Proof. The proof is straightforward by noticing that k 1 is the only tunable parameter in (4.8) and, to assure positiveness of S in , the term k 1 sign[∆(y(t))] must not compensate all other (positive) terms. Since sign[∆(y(t))] ∈ [-1, 1] it is obvious that the selection of k 1 such as

k 1 ≤ DS + 2 i=1 ρ i x i μi (S) (4.12)
renders S in ≥ 0 for all t > 0, as claimed.

As discussed previously, since the control S in is multiplied by the term ∂ξ ∂S in ( 4.3), it implies that if S reaches the interval S ≤ S o , the stabilization of surface (4.3) might be no longer possible (actually the system loses its controllability at S o but, since this value is uncertain, it is prudent to keep S > S o ). Clearly, control law (4.8) does not guarantee a priori that S will not reach this interval.

In order to state conditions to overcome such a problem and also to ensure a proper selection of k 1 , a deeper understanding of what happens in such a control phase is discussed in the sequel. Assume that k 1 is properly tuned and that t 2 is the time instant in which this control phase is activated. Since the previous control phase has stabilized S at a level S m > S o , then ẋ2 (t 2 ) > 0 and hence two different scenarios are possible:

1. If x 2 (t 2 ) > x d
2 , then S must increase in order to have ẋ2 < 0;

2. If x 2 (t 2 ) < x d 2 , then S must decrease in order to have ẋ2 > 0;

Obviously, the second case is troublesome since, depending on the selected gain, S might reach the region in which the stabilization is compromised. Hence, the stabilization of surface (4.3) will proceed for x d 2 = x 2 (t 2 ) is assumed. With this assumption, the problem of having S ≤ S o may be alleviated.

In the following, an interval

I = [S, S] × x d 1 × [x d 2 -δ, x d 2 +
δ] for all t ≥ t 2 and some S o < S < S < +∞ and δ > 0 is considered. Due to functioning of the suboptimal control, the dynamics of the closed-loop system will be governed by the equation ÿ

(t) = a(t) -b(t)k 1 .
If the choice of k 1 imposes that a(t) < b(t)k 1 , then ÿ(t) = -r gives the worst-case trajectory estimate. Since in this case y(t 2 ) = 0 and ẏ(t 2 ) > 0, let us define, for (S(t), x 1 (t) x 2 (t)) ∈ I, the minimum value of ÿ(t) as r = a max (x 1 , x 2 , S) -b min (x 2 , S)k 1 for the time interval [t 2 , t 3 ],

where t 3 > t 2 is the instant of time that ẏ(t 3 ) = 0. Omitting all arguments for the sake of readability, this estimate can be expressed as follows:

r = -(x d 2 + δ) ξ(S) 2 + ξ ′ max x d 1 ρ 1 µ 1 (S)-µ 1 (S) 2 + (x d 2 + δ)ρ 2 µ 2 (S)-µ 2 (S) 2 -(x d 2 -δ)ξ ′ min k 1 (4.13)
where ξ ′ min = min I ∂ξ ∂S and ξ ′ max = max I ∂ξ ∂S (this notation will be also used in the following). In this light, let us state the following result: Lemma 4.1. If there exists a compact interval I in which the following inequalities hold:

(a) 1 ξ ′ min ξ 2 (S) + ξ ′ max 2 i=1 ρ i x i µ i (S)-µ i (S) 2 < DS + 2 i=1 ρ i x i μi (S), (b) |ξ(S)| ≤ |x d 2 ξ(Sm)| x d 2 -δ and (c) 1 2r (x d 2 ξ(S m )) 2 ≤ δ for all (S(t), x 1 (t), x 2 (t)) ∈ I, x 1 = x d 1 and x 2 = x d 2 +
δ, then there is a choice of constant gain k 1 in (4.8) that assures both positiveness of S in and the permanence of trajectories of S and x 2 inside the domain of attraction of the origin of (4.3).

Proof. The idea here is to show that, under the conditions presented in this lemma, it is possible to select a constant value of k 1 that satisfies all constraints of non-negativeness and stabilizability of (4.3), during all transients.

First, the right-hand side of inequality (a) is immediate from Proposition 4.4, corresponding directly to the choice of k 1 that yields S in non-negative. The left-hand side, however, relates to the choice of k 1 that stabilizes these dynamics. Indeed, one can easily see that control law (4.8) renders the dynamics (4.3) as

ÿ = a(t) + b(t)k 1 sign[∆(y(t))]
where a(t) and b(t) are given by (4.10). Noticing that this last expression has the same form as (B.4), the stabilization by the suboptimal control is therefore possible

if b min k 1 > |a max |, hence k 1 > |a max | b min
which, recalling (4.11) and the constraint of set I, leads us to

k 1 > 1 ξ ′ min ξ 2 (S) + ξ ′ max 2 i=1 ρ i x i µ i (S) -µ i (S) 2 (4.14)
where, for simplicity,

x 1 = x d 1 and x 2 = x 2 + δ.
Inequalities (b) and (c) are imposed in order to guarantee that (S(t), x 1 (t), x 2 (t)) ∈ I. According to the design of the control on the first step, at the instant of commutation

x d 2 = x 2 (t 2 ), y(t 2 ) = 0 and ẏ2 (t 2 ) = -x d 2 ξ(S m ) > 0.
Hence, estimates on the behaviour of trajectories of S and x 2 can be obtained by solving analytically the equation ÿ(t) = -r (note that r > 0) for t ∈ [t 2 , t 3 ], where the instant t 3 corresponds to the first time instant in which ẏ(t 3 ) = 0.

Consequently, due to properties of the suboptimal control algorithm as shown in Corollary B.1, the following relations hold true:

|y(t)| ≤ y(t 3 ) | ẏ(t)| ≤ y(t 2 ) (4.15)
and hence, by integration of ÿ(t) = -r in the aforementioned time window, one has that ẏ(t) = ẏ(t 2 ) -rt and, as ẏ(t 3 ) = 0 and ẏ(t 2 ) = -x d 2 ξ(S m ), it is immediate that t 3 = -1 r x d 2 ξ(S m ). Integrating once again in the same interval:

y(t) = y(t 2 ) + ẏ(t 2 ) - r 2 t 2
which, if evaluated at t 3 , results in inequality (c), i.e. y(t 3 ) = Proof. This proof is straightforward by the preliminaries given on the suboptimal control (see Appendix B). From (B.10), then (omitting arguments for readability):

T 3 ≤ t 3 + αp |y(t 3 )| where α(S, x 1 , x 2 ) = max I bmax(x 2 ,S)k 1 +amax(x 1 ,x 2 ,S) -r
and p is defined in Subsection 3.2. Hence, profiting the computations performed in the proof of Lemma 4.1, the relation for T 3 above becomes

T 3 ≤ -x d 2 ξ(S m ) r + αp 1 2r x d 2 ξ(S m ) 2
as stated.

Remark 4.3. For the sake of simplicity, the possibility of having time-varying parameters a and b has not been made explicit in the development. However, it can be tackled by the same control architecture in a very natural manner. Recalling formula (4.3), the influence of these S m = 7.5. Also, in accordance with Lemma 4.5, the selection of gain k 1 as

k 1 = 1 4   ξ 2 (S m ) ξ ′ min (S m ) + µ 1 (S m )  
satisfies all constraints on stabilizability, non-negativeness of S in and transient behavior of all trajectories will be satisfied for this simulated scenario.

First, let us analyze the whole performance of the control law concerning the parameters given in Table 4.1, as illustrated in Figure 4.2. As can be seen, both phases successfully stabilizes both species in a finite-time despite of the parametric uncertainties. Immediately after this step, the second phase is initiated and successfully stabilizes x 2 at x 2 (t 2 ) = 1.2. Also, one can see that the control architecture successfully stabilized S = S a,b , even though S a,b belongs to an uncertain interval.

This fact shows the usefulness of this methodology since it allows a robust solution to the stated stabilization problem.

In Figure 4.4, one can see the time evolution of the control inputs S in (t) and D(t). Clearly, one can see the high frequency switching on both control laws, which is natural to the slidingmode control. This phenomenon is called ideal sliding-mode, where the actuators are supposed to allow this high frequency. In order to investigate this control architecture in a numerical experiment that is closer to the real-world, the same scenario will be simulated considering dynamics on the actuators: a first-order filter and a sampler will be introduced to the output of the proposed controllers.

This means that the pumping rates are not perfect and can only be changed at sampling times τ ≈ 30 seconds.

As it can be seen in Figure 4.5, even under these new circumstances, the system can still be stabilized at the desired positions. However, a larger discrepancy is observed for the variable x 2 (t). Figure 4.6 illustrates the control inputs after the actuator dynamics. Clearly, it still oscillates but at much smaller frequencies. It is worth noticing that this is possible thanks to the fact that the actuator is much faster than the dynamics of the controlled system.

Control of the COSY model

In this section, particular attention is drawn to the control problem concerning the COSY model (2.4). For all simulations and numerical results throughout this section, the model parameters will be selected as the ones given by [START_REF] Mauri | Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs[END_REF] and, arbitrarily (but still realistically), k h = 0.2 and m = 1.

For simplicity, the (multiplying) control inputs will be denoted u 1 = D g G in and u 2 = D a A in . Note that these inputs are constrained such as

u 1 ≥ 0, u 2 ≥ 0. (4.16)

Ensuring coexistence through state feedback

Although the ultimate objective is to stabilize the production of H, the first challenge tackled in this section is how to ensure the coexistence of both producer and cleaner strains.

As it will be shown, this ultimate objective will be attained by using the same control law proposed in this subsection.

First, a control law for the dilution rate D will be proposed. As clearly seen from ( 2.4), controlling the concentration of B p allows one to control the concentration of H. Defining an error variable σ p as

σ p = B p -B ⋆ p (4.17)
where B ⋆ p is a desired concentration level for B p , the control law for the dilution rate is selected as follows:

D = max{D min , (1 -k h )µ p (G, A) -k + φ(σ p )} (4.18)
where D min > 0 is a (predefined) minimum value for the dilution rate, φ(σ p ) is any function satisfying φ(σ p )σ p > 0 for σ p ̸ = 0 and µ p (G, A) is as given in (2.5). As this control law depends on the concentration of the substrates, before stating the result concerning stability of σ p , several conditions (in terms of G and A) will be imposed.

First, as readily seen by computing the dynamics of σ p , the saturation D = D min does not ensure the stability of the origin of (4.17). To remedy this issue, the following constraint is needed:

(

1 -k h )µ p (G ⋆ , A ⋆ ) ≥ k + D min
which basically imposes that, by selecting proper values of required levels G ⋆ and A ⋆ for G and A, respectively, the saturation shown in (4.18) will not be reached at the desired equilibrium.

Ḃc = h(G ⋆ , A ⋆ )B c where h(G, A) = µ c (G, A) -(1 -k h )µ p (G, A).
Therefore, in order to keep the value of B c constant in a positive value in the equilibrium, the choice of G ⋆ and A ⋆ also needs to satisfy Consider two other error variables given by The idea is then to drive these variables to zero. Since D is already fixed (i.e., selected as (4.18) aiming at the stabilization of B p ), this can only be done through the remaining inputs u 1 and u 2 (i.e., through the actuation of G in and A in , respectively).

h(G ⋆ , A ⋆ ) = 0.
σ g = G -G ⋆ , σ a = A -A ⋆ .
In this light, consider the following control laws for u 1 and u 2 :

u 1 = max 0, r 1 (G, A)B p + r 2 (G, A)B c + DG -φ(σ g ) (4.19) u 2 = max 0, r 3 (G, A)B p + r 4 (G, A)B c + DA -φ(σ a ) (4.20)
Proposition 4.4. Consider system (2.4) with ouputs σ g and σ a . If u 1 and u 2 can be chosen arbitrarily large, then the surfaces σ g = 0 and σ a = 0 are asymptotically stable under control laws (4.19), provided that r 1 (G, A) < l and r 2 (G, A) < l.

Proof. The restrictions r 1 (G, A) and r 2 (G, A) imply that the terms proportional to B p and B c in the dynamics of A are sign-defined (i.e., the extra acetate input discussed in Section 4.2 will not disrupt controllability of A through the positive control). Under this condition, the dynamics of G and A are similar, therefore this proof will consider only the former. For simplicity, the dynamics of G is rewritten as

Ġ = -f 1 -DG + u 1 ,
where

f 1 = f 1 (G, A, B p , B c ) = r 1 (G, A)B p + r 2 (G, A)B c ≥ 0. Considering a candidate Lya- punov function V g = 1 2 σ 2 g
, whose time derivative is computed as follows:

Vg = σ g σg = σ g (u 1 -f 1 -DG) .
saturation in (4.18) is selected as D min = 0.05 and, for simplicity, the functions φ(•) where simply selected as φ(s) = s.

The desired steady-states for B p and G are, respectively, B ⋆ p = 0.5 [gL -1 ] and

G ⋆ = 3 [gL -1 ]. According to

Part conclusion

In this chapter, the control problem for systems (2.1) and (2.4) has been investigated. As it was discussed, controllability of such systems is a complicated task, due to their non-negative, nonlinear nature and bounded control inputs.

First, for the competition model (2.1), a gentle motivation for the need of a robust approach is presented. In this light, a robust control architecture composed of two phases (each one with a control objective concerning the input S in ) is proposed, allowing the explicit consideration of parametric uncertainty in the reaction rates, and its stabilizing features are proven. These control laws are based on switching controllers, guaranteeing both non-negativity of the inputs and controllability of the system. A numerical example illustrates the methodology.

Considering the COSY model (2.4), the control problem becomes more complicated. Indeed, the excretion of acetate by the producer strain can be considered as an extra input, which further restrain the controllability of the system. Therefore, in order to stabilize H (which can be done simply by stabilizing the concentration of the producer strain B p ), a control law using the dilution rate was proposed. However, to guarantee controllability and also the permanence of the cleaner strain, several constraints concerning the concentrations of glucose and acetate were imposed. Therefore, two more control laws using the inflow concentration of each substrate were proposed, allowing the values of glucose and acetate to be stabilized at proper values that respect such constraints. A numerical experiment illustrate the efficacy of this methodology.

Part 3:

Robust Output Feedback Model Predictive Control

Chapter MPC using interval estimators

Introduction Problem statement

Consider a general discrete-time system given by

x k+1 = f (x k , u k , w k , θ k ) y k = Cx k + v k (5.1)
where x k ∈ R n is the state vector, u k ∈ R m is the input vector and y k ∈ R p is the available measurement vector. The measurement matrix C ∈ R p×n is known and constant and f :

R 2n+m+r + → R n . The signals w ∈ ℓ n ∞ and v ∈ ℓ p ∞ are
, respectively, process and measurement noise. The time-varying signal θ k ∈ Θ ⊂ R r represent some parametric dependence (such as scheduling or uncertain parameters). Assumption 5.1. Initial conditions of (5.1) are bounded such as x 0 ≤ x 0 ≤ x 0 , for some known x 0 , x 0 ∈ R n . Furthermore, the additive perturbations Roughly speaking, Problem 5.1 requires robust constraint satisfaction, which means that the constraints imposed in both state and control input must be respected even if the system is plagued by disturbances and uncertainties. The main rationale of the solutions proposed for such a problem relies on interval estimators. These estimators, while having a simple structure, provide information on the set-membership of the states (i.e., in the form of intervals), with guaranteed convergence and low computational complexity. Their design is achieved by using the bounds on the uncertainties, and under some conditions on cooperativity (or non-negativeness).

w k ∈ [w k , w k ] and v k ∈ [v k , v k ] for all k ∈ Z + ,
Roughly speaking, an IO uses the available measurement to generate, at every instant k, signals x k , x k ∈ R n that satisfy the following relation:

x k ≤ x k ≤ x k , ∀k ∈ Z + .
(5.2)

On the other hand, an interval predictor (IP, also called a framer [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF]) can be seen as an open-loop estimator, since it does not depend on the measurements y k . For this reason, it can be used to predict an envelope, where all the trajectories of the considered system are enclosed. In this light, the IP generates signals z t , z t ∈ R n such that .3) where N ∈ Z + is the finite prediction horizon. Therefore, the main idea is to use this information to check the fulfillment of constraints in the MPC, since

z t ≤ x t ≤ z t , ∀t ∈ {k, k + 1, . . . , k + N }. ( 5
[z k , z k ] ⊂ X ⇒ x k ⊂ X,
and, obviously, the same property can be checked from (5.2).

Preliminaries Preliminaries on MPC

A standard MPC algorithm solves a constrained optimization problem, by predicting the behaviour of the system over a prediction horizon N (see Figure 5.1 below), to determine a control input sequence S N = {s 0 , . . . , s N -1 } that minimizes a cost function V N (assuming, for simplicity, that the whole state vector is available as exact measurements). This algorithm is summarized as the following minimization problem:

S N := arg min S N V N (x k,0 , . . . , x k,N , S N ) (5.4)
subject to the following constraints (5.5d) where X f ⊂ X is the terminal set.

x k,0 = x k , (5.5a) x k,i+1 = f (x k,i , s i , 0, 0), (5.5b) x k,i+1 ⊂ X, s i ⊂ U, (5.5c) x k,N ∈ X f ,
Constraint (5.5a) imposes the initialization of the prediction with the measured x k , constraint (5.5b) states that the predicted states are obtained using the available dynamical model, constraint (5.5c) requires that the prediction respects the constraints on state and control, and (5.5d) is a stabilizing constraint, imposing the endpoint of the prediction to be in a terminal set (the features of this set will be discussed in the sequel).

The cost function V N is defined as follows

V N = V f (x k,N ) + N -1 i=0 ℓ(x k,i+1 , s i ).
where V f (x) is the terminal cost and ℓ(x, s) is the stage cost.

Applying the philosophy of the receding-horizon control, a sequence S N is obtained at every decision instant k and only its first element, s 0 , is applied to the system from the instant k to k + 1. This algorithm is then applied for all subsequent time instants.

Concerning the stability of an MPC algorithm, an important requirement is recursive feasibility:

Definition 5.1. The MPC algorithm is recursive feasible iff, for all initially feasible x 0 , it remains feasible for all subsequent times.

Indeed, the stability of an MPC controller can be evoked by the classic results reviewed in [START_REF] Mayne | Constrained model predictive control: stability and optimality[END_REF] to establish stabilizing ingredients for recursive feasibility. These ingredients are the terminal set X f , the terminal cost V f and the terminal controller κ f (x), which satisfy the following definition:

Definition 5.2. The stabilizing ingredients are such that the following axioms are verified:

1. X f ⊂ X, closed and 0 ∈ X f : the state constraint is satisfied in X f ; 2. κ f (x) ∈ U, ∀x ∈ X f : the control constraint is satisfied in X f ; 3. f (x, κ f (x)) ∈ X f , ∀x ∈ X f : X f is positively invariant under κ f (x); 4. [V f + ℓ](x, κ f (x)) ≤ 0, ∀x ∈ X f : V f is a local Lyapunov function.

Preliminaries on interval arithmetic

The design of the interval estimators presented in this chapter will require some results on interval arithmetic. To this end, we evoke the following:

Lemma 5.1. [Efimov et al., 2013] Let x ∈ R n be a vector variable, x ≤ x ≤ x for some

x, x ∈ R n . Then, (1) if A ∈ R m×n is a constant matrix, then

A + x -A -x ≤ Ax ≤ A + x -A -x.
(5.6)

(2) if A ∈ R m×n is a matrix variable and A ≤ A ≤ A for some A, A ∈ R m×n , then

A + x + -A + x --A -x + + A -x -≤ Ax ≤ A + x + -A + x --A -x + + A -x -.
(5.7) Lemma 5.2. [START_REF] Efimov | Design of interval observers for uncertain dynamical systems[END_REF] For A ∈ R n×n + , the system

x k+1 = Ax k + ω k , ω : Z + → R n + , ω ∈ ℓ n ∞ , k ∈ Z +
has a non-negative solution x k ∈ R n + for all k ∈ Z + provided that x 0 ≥ 0. Definition 5.3. A system satisfying Lemma 5.2 is called cooperative or monotone. Lemma 5.3. [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF] A matrix A ∈ R n×n + is Schur stable iff there exists a diagonal matrix P ∈ R n×n , P ≻ 0, such that A ⊤ P A -P ≺ 0.

OF-MPC for LTI systems

In this section, a solution for Problem 5.1 will be investigated for the case in which system (5.1) is an LTI one [START_REF] Souza | Robust output feedback model predictive control for constrained linear systems via interval observers[END_REF], i.e., (5.8) in which the (constant) matrices A ∈ R n×n and B ∈ R n×m are known.

x k+1 = Ax k + Bu k + w k

Design of interval estimators

Interval observer

First, by using output injection of the available measurement, system (5.8) can be rewritten as

x k+1 = (A -LC)x k + Bu k + Ly k -Lv k + w k (5.9)
for any L ∈ R n×p . Under Assumption 5.1, an IO for such a system can be constructed by introducing a gain L o to be determined (which replaces L above), and using (5.6) to replace the uncertain terms (i.e., w k -Lv k ) by their interval bounds. The resulting IO [Efimov et al., 2013] is given as follows:

x k+1 = D o x k + Bu k + L o y k -L + o v k + L - o v k + w k x k+1 = D o x k + Bu k + L o y k -L + o v k + L - o v k + w k (5.10)
where D o = A -L o C, and having x 0 , x 0 as specified in Assumption 5.1. Note that the precision of (5.10) can be evaluated by the width of its interval, i.e., δx k = x k -x k , whose dynamics are given by: (5.11) where δw k = w k -w k and δv k = v k -v k determine the uncertainty size of the state and the output disturbances, respectively. Since δx 0 = x 0 -x 0 , the solution of (5.11) is obtained by

δx k+1 = D o δx k + δw k + |L o |δv k ,
δx k = D k o (x 0 -x 0 ) + k-1 i=0 D k-1-i o (δw i + |L o |δv i )
for all k ∈ Z + . Hence, the values of δx k are completely determined by the choice of the gain L o and the uncertainty levels given in Assumption 5.1. Our first result shows a procedure to compute this gain.

Theorem 5.1. Let Assumption 5.1 be satisfied. If the following inequalities are verified for a scalar ρ > 0, a diagonal matrix P ∈ R n×n and matrices W + , W -∈ R n×p + :

min ρ,P,W + ,W -ρ P > 0, P A -(W + -W -)C ≥ 0,      P -I n 0 A ⊤ P -C ⊤ (W + -W -) ⊤ ⋆ ρI n E ⊤ ⋆ ⋆ P      ⪰ 0, E = P W + + W -,
then system (5.10) with gains L + o = P -1 W + and L - o = P -1 W -is an IO for system (5.8), i.e., relation (5.2) holds and δx k ∈ ℓ n ∞ . Furthermore, the transfer (δw k , δv k ) -→ δx has a gain lesser than √ ρ.

Proof. Let the estimation errors be given by e k = x k -x k and e k = x k -x k , whose increments have the following form:

e k+1 = D o e k + w k -w k + L o v k -L + o v k + L - o v k , e k+1 = D o e k + w k -w k -L o v k -L + o v k + L - o v k .
(5.12) Then, under Assumption 5.1 and Lemma 5.1, all exogenous inputs (i.e., the independent right-most terms) in (5.12) are non-negative. Consequently, if D o is also non-negative, we

have that e k , e k ≥ 0 for all k ∈ Z + under Lemma 5.2. This requirement is imposed by the first inequality on this theorem and implies relation (5.2).

Note that, under Assumption 1, the asymptotic stability of (5.11) and ( 5.12) sums up to the Schur stability of D o . In this light, let ν k = vec(δw k , δv k ) and consider a candidate Lyapunov function V k = δx ⊤ k P δx k , whose increments for (5.11) are given by

V k+1 -V k =   δx k ν k   ⊤   D ⊤ o P D o -P D ⊤ o E E ⊤ D o E ⊤ P -1 E   Π   δx k ν k   .
(5.13) Thus, a sufficient condition for Schur stability of D o is Π ⪯ 0. Now, by introducing a performance index given by J

= ∞ k=0 δx ⊤ k δx k -γ 2 ν ⊤ k ν k + (V k+1 -V k )
, the desired ℓ ∞ performance is achieved by rendering J < 0, while minimizing γ. Taking (5.13) into account, this condition is equivalent to

  D ⊤ o P D o -P + I n D ⊤ o E E ⊤ D o E ⊤ P -1 E -γ 2 I n   ⪯ 0.
Since P is diagonal, we can decompose the relation above as follows:

  D ⊤ o P E ⊤   P -1 P D o E -   P -I n 0 0 γ 2 I n   ⪯ 0
and hence, by applying the Schur complement in the inequality above, and by introducing the variables ρ = γ 2 , W + = P L + o and W -= P L - o , the inequalities given in this theorem are obtained, proving the claim.

Interval predictor

To avoid confusion with the IO developed in the previous section, z k and z k will denote the predicted upper and lower bounds for x k , respectively, for all k ∈ Z + . Also, the predictor gain (which replaces L in (5.9)) will be denoted as L p .

By definition, L

p = L + p -L - p for some L + p , L - p ∈ R n×p + .
Then, under (5.2), Assumption 5.2 and Lemma 5.1, we can state that

L + p Cx k -L - p Cx k ≤ L p Cx k ≤ L + p Cx k -L - p Cx k
which allows us to substitute the terms which are unavailable for prediction in (5.9) (i.e.,

Ly k -Lv k + w k = LCx k + w k ),
with their respective bounds. The following IP is then obtained: (5.14) where D p = A -L p C. As readily seen, (5.14) is composed only by known terms (under Assumption 5.1 and assuming that u k is to be computed in the control algorithm).

z k+1 = D p z k + Bu k + L + p Cz k -L - p Cz k + w k , z k+1 = D p z k + Bu k + L + p Cz k -L - p Cz k + w k ,
Introducing a change of coordinates describing the center of the interval,

z ⋆ k = z k +z k 2
, and the width of the interval, δz k = z k -z k , one gets the following equivalent representation of dynamics of (5.14): (5.15)

z ⋆ k+1 = Az ⋆ k + Bu k + w ⋆ k ,
δz k+1 = (A + 2L - p C)δz k + δw k (5.16)
where, similarly,

w ⋆ k = w k +w k 2
and δw k = w k -w k . Two interesting features are to be noted:

(i) the dynamics of the center of the predicted interval z ⋆ k is independent of L p , but it is controlled by u k and, hence, (5.15) can be used in MPC algorithm, and (ii) the dynamics of the interval width δz k is governed by the (known) interval width of the state disturbance δw k , and its stability can be assured by a proper choice of the gain L p . The following theorem offers a procedure to compute such a gain. Theorem 5.2. Let assumptions 5.1-5.2 be satisfied, and there exist a diagonal matrix P ∈ R n×n , matrices Q, Γ ∈ R n×n and U -, U + ∈ R n×p + such that the following linear matrix inequalities are verified:

P A -U + C + U -C ≥ 0, (5.17)      P -Q 0 A ⊤ P + 2C ⊤ U -⊤ 0 Γ P P A + 2U -C P P      ⪰ 0 P ≻ 0, Q ≻ 0, Γ ≻ 0.
Then, system (5.14) with gains L - p = P -1 U -, L + p = P -1 U + and z 0 = x 0 , z 0 = x 0 is an IP for system (5.8), i.e., the relation (5.3) is satisfied and δz ∈ ℓ n ∞ .

Proof. Assume that the gain L p is selected such that the matrix D p is non-negative. Then, the realization of relation ( 5.3) for system (5.8), (5.14) (under substitution of z k , z k instead of

x k , x k , respectively) can be proven following the same arguments as in Theorem 5.1. Recalling

that L p = L + p -L - p with L + , L -∈ R n×p + , then A -L p C = A -L + p C + L - p C ≥ 0
is the condition to verify under Assumption 5.2. The stability of (5.16) implies the stability of the IP (5.14). In this sense, note that

0 ≤ A -L + p C + L - p C ≤ A + 2L - p C since L - p C ≥ 0 and L + p C ≥ 0 under Assumption 5.2. Then A + 2L - p C ∈ R n×n +
and according to Lemma 5.3, we can select a diagonal matrix P ∈ R n×n considering V k = δz ⊤ k P δz k as a Lyapunov function candidate (with P ≻ 0), whose increment takes, for any Q ∈ R n×n and Γ ∈ R n×n , the following form:

V k+1 -V k =   δz k δw k   ⊤ Σ   D ⊤ P D -P + Q D ⊤ P P D P -Γ     δz k δw k   -δz ⊤ k Qδz k + δw ⊤ k Γδw k where D = A + 2L -C. If Q ≻ 0, Γ ≻
0 and Σ ⪯ 0, then the stability conditions are fulfilled and the system (5.16) is ISS with respect to the input δw k :

V k+1 -V k ≤ -δz ⊤ k Qδz k + δw ⊤ k Γδw k .
The first LMI in this theorem relates to the condition A -(L + -L -)C ≥ 0 and is obtained by defining U -= P L -∈ R n×p + and U + = P L + ∈ R n×p + . The second condition concerns the Schur stability of Σ. Decomposing the matrix Σ, similarly as in Theorem 5.1, and taking its Schur complement, the second LMI is readily obtained by recalling that P D = P A + 2P L -C. This proves the theorem.

Remark 5.1. Similarly to Proposition 1, an additional optimization problem can be posed for the maximization of Q and minimization of Γ while solving the LMIs of Theorem 5.2,aiming an optimal accuracy in the interval prediction.

Control design

Before discussing the MPC algorithm to be proposed, we stress the stabilization of the IP through a (static) state feedback. Note that IP (5.14) can be stabilized by designing a feedback for its completely known center dynamics (5.15)

z ⋆ k+1 = (A + BK f )z ⋆ k + (I n + BS f )w ⋆ k , (5.18)
using the following control law:

u k = S f w ⋆ k + K f z ⋆ k (5.19)
with the gains K f ∈ R m×n and S f ∈ R m×n guaranteeing input-to-state stability of the closedloop system and minimizing the influence of w ⋆ k , respectively. These gains can be selected as follows:

Theorem 5.3. Let S f = ΣP -1 and K f = ΥP -1 , where P ∈ R n×n and Σ, Υ ∈ R m×n are solutions of a linear optimization problem:

max Q,Γ 1 ,Γ 2 ∈R n×n Q -Γ 1 -Γ 2 , Q = Q ⊤ ≻ 0, Γ 1 = Γ ⊤ 1 ≻ 0, Γ 2 = Γ ⊤ 2 ≻ 0, P = P ⊤ ≻ 0, Π ⪰ 0 Π =         P -Q 0 0 P A ⊤ + Υ ⊤ B ⊤ 0 Γ 1 0 P + Σ ⊤ B ⊤ 0 0 Γ 2 P AP + BΥ P + BΣ P P        
.

Then the center dynamics (5.15) with the control (5.19) is ISS (from the input (w ⋆ k ) to the state z ⋆ k ) with the optimal attenuation of the disturbances and with an ISS-Lyapunov function (5.19) into the (5.15) yields the closed-loop dynamics (5.18). Consider the increment of the candidate Lyapunov function V (z ⋆ k ):

V (z ⋆ k ) = z ⋆ k ⊤ P -1 z ⋆ k . Proof. Substituting the control
V (z ⋆ k+1 ) -V (z ⋆ k ) =   z ⋆ k w ⋆ k   ⊤ Π   Ã⊤ P -1 à -P -1 + Q Ã⊤ P -1 D D⊤ P -1 à D⊤ P -1 D -Γ1     z ⋆ k w ⋆ k   -z ⋆ k ⊤ Qz ⋆ k + w ⋆ k ⊤ Γ1 w ⋆ k , where à = A + BK f , D = I n + BS f and for any Q, Γ1 ∈ R n×n . If Π ⪯ 0, then V (z ⋆ k+1 ) -V (z ⋆ k ) ≤ -αV (z ⋆ k ) + w ⋆ k ⊤ Γ1 w ⋆ k (5.20)
meaning that the system is ISS, provided that Q ⪰ αP -1 and Γ1 ≻ 0, where α > 0 always

exists if Q ≻ 0. Note that Π = -   P -1 -Q 0 0 Γ1   +   Ã⊤ P -1 D⊤ P -1   P P -1 Ã P -1 D
is such that, by applying the Schur complement, the LMIs given in the formulation of this

proposition are verified for Σ = S f P , Υ = K f P , Q = P -1 QP -1 , Γ1 = P -1 Γ 1 P -1 .
Remark 5.2. Additional constraint can be imposed on P , in order to guarantee some optimal performance for the terminal cost under the control (5.19).

Remark 5.3. If w ⋆ k = 0 for k ∈ Z + then an obvious choice is S f = 0, which can be imposed as a constraint on the optimization problem given by Proposition 5.3.

As a consequence of property (5.20), the following ellipsoid is a positively invariant set for (5.8) under control (5.19):

X = z ⋆ ∈ R n : z ⋆⊤ P -1 z ⋆ ≤ α -1 sup k≥0 w ⋆ k ⊤ Γ1 w ⋆ k , ( 5.21)
for the α obtained in Theorem 5.3. The following assumption, which is conventionally imposed in MPC [START_REF] Mayne | Constrained model predictive control: stability and optimality[END_REF], is introduced:

Assumption 5.3. Let the terminal set X f ⊆ X ⊆ X and u k ∈ U for all z ⋆ k ∈ X f in (5.19).
Assumption 5.3 may be relaxed by further constraining the conditions given in Theorem 5.3, as suggests the following corollary:

Corollary 5.1. Let there exist symmetric and positive definite matrices U ∈ R m×m and

Z ∈ R 2n×2n such that U = {u ∈ R m : u ⊤ Uu ≤ 1} and W k ∈ {W ∈ R 2n : W ⊤ ZW ≤ 1}
, and the conditions of Theorem 5.3 be satisfied with additional inequalities:

η ακ Γ ≤ min{κ -1 Z, P }, P ≥ κZ -1 ,      η 3 P 0 Υ ⊤ 0 κ 3 P Σ ⊤ ⋆ ⋆ U -1      ≥ 0 (5.22)
for some constants η > 0 and κ > 0, then control (5.19)

satisfies the constraint u k ∈ U for all z ⋆ k ∈ X f .
Proof. First of all, the imposed inequalities on Z and Γ imply:

Z ≥ κP -1 , η α Γ 1 ≤ Z, η α Γ 1 ≤ κP ≤ P ZP ⇒ η α Γ1 = η α P -1 Γ 1 P -1 ≤ Z, consequently, ηz ⋆ k ⊤ P -1 z ⋆ k ≤ η α w ⋆ k ⊤ Γ1 w ⋆ k ≤ w ⋆ k ⊤ Zw ⋆ k ≤ 1 for z ⋆ k ∈ X.
Next, note that the condition u k ∈ U takes the form:

  z ⋆ k w ⋆ k   ⊤   K f S f   U   K f S f   ⊤   z ⋆ k w ⋆ k   ≤ 1.
If z ⋆ k ∈ X, then the ellipsoid (5.21) can be evoked to rewrite the previous inequality as

  K f S f   U   K f S f   ⊤ ≤ 1 2   ηP -1 0 0 Z   .
Applying Schur complement the latter property is equivalent to

     η 2 P -1 0 K ⊤ f 0 1 2 Z S ⊤ f K f S f U      ≥ 0,
then multiplying this inequality from both sides by diag{P, P, I m } and taking into account that P ZP ≥ κP we get the LMI of the corollary.

An extension to linear time-delayed systems

The idea developed in this section was also applied to linear time-delayed systems [START_REF] Souza | Robust output feedback mpc of time-delayed systems using interval observers[END_REF], whose dynamics are given by the following retarded difference equation:

x k+1 = A 0 x k + A 1 x k-h + Bu k + w k ,
where x k-h ∈ R n is the delayed state, with a known and fixed delay h > 0.

For brevity of exposition, these results are not reported in this thesis. In such a case, each interval estimator has two gains (for illustration, L 0 and L 1 ), one relating to the actual state, and another relating to the delayed one. Due to this new term, the stability analysis must be then carried out using the Lyapunov-Krasovskii framework [Fridman, 2014].

Interestingly, the proposed methodology easily encompasses delays in states, inputs and outputs, features that were not entirely reported on the robust MPC literature. The challenge in this scenario, however, remains in the fact the two gains must be designed to render both A 0 -L 0 C and A 1 -L 1 C non-negative (in order to ensure non-negativity of the estimation errors), instead of only one as discussed in this section.

Numerical illustration

Consider the (linearized) continuous stirred tank reactor (CSTR), where an exothermic and irreversible reaction S → P occurs [START_REF] Henson | [END_REF]]. The model is composed of two states, the concentration of reactant C S and the temperature of the reactor T , and one controlled input, the coolant stream T C . The resulting linear system is as follows:

x k+1 =   0.745 -0.002 5.610 0.780

  x k +   5.6 × 10 -6 0.464   u k + w k , y k = 0 1 x k + v k .
( 5.23) where the disturbances are considered to be enclosed in W = [-0.02, 0.02]×[-0.2, 0.2] and V = [-0.2, 0.2]. The initial conditions are assumed to be bounded such as

x 0 ∈ [-1.1, -1] × [6, 7].
Solving the conditions given in this section, the following gains are obtained for the IO and IP:

L o =   -0.02 0.78   , L p =   -0.02 0.39  
while for the control (considering S f = 0 due to the symmetry of the disturbances bounds), the feedback gain obtained was K f = -6.99 -0.50 . The IO/IP are initialized according to the initial conditions given above.

In the following, Figures 5.2 and 5.3 depict, respectively the evolution of the states of the real system and the pair IO/IP, and the evolution of the computed control signal. As it can be seen, the stabilization of the interval center is successful. Since the pair IP/IO has stable widths and thanks to relations (5.2) and ( 5.3), the real trajetories of (5.23) is also stabilized in a vicinity of the origin (which evidences the practical ISS).

However, in Figure 5.2, one can see that the IP is much more conservative than the IO. by their interval bounds leads to the following IO: (5.26) where Proof. First, by applying (5.7) in λ i x k , one obtains

x k+1 =D o x k + ∆A + x + k + ∆A -x - k + B 0 u k + ∆Bu + k + L o y k -L + o v k + L - o v k + w k x k+1 =D o x k -∆A + x - k -∆A -x + k + B 0 u k -∆Bu - k + L o y k -L + o v k + L - o v k + w k
D o = A 0 -L o C, L o ∈ R n×p
λ + i x + k -λ + i x - k -λ - i x + k + λ - i x - k ≤ λ i x k ≤ λ + i x + k -λ + i x - k -λ - i x + k + λ - i x - k Since λ i ∈ [0, 1], then λ + i = 1 and λ + i = λ - i = λ - i = 0, by definition. Under this, the relation -x - k ≤ λ i (θ k )x k ≤ x + k holds.
Then, since the vertices of the polytopic system are known, relation (5.6) can be applied to obtain the following inequalities

-∆A + x - k -∆A -x + k ≤ N i=1 λ i (θ k )∆A i x k ≤ ∆A + x + k + ∆A -x - k .
(5.27)

Finally, the same idea applies to the term proportional to u k , leading to the following relation: where

-∆B i u - k ≤ -(∆B i u k ) -≤ λ i ∆B i u k ≤ (∆B i u k ) + ≤ ∆B i u + k ( 5
r 1,1 = ∆A + x + k + ∆A -x - k - N i=1 λ i (θ k )∆A i x k + ∆Bu + k - N i=1 λ i (θ k )∆B i u k , r 1,2 = L o v k -L + o v k + L - o v k + w k -w k , r 2,1 = N i=1 λ i (θ k )∆A i x k -(-∆A + x - k -∆A -x + k ) + N i=1 λ(θ k )∆B i u k + ∆Bu - k , r 2,2 = L + o v k -L - o v k -L o v k + w k -w k .
(5.29)

From assumptions 5.1 and 5.4 and relations (5.27)-(5.28), the quantities given in (5.29) are positive. Hence, if D o is non-negative, then e k , e k > 0 for all k > 0 under Lemma 5.2, and thus satisfying relation (5.3). Now, it is needed to derive stability conditions for IO (5.26). First, let us denote χ k = vec(x k , x k ) and rewrite (5.26) as 5.30) where

χ k+1 = A 0 -Lo C 1 χ k + A + χ + k + A -χ - k + δ k (
A 0 = diag (A 0 , A 0 ) ∈ R 2n×2n , Lo = diag(L o , L o ) ∈ R 2n×2p , C 1 = diag (C, C) ∈ R 2p×2n , δ k = vec(δ k , δ k ),
and

A + =   ∆A + 0 -∆A -0   , A -=   0 ∆A - 0 -∆A +   , δ k = B 0 u k + ∆Bu + k + L o y k -L + o v k + L - o v k + w k , δ k = B 0 u k -∆Bu - k + L o y k -L + o v k + L - o v k + w k .
For ease of notation in the sequel, let us denote Ũ = diag (U, U ) and P = diag (P, P ) for some decision variables P ∈ R n×n and U ∈ R n×p . A gain L o that stabilizes (5.26) and satisfies the restrictions of Lemma 4 can be computed by verifying the following conditions:

Theorem 5.4. Let assumptions 5.1-5.4 be satisfied. If there exist diagonal matrices P , Q 1 ,

Q 2 , Q 3 , Ω + , Ω -, Ψ ∈ R 2n×2n , matrices Γ ∈ R 2n×2n and Ũ ∈ R 2n×p
, such that the following LMIs are verified:

P A 0 -Ũ C 1 ≥ 0            P -Q 1 -Ω + -Ω -0 A ⊤ 0 P -C ⊤ 1 Ũ ⊤ ⋆ -Q 2 -Ψ 0 A ⊤ + P ⋆ ⋆ -Q 3 0 A ⊤ ⋆ P            ⪰ 0 P > 0, Γ ≻ 0, Q 1 , Q 2 , Q 3 , Ω + , Ω -≥ 0, Q 1 + min{Q 2 , Q 3 } + 2 min{Ω + , Ω -} > 0
(5.31)

then system (5.26) with a gain L o = P -1 U is an IO for system (5.24), i.e., relation (5.2) is satisfied and, in addition, χ ∈ ℓ 2n ∞ provided that δ ∈ ℓ 2n ∞ .

Proof. Let us consider a Lyapunov function candidate V k = χ ⊤ k P χ k , whose increments are given by

V k+1 -V k =         χ k χ + k χ - k δ k         ⊤ Σ         D ⊤ o P D o -P + Q 1 D ⊤ o P A + + Ω + D ⊤ o P A -+ Ω -D ⊤ o P ⋆ A ⊤ + P A + + Q 2 A ⊤ + P A -+ Ψ A ⊤ + P ⋆ ⋆ A ⊤ -P A -+ Q 3 A ⊤ - P ⋆ ⋆ ⋆ P -Γ                 χ k χ + k χ - k δ k         -χ k ⊤ Q 1 χ k -χ + k ⊤ Q 2 χ + k -χ - k ⊤ Q 3 χ - k -2χ ⊤ k Ω + χ + k -2χ ⊤ k Ω -χ - k -2χ +⊤ k Ψχ - k + δ ⊤ k Γδ k . If Q = Q 1 + min{Q 2 , Q 3 } + 2 min{Ω +
, Ω -} ≻ 0, Γ ≻ 0 and Σ ⪯ 0 and provided that δ ∈ ℓ 2n ∞ , then the stated stability conditions are fulfilled and system (5.26) is input-to-state stable (ISS) with respect to the input δ k (the diagonal matrix Ψ can be sign indefinite since χ +⊤ k Ψχ - k = 0 by definition):

V k+1 -V k = -χ ⊤ k Qχ k + δ ⊤ k Γδ k
Hence, it is needed to show that the above stability conditions can be formulated as LMIs.

First, since P > 0, then Σ can be decomposed as

Σ = D ⊤ o P A ⊤ + P A ⊤ - P P P -1         P D o P A + P A - P         ⊤ -         P -Q 1 -Ω + -Ω -0 -Ω + -Q 2 -Ψ 0 -Ω - -Ψ -Q 3 0 0 0 0 Γ        
then, by applying the Schur complement, the condition Σ ⪯ 0 can be equivalently written as

           P -Q 1 -Ω + -Ω -0 D ⊤ o P ⋆ -Q 2 -Ψ 0 A ⊤ ⋆ P            ⪰ 0 (5.32) Denote U = P L o ∈ R n×p + . Recalling that P D o = P A 0 -P L o C, then inequality (5.32) becomes linear in P , Q 1 , Q 2 , Q 3 , Γ, Ω + , Ω -, Ψ and U :            P -Q 1 -Ω + -Ω -0 A ⊤ 0 P -C ⊤ 1 Ũ ⊤ ⋆ -Q 2 -Ψ 0 A ⊤ + P ⋆ ⋆ -Q 3 0 A ⊤
Then, according to relations (5.27), (5.28) and (5.33), all inputs (the terms independent on ϵ and ϵ) in (5.35) are non-negative. Thus, if A 0 -L p C is also non-negative, then ϵ k , ϵ k ≥ 0 for all k ∈ Z + under Lemma 5.2. Now, let us address the conditions for the stability of system (5.34). Since this system is nonlinear, these conditions will be derived by denoting Z k = vec(z k , z k ), which allows us to rewrite (5.34) as: (5.36) where A + and A -are the same as in (5.30)

Z k+1 = A 0 + Lp C 2 Z k + A + Z + k + A -Z - k + ϱ k ,
, Lp = diag(L - p , L - p ) ∈ R 2n×2p , ϱ k = vec(ϱ k , ϱ k ) and C 2 =   C -C -C C   , ϱ k = B 0 u k + ∆Bu + k + w k , ϱ k = B 0 u k -∆Bu - k + w k .
The next theorem presents conditions that any gain L p has to fulfill in order to render IP (5.36) stable:

Theorem 5.5. Let assumptions 5.4-5.2 be satisfied and let L p be a given gain. If there exist

a matrix P1 ∈ R 2n×2n , diagonal matrices Q 1 , Q 2 , Q 3 , Ω + , Ω -, Ψ, Γ ∈ R 2n×2n
, such that the following linear matrix inequalities are verified:

           P1 -Q 1 -Ω + -Ω -0 ( P1 A 0 + P1 Lp C 2 ) ⊤ ⋆ -Q 2 -Ψ 0 ⋆ P1            ⪰ 0 P1 ≻ 0, Γ ≻ 0, Q 1 , Q 2 , Q 3 , Ω + , Ω -≥ 0, Q = Q 1 + min{Q 2 , Q 3 } + 2 min{Ω + , Ω -} > 0,
then system (5.36) is ISS with respect to the input ϱ ∈ ℓ 2n ∞ . Proof. The proof follows the same rationale as in Theorem 5.4 and is only sketched. By considering a Lyapunov function candidate given by V k = Z ⊤ k P1 Z k , one can show input-tostate stability, i.e.,

V k+1 -V k = -Z ⊤ k QZ k + ϱ ⊤ k Γϱ k
and the LMI conditions are stated by applying the Schur complement in the resulting stability conditions.

Finally, sufficient conditions for the existence of a gain L p providing non-negativity of the matrix A 0 -L p C and stability of (5.36) simultaneously are obtained by the following corollary:

Corollary 5.2. Let assumptions 5.1-5.4 be satisfied. If there exist diagonal matrices P2 ,

Q 1 , Q 2 , Q 3 , Ω + , Ω -, Ψ, Γ ∈ R 2n×2n
and U + , U -∈ R n×p such that the following linear matrix inequalities are verified:

P2 A 0 -Ũ + C 1 + Ũ -C 1 ≥ 0            P2 -Q 1 -Ω + -Ω -0 ( P2 A 0 + Ũ -C 2 ) ⊤ ⋆ -Q 2 -Ψ 0
Then, if the control u k in (5.36) is selected such as 5.38) where W k = vec(w k , w k ) and for some K, K -, K + , R ∈ R m×2n , the resulting dynamics is given by 5.39) where

u k = KZ k + K + Z + k + K -Z - k + RW k (
Z k+1 = KZ k + K + Z + k + K -Z - k + DW k (
D = I 2n + B 0 R, K = D z + B 0 K, K + = A + + B 0 K + and K -= A -+ B 0 K -, in which,
for ease of notation, denote B 0 = [B ⊤ 0 , B ⊤ 0 ] ⊤ . This brings us to the following result:

Theorem 5.6. Let assumptions 5.4-5.6 be satisfied. If there exist diagonal matrices

P , Q 1 , Q 2 , Q 3 , Γ, Ω + , Ω -, Ψ ∈ R 2n×2n and W 1 , W 2 , W 3 , W 4 ∈ R m×2n such that the following inequalities are verified            P -Q 1 -Ω + -Ω -0 W ⊤ 1 B ⊤ 0 + P D ⊤ z ⋆ -Q 2 -Ψ 0 W ⊤ 2 B ⊤ 0 + P A ⊤ + ⋆ ⋆ -Q 3 0 W ⊤ ⋆ P            ≻ 0 P > 0, Γ ≻ 0, Q 1 , Q 2 , Q 3 , Ω + , Ω -≥ 0, Q = Q 1 + min{Q 2 , Q 3 } + 2 min{Ω + , Ω -} > 0
then system (5.36) under control (5.38)

with gains K = W 1 P -1 , K + = W 2 P -1 , K -= W 3 P -1 , R = W 4 P -1
is input-to-stable stable with respect to the inputs W k .

Proof. Let us consider a Lyapunov function candidate given by V k = Z ⊤ k P -1 Z k , whose increments are given by

V k+1 -V k = =        Z k Z + k Z - k W k        ⊤ Π        K ⊤ P -1 K -P -1 + Q1 K ⊤ P -1 K + + Ω+ K ⊤ P -1 K -+ Ω- K ⊤ P -1 D ⋆ K ⊤ + P -1 K + + Q2 K ⊤ + P -1 K -+ Ψ K ⊤ + P -1 D ⋆ ⋆ K ⊤ -P -1 K -+ Q3 K ⊤ -P -1 D ⋆ ⋆ ⋆ D⊤ P -1 D - Γ              Z k Z + k Z - k W k        -Z ⊤ k Q1 Z k -Z + k ⊤ Q2 Z + k -Z - k ⊤ Q3 Z - k + W ⊤ k ΓW k -2Z ⊤ k Ω+ Z + k -2Z ⊤ k Ω-Z - k -2Z + k ΨZ - k ,
for any Ψ ∈ R 2n×2n , non-negative definite Qi , Ω+ , Ω-∈ R 2n×2n , i = 1, 3, and positive definite Γ. Then, if Π ⪯ 0, we have that As a direct implication of property (5.40), the following invariant set can be estimated:

V k+1 -V k = -αV k + W ⊤ k ΓW k (5.40) meaning that V k is an ISS Lyapunov function, provided that Q = Q1 + min{ Q2 , Q3 } + 2 min{ Ω+ , Ω-} ⪰ αP -1 (such an α always exists if Q ≻ 0).
X = x ∈ R 2n : x ⊤ P -1 x ≤ α -1 sup k≥0 W ⊤ k ΓW k
for α > 0 given in the proof of Theorem 5.6, and is an invariant set for (5.39).

In accordance with item 4 in 5.2, the following assumption (which is conventional in MPC) is imposed:

Assumption 5.7. Let X f ×X f ⊆ X ⊆ X×X and in (5.38) u k ∈ U provided that Z k ∈ X f ×X f .
However, if the set U is ellipsoidal, one may relax Assumption 5.7 by adding additional LMIs:

Corollary 5.3. Let there exist symmetric and positive definite matrices U ∈ R m×m and

Z ∈ R 2n×2n such that U = {u ∈ R m : u ⊤ Uu ≤ 1} and W k ∈ {W ∈ R 2n : W ⊤ ZW ≤ 1}
, and the conditions of Theorem 4 be satisfied with additional inequalities:

η ακ Γ ≤ min{κ -1 Z, P }, P ≥ κZ -1 ,         η 3 P 0 0 W ⊤ 1 + W ⊤ 2 0 η 3 P 0 W ⊤ 3 -W ⊤ 1 0 0 κ 3 P W ⊤ 4 W 1 + W 2 W 3 -W 1 W 4 U -1         ≥ 0 (5.41)
for some constants η > 0 and κ > 0, then control (5.38)

satisfies the constraint u k ∈ U for all Z k ∈ X f × X f .
Proof. This proof follows the same lines as the one of Corollary 5.1. However, in this case, the control law is designed for Z k instead of the interval center. To highlight the different steps, note that Z k = Z + k -Z - k . Therefore, the control (5.38) can be rewritten as

u k = (K + K + )Z + k + (K --K)Z - k + RW k ,
and thus the condition u k ∈ U takes the form:

     Z + k Z - k W k      ⊤      K + K + K --K R      U      K + K + K --K R      ⊤      Z + k Z - k W k      ≤ 1. If Z k ∈ X, then obviously Z + k , Z - k ∈ X.
The remainder of the proof remains identical.

Numerical illustration

To illustrate the developments of this section, consider the following LPV system: As it can be seen, the stabilization of (5.42) is successful and the computed bounds, for both IO and IP, respect relations (5.2) and ( 5.3). 5.42), in the case of higher measurement noises. In this simulated scenario, the state variable x 2 (left plot) is considered as available for measurement. As it can be seen in several time instants, due to this bigger measurement noise, the IP offers less conservative bounding estimates.

x k+1 =   0.5 0.6 + θ k θ k 0.3   x k +   0 1   u k + w k y k = 0 1 x k + v k (5.42) where x k = vec(x 1 , x 2 ) ∈ R 2 , θ k ∈ Θ = [-0.2,
Remark 5.5. Note that for the LTI case, constraint (5.44d) can be simplified and, instead of requiring the whole predicted interval z k,i+N to belong to X f × X f , only that the interval center must do so, i.e., z ⋆ k,i+N ∈ X f .

The cost function V N is defined as

V N = V f (x k,N ) + N -1 i=0 x ⊤ k,i+1 Ψ 2 x k,i+1 + s ⊤ i Ψ 3 s i ,
where the terminal cost V f and the terminal set X f are the stabilizing ingredients discussed, for each case, in the previous sections. The positive-definite matrices Ψ 2 , Ψ 3 represent the stage costs on the state and control input, respectively. Finally, the proposed IO-MPC algorithm is described in Algorithm 1:

Algorithm 1 IO-MPC Offline: Solve the appropriate LMIs to design the pair IO/IP and the respective controller. Estimate the terminal set X f and select Ψ 1 = P -1 , Ψ 2 ≤ α 2 P -1 , and Ψ 3 ≤ α 8 P -1 . Input: Initial conditions x 0 , x 0 , prediction horizon N . Online:

1: for each decision instant k ∈ Z + do 2: Measure y k and update the IO. 3: Initialize the IP such as z k,0 = vec(x k , xk ). 4: Solve OCP (5.43) under constraints (5.44a)-(5.44d). 5: Assign u k = s k 0 and apply it to system (5.1). 6: end for Theorem 5.7. Let [x 0 , x 0 ] ⊂ X and assumptions 5.1-5.7 be satisfied with [w k+1 , w k+1 ] ⊆ [w k , w k ] for all k ∈ Z + . Then, following Algorithm 1, the closed-loop system composed by (5.1), O(z k , y k , u k ) and P(z k , u k ) has the following features:

1. Recursive feasibility of reaching the terminal set in N steps; 2. Constraint satisfaction.

3. ISS of the dynamics of P in X f and practical ISS for system (5.1); Proof. Suppose that for any [x 0 , x 0 ] ⊂ X a solution of OCP (5.43) exists, i.e., there is a sequence of inputs S 0 N that leads the trajectories of P to X f . This means that for k = 0, by applying u k = s k 0 , we ensure that [z k+1,i , z k+1,i ] ⊂ z k,i+1 , z k,i+1 ⊂ X at least for i = 0 and also z k,N , z k,N ⊂ X f . Then, following Algorithm 1, the procedure can be iteratively repeated for k ∈ Z + since [w k+1 , w k+1 ] ⊆ [w k , w k ]. Moreover, the control sequence S N that steers z k,0 , z k,0 to X f also steers x k (as a consequence of (5.3)). This implies point (1).

For point (2), the ISS property in X f follows directly by the selection of the terminal ingredients (the choice of Ψ j , j = {1, 2, 3}, as given in Algorithm 1) which, due to the results of theorems 5.3 and 5.6 (for LTI and LPV systems, respectively), guarantees that

V (z k+1,N ) -V (z k,N ) ≤ d ⊤ k Γd k -ℓ(z k,N , s N )
where d k is the concerned disturbances depending on the case (either w ⋆ k or W k for the LTI or LPV case, respectively). Note that the matrix Γ are the same appearing in results (5.20) and (5.40). Furthermore, the practical ISS follows from the fact that

|x k | ≤ |z k |, |z 0 | ≤ |z 0 | + |z 0 | ≤ |x 0 | + c where c = |x 0 | + |x 0 |.
several relaxations for numerical tractability and the price is, obviously, an increased number of variables and conservativeness.

Numerical examples

In this section, numerical simulations are presented to illustrate the MPC algorithms proposed in this chapter. All simulations have been implemented in MATLAB 2017a, using an Intel i7-8565U processor (1.8GHz) and 16GB RAM. The toolboxes YALMIP [Löfberg, 2004], MPT3 [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF] and PnPMPC [START_REF] Riverso | Pnpmpc toolbox[END_REF] were used. The solvers employed in the optimization routines are discriminated in each case study.

The LTI case

To illustrate the proposed MPC algorithm, consider the CSTR model given as an example in Section 5.2.4. In this scenario, the constraints on state and control are given by X = [-2, 2] × [-10, 5] and U = [-4.5, 4.5], respectively. The disturbance sets are assumed to be given by W

= [-0.02, 0.02] × [-0.2, 0.2] and V = [-0.3, 0.3].
The prediction horizon is selected as N = 10, and weighting matrices H = 1000I 2 and R = 0.001. For the IO-MPC, the terminal set and the controller K f = [-6.99, -0.50] are obtained by solving the conditions given in Proposition 5.3, whereas the terminal cost is defined as its correspondent Lyapunov function.

Comparison with Tube-MPC

For comparison purposes, we have implemented the Tube-based MPC from [START_REF] Mayne | Robust output feedback model predictive control of constrained linear systems: timevarying case[END_REF]. For this implementation, a Luenberger observer was designed by pole placement (for the pair (A ⊤ , C ⊤ )). The terminal set and the static controller K LQ = [5.58, 0.45] are obtained by computing the associated LQR with weighting matrices Q LQ = 0.1I n and R LQ = 0.1. Figure 5.7 shows the comparison of the feasible sets for both techniques. As expected, due to the constraint tightening used to guarantee robust constraint satisfaction in the Tube-MPC, the obtained feasible region for the IO-MPC is much wider (and similar to the nominal MPC). Furthermore, Figure 5.7 also illustrates the impact of a naive design of the gains for the Tube-MPC. Clearly, the obtained regions are much different depending where the closedloop poles of the terminal controller/observer are placed, whereas our method has all gains readily obtained by conditions in the form of LMIs, which reduces the number of parameters to be tuned. Figure 5.7: Comparison of the feasible regions (approximations obtained by ray-shooting using YALMIP [Löfberg, 2004]). For the Tube-MPC, these regions were obtained by placing different poles when designing the Luenberger observer.

For the Tube-MPC, we selected the Luenberger observer which offered the biggest feasible region (by trial-and-error), which is initialized with x0 = (-0.65, -7) ⊤ , while the IO is initialized with x 0 = (-0.8, -7.8) ⊤ and x 0 = (-0.5, -6) ⊤ . The reactor is then simulated over a time window T simu = 60 time steps, considering several initial conditions satisfying x 0 ∈ [x 0 , x 0 ] (note that this setup is feasible for both techniques) and several realizations of w k and v k . For comparison purposes, let us define a performance index 5.1 illustrates the average results for 25 simulation runs. These results show that, even for a region where both techniques are feasible, the IO-MPC shows a faster solution of the optimization problem, as well as a better performance index.

J p = T simu k=0 x ⊤ k x k + u ⊤ k u k . Table
This latter fact was expected, since the control constraint was not tightened, meaning that the MPC algorithm can use its full range to achieve faster stabilization. However, since the Tube-MPC is designed to control the trajectories of the observer, it stabilizes the system closer to the origin, whereas the IO-MPC can only drive the envelope of trajectories close to it.

Simulation results

Let us initialize the IO with x 0 = (1.2, -8) ⊤ and x 0 = (1, -9) ⊤ . Notice that this region is not feasible at all for the Tube-MPC presented previously. For simplicity, let [z 0 , z 0 ] = [x 0 , x 0 ] and, again, let the reactor be simulated considering several x 0 ∈ [x 0 , x 0 ]. This simulation scenario is shown in Figure 5.8, for several realizations of the perturbations v k and w k . As it can be seen, all trajectories of the perturbed system satisfied the constraints and were stabilized close to the origin. Furthermore, it is worth noticing that the system trajectories are able to get very close to the constraint boundaries, indicating very low conservativeness. Finally, Figure 5.9 shows the input applied to the system.

The LPV case

To illustrate the proposed MPC algorithm, consider the model (5.42), given in Section

A 0 =   0.5 0.6 0 0.3   , A(θ) =   0 θ θ 0   .
and interpolating functions λ 1 = θ-θ θ-θ and λ 2 = θ-θ θ-θ . We considered initial conditions for the IP/IO as x 0 = vec(-7, -12) and x 0 = vec(-6, -10) and several initial conditions for (5. The control input computed by solving OCP (5.43) also fulfilled the constraints, as shown in Fig. 5.11. Finally, Fig. 5.12 shows the estimated feasible regions for the initial conditions of IP (5.34). We did not observe any improvement after N = 10. The mean computation time for solving OCP (5.43) was 0.22 ± 0.0313 second/step, with a maximum of 0.7725 second.

An interior-point method provided by the fmincon solver was used to solve the optimization problem.

Figure 5.10: Trajectories of the system and IP in contrast to the constraint set.

Part conclusion

This part of the thesis was devoted to the robust output feedback MPC problem, applied to both LTI and LPV cases. New interval estimators (both observers and predictors) were ). These regions were computed using YALMIP [Löfberg, 2004] and MPT [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF].

proposed, as well as their corresponding stabilizing controllers. All design conditions are formulated as offine LMIs.

These new estimators were incorporated into the classic MPC algorithm: the interval observer updates the set-membership of the states at every decision instant, and the interval predictor then casts an envelope containing all possible trajectories of the disturbed system into the prediction horizon. Using this feature, constraints on state are verified with a low computational burden.

Due to the form of these estimators, the stability analysis of the MPC algorithm is car-

Conclusions

In this thesis, the control and estimation problems were considered for two different models of microbial growth in continuous culture: one describing the competition of two species for a single limiting substrate, and another describing a mutualistic relationship between a producer and a cleaner strain.

These systems are structurally similar and highly complex: the dynamics of the substrates and the species are coupled by nonlinear functions of both variables, they are constrained both in states (since each state relates to a concentration, which characterizes such systems as positive) and control (the actuating variables take only non-negative values), and are possibly uncertain. These particularities lead to several challenges when analyzing and designing the respective state estimators and controllers, since many of the tools available are not applicable in these scenarios.

Concerning the estimation problem for such systems, the following topics have been addressed:

1. Observability analysis considering the availability of two realistic measurements (the total biomass and a fluorescent reporter). Since the first measurement is more common in real experiments, observability is also studied considering solely this measurement.

2. Design of state observers for both cases.

3. The advantage of having both measurements, in terms of the needed model parameters, is also discussed.

Indeed, if both measurements are available, it has been shown that the system is observable, whereas the case in which only the total biomass is available leads to merely detectability.

The proposed state observers are based on the use of differentiators (for the first scenario), and by using nonlinear estimators with output injection (for the second scenario).

Although the observers based on output injection show better performance (as investigated through numerical experiments), they require knowledge of all parameters of the model. If these parameters are not perfectly known, more errors will be introduced in the resulting estimation. On the other hand, the observers based on the differentiators are intrinsically sensitive to noise.

Concerning the control problem, the following topics have been tackled:

1. Robust control of the competition model (2.1), considering that all parameters are unknown but belonging to bounded interval;

2. Design of a control methodology for the COSY model (2.4), considering that all parameters are perfectly known.

For the first item above, the robust control architecture is composed of several discontinuous control laws that explicitly take into account the known bounds on the parameters. In fact, the dilution rate is actuated by a single control law, aiming to stabilize the concentration of a first species, whereas the one concerning the inflow substrate concentration is switched between to different control goals. This is due to the fact that, depending on the concentration of the substrate, the system might loose its controllability. Then, first this problem is avoided and, once a certain condition is reached, the goal then becomes the stabilization of the remain species.

For the microbial consortium described by the COSY model, however, the problem is more intricate. Indeed, the presence of acetate (excreted by the producer strain) further constrains the control. The concentration of the producer strain is regulated by state-feedback through the dilution rate, whereas the concentrations of glucose and acetate are then steered to specific levels, using the respective inflow concentration. These specific levels are selected in order to guarantee satisfaction of three constraints, which basically guarantee controllability and the permanence of the cleaner strain (i.e., in a concentration higher than zero).

Perspectives

A natural continuation of the theoretical work presented in this thesis is the design of a robust control algorithm for the COSY model. Indeed, the control laws developed in Chapter 4 are highly dependent on the parameters of model (2.4). This could be done, similarly as for the competition model, by using sliding-mode controllers that explicitly take the bounds of the parameters into account. Clearly, a more intricate stability analysis would be necessary.

Another interesting directions are (i) developing also a control law for B c , aiming to stabilize it at a certain level (this would reduce the consumption of glucose by this species and possibly enhance the process yield), (ii) developing an extra layer in the control architecture to encompass also some optimization criteria, (iii) couple the developed estimation algorithms to these control architectures, what would culminate in a methodology closer to the real-world.

Since the developments presented in this manuscript were theoretical, an interesting perspective is the implementation of the proposed algorithms into a real experimental set-up.

This would allow, for instance, a better understanding of the levels of uncertainty present (both parametric and noise-related), as well as the variability of parameters. Also, since real controller would have to undergo discretization schemes, practical experiments would also permit the evaluation of having actual sensors and non-ideal actuators.

Robust output feedback MPC Conclusions

The control of constrained system is a challenging problem in the control literature. This fact led MPC to be an interesting topic and several extensions were reported in the literature.

In this framework, the second part of this thesis dealt with the problem of robust output feedback MPC.

Indeed, by merging the classic MPC algorithm with interval estimators (an observer and a predictor), the algorithm proposed in this thesis guarantees both robust constraint satisfaction and recursive feasibility, at a low computation complexity and ease of implementation. Two broad classes of systems are encompassed by this algorithm: linear and linear parameter varying systems (and also to the case of linear time-delayed systems, whose results were not reported in this manuscript).

For each of the cases mentioned above, several design conditions for these interval estimators, as well as for their respective feedback control, were given in the form of linear matrix inequalities. This is an interesting feature, especially if compared to other solutions available (for instance, Tube-MPC or the ones using estimation techniques such as moving-horizon or zonotopes) which rely on complex set-algebraic computations, since LMIs are widely used and easy to implement and verify.

Although the design conditions given for the interval observer are not restrictive, a drawback of this methodology concerns the design of the interval predictor: due to its structure, it may be hard to find proper gains that guarantee its stability. This fact possibly shrinks the range of possible applications.

Perspectives

The versatility of the interval estimators used in the proposed IO-MPC algorithm offers some interesting research perspectives. For instance, extra estimation techniques could be used to identify parameters and exogenous inputs and enhance the accuracy of the interval Ttwo interesting phases of this methodology can be discriminated: the reaching phase, which is guaranteed if condition (B.3) is satisfied, it describes the motion of the system towards the surface σ(x), and the sliding phase, describing the motion of the system in the surface σ(x) = 0.

Recently, the concept of higher-order sliding-mode control has been widely studied. The techniques developed in this framework extend all interesting properties of standard slidingmode control to systems with a relative degree greater than one. To introduce this concept, consider the case of relative degree 2 and assume that the dynamics of the sliding variable σ satisfies the following non-linear system: order to steer σ and σ to zero in finite time, such as the twisting controller [START_REF] Emelyanov | High-order sliding modes in binary control systems[END_REF], the suboptimal controller [START_REF] Bartolini | Output tracking control of uncertain non linear second order systems[END_REF], quasi-continuous controllers [Levant, 2005] and others.

The properties of the suboptimal control algorithm are substantiated in the following.

This algorithm is given by the control law [START_REF] Bartolini | Output tracking control of uncertain non linear second order systems[END_REF] As one can see, control law (B.5) is a functional, since it requires information of its current and past trajectories of the system states. Also, it is interesting to notice that this control law does not require information on σ(t), but only the detection of an event where σ(t) = 0. This clearly offers an advantage for practical implementations, since information of σ(t) might not be always available and its estimation might suffer from numerical complications.

Denoting σ = y and invoking the work presented in [START_REF] Polyakov | Unified lyapunov function for a finite-time stability analysis of relay second order sliding mode control systems[END_REF] where t ′ = {inf t ≥ 0 : σ(t) = 0}.

Proof. This proof relies on the integration of the expression of the worst-case trajectories for the closed-loop system, i.e., σ(t) ≤ C -b min k 1 sign(∆(t)). Initially, as σ(0) > 0 and σ(0) = 0, we have that σ ≤ -r, for r = b min k 1 -C > 0.

Considering a general time interval [t 0 , t] (where t 0 is a certain initial time instant), the integration of σ(t) results in σ(t) ≤ σ(t 0 ) -r(t -t 0 ) σ(t) ≤ σ(t 0 ) + σ(t 0 )(t -t 0 ) -r (t -t 0 )

4. 1 •

 1 Random Parameter Sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1 Comparison between IO-MPC and Tube-MPC (time simulations) . . . . . . The sets of real and integer numbers are defined by R and Z, respectively, R + = {s ∈ R : s ≥ 0} and Z + = Z ∩ R + .

  |x| ∞ = ess sup t∈R ∥x(t)∥, and define L ∞ (R, R n ) as the set of all such functions with finite norms | • | ∞ ;

Figure 1 . 1 :

 11 Figure 1.1: Schematic of a computed-based control of biological process

2 .

 2 Develop a control algorithm for the co-population, aiming at a certain task. This control should take into account possible uncertainties and robustly achieve the control problem; 3. By developing the interval observers and predictors, propose new robust output feedback MPC algorithms that are less computationally expensive, less conservative and easier to implement. Therefore, this thesis is divided into four parts: the general introduction, the results concerning control and observation of microbial communities, the development of new robust output feedback MPC algorithms using interval estimators and, finally, general conclusions and outlooks. Part 2 is divided in two chapters: Chapter 2 addresses the problem of observation of microbial communities, considering a continuous bioreactor with a community of different microbes competing for a single limiting substrate and also the specific problem statement of the COSY model. Chapter 3 addresses the problem of robust stabilization of competing species at given concentration levels using a hybrid control architecture consisting of different discontinuous control laws. The control problem considering the COSY model is also presented. Part 3 is composed of a single chapter, which discusses the developments on the OF-MPC problem for both LTI and LPV systems. In this chapter, new interval estimators (observers and predictors), as well as their corresponding feedback control laws, are proposed. The MPC algorithms incorporate these estimators and stability is proven using classic arguments. It is shown that robust constraint satisfaction is assured and the complexity of the resulting algorithm is low. Numerical examples illustrate these novel algorithms.

Figure 1 . 2 :

 12 Figure 1.2: Overview of the thesis organization

Figure 2 Figure 2 . 1 :

 221 Figure 2.1: Illustration of the reactions rates (µ max = 1,K m = 1 and K i = 2).

Figure 2 . 2 :

 22 Figure 2.2: Inhibitory effect of A in the growth of the producer strain (parameters: µ max = 1, K g = 0.1, Θ a = 1).

Figure 2 .

 2 Figure 2.3 illustrates the pathways of the complete model.

Figure 2 . 3 :

 23 Figure 2.3: Pathways of the COSY model. The arrows indicate transfers, whose rates are depicted in red and purple. The transfers shown in red describe biological mechanisms (such as uptake, overflow, degradation, etc.), while the ones shown in purple are physical transfers (such as dilution rates). The symbol ∅ means that the refereed concentration is depleted out from the bioreactor (either by degradation or by dilution).

Figure 2 . 4 :

 24 Figure 2.4: Simulation of (2.4) for D = 0.2 [h -1 ] (dashed lines) and D = 0.3 [h -1 ] (solid lines). Left plot: trajectories of B p (in blue) and B c (in red). Center plot: trajectories of G. Right plot: trajectories of A.

FFigure 2 . 5 :

 25 Figure 2.5: Trajectories of the concentrations of the product H (left) and its fluorescent reporter F (right), in two scenarios: D = 0.3[h -1 ] (solid lines), and D = 0.2[h -1 ] (dashed lines).

Figure 2 . 6 :

 26 Figure 2.6: Numerical analysis of coexistence. Legend: green -only producer, red -coexistence, black -washout. This was obtained by numerical integration over a time span T = 6000h, having initial conditions given by G(0) = 0.5, A(0) = 0.1, B p (0) = 0.1, B c (0) = 0.1.

  result, consider the estimation errors given by e p = B p -Bp and e c = B c -Bc . Theorem 3.2. Consider system (2.4), differentiators (3.3)-(3.4) with gains satisfying (3.5), and observers (3.8)-(3.9). Then, the dynamics of the estimation errors e p and e c are inputto-state stable with respect to the measurement noises. Proof. Computing the time derivatives of e p , e c , one gets ėp = -(k + D)e p + (µ p -μp ), ėp = -(k + D)e c + (µ c -μc ). Consider a Lyapunov function candidate given by V = 1 2 (e 2 p + e 2 c ), whose time derivative is given by V = e p ėp + e c ėc = -(k + D)(e 2 p + e 2 c ) + e p (µ p -μp ) + e c (µ c -μc ).

  .10) Finally, recalling the results concerning the differentiator (see the Appendix B, Theorem B.1), the positive terms on the right-hand side of (3.10) are bounded. Then, provided that k + D > ϵ 2 (which is always verified by a proper selection of ϵ), one can conclude that the function V is an ISS Lyapunov function, proving the claim.

Figure 3 . 1 :

 31 Figure 3.1: Measurements y 1 and y 2 (in red), in contrast to the output of the differentiators (3.3)-(3.4) (in black).

Figures 3 .

 3 Figures 3.3 and 3.4 illustrate the real trajectories of B p , B c and G, A, respectively, as well as their respective estimates. Finally, Figure 3.2 shows the estimation of the heterologous protein H. Although the estimate Ĥ shows more sensitivity to noise than the other estimates,

FiguresFigure 3 . 2 :Figure 3

 323 Figures 3.5 and 3.6 illustrate, respectively, the estimation of the nutrients G and A, and biomasses B p and B c . As it can be seen, since all parameters are perfectly known, all estimates converge to the true values. Also, since they undergo integration, the influence of

Figure 3 . 4 :

 34 Figure 3.4: Estimation of G(t) and A(t).

Figure 3 . 5 :Figure 3 Figure 3

 3533 Figure 3.5: Estimation of G(t) and A(t) using observers (3.15).

Figure 4 . 1 :

 41 Figure 4.1: Mixing unit to modulate the concentration S in (t)

Assumption 4 . 1 .Assumption 4 . 2 .

 4142 For any a i ∈ A i and b i ∈ B i , i = 1, 2, let there exist a constant S a,b such that ξ(S a,b ) = 0 and ξ(S) > 0 for S > S a,b . For any a i ∈ A i and b i ∈ B i , i = 1, 2, the equality ∂ξ(So) ∂S = 0 is satisfied for some S o < S a,b .

Remark 4 . 1 .

 41 An option for selection is S m = So+S a,b 2 .

Proposition 4 . 2 .

 42 Consider model (2.1) with an output σ(S) = S -S m . By means of a switching law given by

Theorem 4 . 1 .

 41 Let S(0) ≤ S m . Then controls(4.5), (4.7) with χ ∈ (0, µ 1 (S(0))) and ϵ > 0 ensure that x 1 (t) = x d 1 and S(t) = S m for all t ≥ T , where T > 0 is a finite time dependent on the initial conditions.Proof. First, note that such a selection of χ induces D(t) > 0. Furthermore, as proven in Proposition 4.2 by means of Lyapunov function (4.6), the origin of (4.4) is globally finite-time stable after a time T 1 > 0 bounded by

Figure 4 . 2 :

 42 Figure 4.2: Convergence of the error (logarithmic scale). Legend: magenta -parameter set 1, green -parameter set 2, red -parameter set 3, blue -parameter set 4. The second-order nature of the control used for the variable x 2 is clearly seen in the behaviour of the convergence error.

Figure 4 . 4 :

 44 Figure 4.3: Stabilization at x 1 = 1 and x 2 = 1.2

Figure 4 . 5 :

 45 Figure 4.5: Stabilization considering actuator dynamics.

Figure 4 . 6 :

 46 Figure 4.6: Evolution of the control inputs (output of the actuator). The small boxes zooms the signal on a 1 hour time-window.

Figure 4 .

 4 Figure 4.7 illustrates this constraint as a function of G and A, highlighting that it is only transgressed for very low values of A.

Figure 4 . 7 :

 47 Figure 4.7: Illustration of the constraint (1 -k h )µ p (G, A) ≥ k + D min as a function of G and A and for D min = 0.05.

Figure 4 . 4 .

 44 Figure 4.4.1 illustrates the surface obtained by plotting the function h(G, A). The line rendering h(G, A) = 0 (and, consequently, the appropriate selection of G and A) is also shown. Now, let us investigate how steer the concentrations of G and A through the state-space.

Figure 4 .

 4 Figure 4.8: Surface of h(G, A). The white line represents the values of G and A rendering h(G, A) = 0.

Figure 4 . 4 . 1 ,

 441 in order to have h(G, A) = 0 in the steady-state, one must select A ⋆ = 0.0175 [gL -1 ].

Figure 4 .

 4 Figure 4.10 shows the trajectories of the concentrations of G and A. As it can be seen, both desired levels were attained roughly after t = 20h.Figure 4.11 then illustrates the trajectories

Figure 4 .

 4 11 then illustrates the trajectories of B p and B c , where one sees that the desired steady-state for B p = 0.3 [gL -1 ] was attained, while the B c converged to a positive level (≈ B c = 1.73 [gL -1 ]), as expected. Finally, Figure 4.12 shows the computed control inputs, reaching their steady-states at D = 0.414 [h -1 ],G in = 6.22 and A in = 3.25 (both in [gL-1 ]).

Figure 4 .

 4 Figure 4.10: Control of the concentrations of G and A.

Figure 4 .Figure 4 .

 44 Figure 4.11: Control of the concentrations of B p at 0.5 [gL -1 ] and the consequent evolution of B c .

  where w, w ∈ ℓ n ∞ and v, v ∈ ℓ p ∞ are known signals. Assumption 5.2. Let C ≥ 0. Assumption 5.1 imposes that the three sources of uncertainties in (5.8), i.e., x 0 , w k and v k are enclosed in given intervals, which is a classic hypothesis on the design of interval observers (IOs), whereas Assumption 5.2 is technical and can always be achieved by a proper change of coordinates. The objective of this chapter is to address the following problem statement: Problem 5.1. (OF-MPC) Let [x 0 , x 0 ] ∈ X and assumptions 5.1 be satisfied. The objective is to design an output feedback controller stabilizing system (5.1) in a vicinity of the origin while robustly satisfying state and control constraints x k ∈ X, u k ∈ U, ∀k ∈ Z + for any admissible realization of θ k , w k and v k , where X ⊂ R n and U ⊂ R m are given bounded convex sets.

Figure 5 . 1 :

 51 Figure 5.1: Illustration of the MPC algorithm.

Figure 5 . 2 :Figure 5 . 3 :

 5253 Figure 5.2: Evolution of the states. Legend -continuous lines: IO, dashed lines: IP, black lines: real system, blue lines: upper estimates, red lines: lower estimates, green lines: interval center.

Lemma 5 . 4 .

 54 being the observer gain to be determined. The fulfillment of relation(5.2) follows the cooperativity of the estimation errors e k = x k -x k and e k = x k -x k , the respective conditions are given in the following lemma: Let assumptions 5.4-5.5 be satisfied. Then, provided that A 0 -L o C is nonnegative, the estimation errors are non-negative, i.e., e k , e k ≥ 0 for all k > 0.

  .28) Now, computing the increments of the estimation errors e k , e k and taking (5.27)-(5.28) into account, one gets e k+1 = D o e k + r 1,1 + r 1,2 , e k+1 = D o e k + r 2,1 + r 2,2

  Decomposing Π and applying the Schur complement (similarly as done for matrix Σ in the proof of Theorem 5.1), the LMIs presented in this theorem are obtained by multiplying the resulting inequality by diag{P, P, P, P, P } from the left and right, and by introducing slack variables Γ = P ΓP , Ω + = P Ω+ P , Ω -= P Ω-P , Ψ = P ΨP , Qi = P Q i P , for i = 1, 3, and by introducing new decision variables W 1 = KP , W 2 = K + P, W 3 = K -P and W 4 = RP .

  0] and the disturbances are assumed to beenclosed in W = [-0.1, 0.1] 2 , V = [-0.5, 0.5].Solving the conditions presented in this chapter, the following gains are obtained for the pair IO/IP: gains were obtained as K 0 = 0.0004 -0.3034 , K 1 = 0.0999 0.0008 , K 2 = -0.1011 0.0015 .For simulation purposes, the IO is initialized with x 0 = [-10, 9] and x 0 = [-9, 10] (for simplicity, the IP is initialized with the same values), and θ k = -|0.2 sin(k)|.Figures 5.4 and 5.5 illustrate, respectively, the evolution of the states of (5.42) and the pair IO/IP, and the evolution of the control input.

Figure 5 . 4 :Figure 5 . 5 :

 5455 Figure 5.4: Evolution of the states. Legend -continuous lines: IO, dashed lines: IP, black lines: real system, blue lines: upper estimates, red lines: lower estimates.

Figure 5 . 6 :

 56 Figure 5.6: Comparison of the trajectories of the IO and IP, for model(5.42), in the case of higher measurement noises. In this simulated scenario, the state variable x 2 (left plot) is considered as available for measurement. As it can be seen in several time instants, due to this bigger measurement noise, the IP offers less conservative bounding estimates.

Figure 5 . 8 :

 58 Figure 5.8: Evolution of the states under the IO-MPC.

5 .Figure 5 . 9 :

 559 Figure 5.9: Control input under the IO-MPC

  42) satisfying x 0 ∈ [x 0 , x 0 ]. We ran 100 simulations of this scenario, each with a time span of T = 20 steps, considering several realizations of θ k , w k , and v k . The gains for the IP and the IO, obtained by solving the offline LMIs of Section III, are L o = [0.489, 0.1945] and L p = [0.232, 0.122]. For the MPC algorithm, we solve the conditions of Theorem 5.6 to obtain α = 1can estimate X f and select Ψ 2 . Finally, we select N = 10 and Ψ 3 = 0.001.

Fig. 5 .

 5 Fig.5.10 illustrates the trajectories of system(5.42) and the IP(5.36), in contrast to the constraint set. It is worth noticing that all constraints were respected, including those on the unmeasured state. The IP even reaches the boundary of constraint, indicating low conservativeness.

Figure 5 .

 5 Figure 5.11: Evolution of the control inputs

Figure 5 .

 5 Figure 5.12: Comparison of the feasible regions for prediction horizon with different lengths (N ∈ {4, 6, 8, 10}). These regions were computed using YALMIP[Löfberg, 2004] and MPT[START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF].

  ried out following the well-known stabilizing ingredients. Numerical examples were used to illustrate both methodologies.

  σ = a(t, x(•)) + b(t, x(•))u(t, σ, σ) (B.4) where |a(t, x(•))| ≤ C and 0 < b min ≤ b(t, x(•)) ≤ b max , for all x ∈ R. Also, constants C ≥ 0 and 0 < b min ≤ b max are supposed to be known. Many controllers have been proposed in

u = -k 1

 1 sign [∆(t)] (B.5) where k 1 , λ > 0, and∆(t) = ∆(t, σ(•)) := σ(t) -λσ(t * ) (B.6)where t * is the last instant of time in which σ(t * ) = 0, i.e.,t * = t * (t, σ(•)) := sup τ ≤t: σ(τ

1 Corollary B. 1 .

 11 , a candidate Lyapunov function for such a controller is given byV sub (t, σ(•), y(•)) =p σ(t) -y 2 (t) 2(-Csign(∆(t)) -b max k 1 sign(∆(t)) -y(t) -Csign(∆(t)) -b max k 1 sign(∆(t)) (B.8) where 0 < λ < 1, b min k 1 > C and λ estimate of the reaching time can be determined by t reach ≤ t ′ + kp |σ(t ′ )| (B.10) where t ′ is the first moment in time such that σ(t) = 0. Consider the suboptimal control (B.5) and dynamics (B.4). Let σ(0) = 0 and σ(0) > 0, then the following estimates are satisfied for all t ≥ 0: | σ(t)| ≤ | σ(0)| |σ(t)| ≤ |σ(t ′ )|

  2. S in (t): describing the concentration of substrates diluted in the inflow of medium (in unit gL-1 ). In the COSY model(2.4), this variable is represented by G in and A in are, respectively, the concentration of glucose and acetate.It is plausible to have the inflow of different media, for instance, containing or not diluted nutrients (D s (t) and D 0 (t), respectively). By considering a continuous bioreactor of constant volume, the outflow is equal to total inflow, i.e., D(t) = D s (t) + D 0 (t).The concentration S in (t) (or, analogously, G in (t) and A in (t)) can be modulated by properly selecting the rates D 0 and D s , which indicates that it also might be actuated by electromechanical pumps, as illustrated in Figure4.1. Considering that D s relates to the transfer of medium with a high concentration S max in , the inflow concentration S in can be determined by

Table 5 . 1 :

 51 Comparison between IO-MPC and Tube-MPC (time simulations)

		IO-MPC	Tube-MPC
	Perf. index J p	232.650	253.560
	OCP sol. time	0.0265	0.0524
	Final values (around origin)	[-0.070, 0.074] × [-0.993, 0.993] [-0.047, 0.057] × [-1.792, 1.743]

https://project.inria.fr/iplcosy/

1.4. Thesis outline

≤ δ.

-P ⋆ ⋆ ⋆ Γ P ⋆ ⋆ ⋆

+ P ⋆ ⋆ -Q 3 0 A ⊤ -P ⋆ ⋆ ⋆ Γ P ⋆ ⋆ ⋆

( P1 A + ) ⊤ ⋆ ⋆ -Q 3 0 ( P1 A -) ⊤ ⋆ ⋆ ⋆ Γ P1 ⋆ ⋆ ⋆

( P2 A + ) ⊤ ⋆ ⋆ -Q 3 0 ( P2 A -) ⊤ ⋆ ⋆ ⋆ Γ P2 ⋆ ⋆ ⋆

B ⊤ 0 + P A ⊤ -⋆ ⋆ ⋆ Γ W ⊤

4 B ⊤ 0 + P ⋆ ⋆ ⋆

Acknowledgments

First and foremost, I wholeheartedly thank the Lord for guiding me every single day of this journey. Under Your light, I have never felt truly alone along this journey. I'd like to express my eternal gratitude to my parents, Sebastião and Marilene, who have always made everything possible so I could become what I am today. You are the reason for all my efforts and your love always fueled my will of becoming a good and better man. I extend this gratitude to my supervisors Dr. Denis Efimov, Dr. Jean-Luc Gouzé and

time-varying parameters would appear such as

where the two right-most new terms can be seen as a perturbation (especially because ȧ and ḃ are supposed to be small, given the slow nature of the system) and, taking the uncertainty into account, they can be added to a max . The design of the control law then proceeds as explained in this subsection.

Numerical example

In this section, a numerical experiment is presented to illustrate the proposed methodology.

Consider model (2.1) with the reaction rates given as the Monod equation ( 2.2) as follows:

, µ 2 (S) = 2S 6 + S and with initial conditions S(0) = 1 [gL -1 ] and x 1 (0) = x 2 (0) = 1.5 [gL -1 ]. The control inputs are constrained such as D ∈ [0, 2] (in h -1 ) and S in = [0, 30] (in g L ). Also, for simplicity,

The objective is to stabilize x 1 at x d 1 = 2 [gL -1 ] and x 2 at x d 2 = 1.5 [gL -1 ]. The uncertainty will be considered as ±20% on each parameter. However, for the control synthesis, it is assumed that the parameter uncertainties satisfy the following intervals: 40, 4.60] ,

A 2 = [1.70, 2.30] [17.00, 23.00] , B 2 = [4.80, 7.20] To assess the robustness of the controller, four parameter sets will be considered for simulation purposes (i.e., the controller does not have such an information). These sets are shown in Table 4.1 as follows. [4.00, 20.00] [2.00, 6.00] Parameter set 2 [3.80, 24.00] [2.40, 7.20] Parameter set 3 [3.60, 22.00] [2.20, 5.40] Parameter set 4 [4.60, 23.00] [1.80, 4.80] Analyzing the impact of these uncertainties in S a,b and S o , one can arbitrarily select Consider the control u 1 as given by (4.19). Recalling that φ(σ g )σ g > 0 by definition, the derivative VG becomes Vg = -φ(σ g )σ g ≤ 0.

provided that σ g ≤ 0, or φ(σ g ) ≤ DG + f 1 .

However, σ g > 0 may lead to u 1 = 0 and the time-derivative Vg = σ g (-f 1 -DG) < 0.

This implies that V g is a Lyapunov function and the claimed stability property for σ g .

Theorem 4.3. Consider system (2.4) with an output σ p as given by (4.17).

and the conditions of Proposition 4.4 hold true. Then, by selecting D as given by (4.18), the closed-loop system is locally asymptotically stable with respect to the output σ p . Furthermore, if B c (0) > 0, then B c > 0 for all time.

Proof. Consider a Lyapunov function candidate given by V p = 1 2 σ 2 p and assume that G and A have converged to the desired steady-states G ⋆ and A ⋆ , respectively. Computing the time derivative of V p , one gets

By selecting D as given by (4.18) and thanks to the selection of the pair (G ⋆ , A ⋆ ), this derivative becomes Vp = σ p φ(σ p )B p < 0. which characterizes V p as a Lyapunov function. Assume now that G and A have not converged (due to, for instance, some transient response) and, therefore, the condition (1-k h )µ p (G, A)-(k + D min ) ≥ 0 is not satisfied. This would enforce D = D min and lead the time-derivative of V p to be Vp < -D min σ p ≤ 0, indicating that σ p is, at least, bounded. Once G and A reaches their steady-state, the previous conclusion concerning Vp will be valid and asymptotic stability is assured.

An interesting fact emerges when applying control law (4.18): the concentration of the product H, in its steady-state, becomes a monotonic function of µ p (G, A)B p . Indeed, plugging control law (4.18) into the dynamics of H and setting it to zero, the steady-state is readily computed as:

Therefore, stabilizing H at a desired concentration level basically sums up to choosing proper values of B ⋆ , G ⋆ and A ⋆ . Clearly, the control laws proposed in this section are suitable for such an objective. Figure 4.9 illustrates the heatmap of H ⋆ as a function of G and A.

Numerical example

In this section, numerical experiments will be presented to illustrate the control methodology developed in this chapter. The COSY model (2.4) will be simulated considering the initial conditions for the biomasses as B p (0) = B c (0) = 0.1 [gL -1 ], for the nutrients as

] and for the heterologous protein as H(0) = 0 [gL -1 ]. Arbitrarily, the

OF-MPC for LPV systems

In this section, a solution for Problem 5.1 will be investigated for the case in which system (5.1) is described by linear parameter-varying model [START_REF] Souza | Robust output feedback mpc for lpv systems using interval observers[END_REF], i.e.,

It is assumed that θ k is not measured, but its set of admissible values Θ is known. Furthermore, the matrix functions A : Θ → R n×n and B : Θ → R n×m are locally bounded and known. The measurement matrix C ∈ R p×n is assumed to be known.

Assumption 5.4. There exist matrices

. . , ν for some ν ∈ Z + , such that the following relations are satisfied for all θ ∈ Θ:

Assumption 5.4 (which is technical and it is introduced to simplify the writing) states that system (5.24) admits a convex embedding in a polytope defined by ν known vertices ∆A i and ∆B i with known centers A 0 , B 0 . Note that, since functions A, B and the set Θ are known, then there exists matrices A, A, B, B such that

Design of interval estimators

Interval observer

In this section, the objective is to design an IO for (5.24) by exploiting the available measurement. To this end, let us evoke Assumption 5.4 and rewrite (5.24) as

Then, under assumptions 5.4-5.5 and Lemma 5.1, replacing the uncertain terms in (5.24

Finally, since P > 0 and diagonal, the constraint

0, which is also linear in P and Ũ , finalizing the proof.

Interval predictor

Let the predictor gain (which replaces L in (5.24)) be L p . By definition,

Then, under Assumption 5.2, let us denote z k , z k as, respectively, the upper and lower predictive bounds of x k , and evoke Lemma 5.1 and Assumption 5.2 to write

Hence, the relation above allows us to rewrite (5.26) by replacing the terms unavailable for prediction (i.e., Ly k -Lv k + w k = LCx k + w k ) by their respective bounds:

where D p = A 0 -L p C. Hence, under assumptions 5.4-5.2, the IP (5.34) is composed solely by known terms.

In order to show that system (5.34) is cooperative, let us consider the prediction errors Proof. The increments of the prediction errors ϵ k and ϵ k are given by

where r 1,1 and r 2,1 are as in (5.29) and

(5.37)

then system (5.36) with gains L - p = P -1 2 U -and L + p = P -1 2 U + , is an IP for system (5.24), i.e., the relation z k ≤ x k ≤ z k is satisfied for all k ∈ Z + . Furthermore, the system (5.36) is ISS with respect to the input ϱ ∈ ℓ 2n ∞ . Proof. This proof follows directly from the implications of Lemma 5.5 and Theorem 5.5. Indeed, evoking [Efimov et al., 2013], the first constraint makes A 0 -L p C non-negative, thus satisfying the conditions on cooperativity imposed on Lemma 5.5.

Then, by introducing a new variable U = P 2 L p and recalling that U = U + -U -by definition, this search for L p is additionally constrained by the conditions stated in Theorem 2. This guarantees stability of (5.34), concluding the proof.

Control design

In this section, the control design for the IP (5.34) is addressed. To this end, the following simplifying assumption is needed: Assumption 5.6. Let ∆B = 0. Remark 5.4. Assumption 5.6 is imposed to streamline the presentation. Indeed, if system (5.24) is polytopic with ∆B ̸ = 0, the design conditions given in the following are affine and have to be checked over all of its vertices. This scenario requires an intricated presentation, which is to be avoided here.

Design of the predictive controllers

In the previous sections, the interval estimators that will be used in the MPC design have been presented. Also, conditions for the existence of a feedback control ensuring ISS, for each one of the cases considered, was derived. To synthesize the results for the both cases presented, the interval observer and predictor routines will be abbreviated as O(z k , y k , u k ) and P(z k , u k ), where z k = vec(x k , x k ) or z k = vec(z k , z k ), respectively.

Initialization

Consider the initialization of the MPC algorithm (analogously to constraint (5.5a)). As it was discussed in the previous sections, the IO uses the available measurement to update the set-membership of the states at every decision instant k, and thus can be used to initialize the predictor. However, it is worth noticing that the IO depends on the bounds of the measurement noise: if these bounds are too large, the obtained envelope of estimates can be more conservative. In this sense, the first step of the prediction can offer a more accurate estimation, as shown in Figure 5.6. Therefore, the MPC algorithm will be initialized using the following combined estimates:

The optimal control problem

Analogously to (5.4), the optimal control problem (OCP) for the IO-MPC can be formulated as follows:

where z k,i ∈ R 2n are the predicted state variables (containing both upper and lower bounds).

This OCP must be solved under the following constraints:

where P is the predictor used: IP (5.14) for the LTI case, and IP (5.36) for the LPV case.

For point (3), we have that the solution of OCP (5.43) implies that x k ∈ z k,k+1 , z k,k+1 ⊂ X due to relation (5.3) and u k = s k 0 ∈ U, under Assumption 5.3 (or Assumption 5.7 for the LPV case) and the discussed features of each possible set X f .

Complexity and performance

The constraints imposed on the OCP (5.43) is similar for both LTI and LPV cases, and scales with O(N n). Indeed, assuming that the number of hyperplanes needed to define sets X, U and X f depends linearly on n and m = n, the (worst-case) number of variables describing constraints (5.44a)-( 5.44d) is 8N n for the LTI case, and 10N n for the LPV case.

It is worth noticing that, for both LTI and LPV cases, the aforementioned complexity is fixed, which is an interesting feature of the proposed methodology.

The LTI case

For the linear case, the OCP (5.43) is a quadratic programming (QP) problem, being very similar to the nominal MPC. If compared to solutions based on Tubes [START_REF] Mayne | Robust output feedback model predictive control of constrained linear systems[END_REF], Mayne et al., 2009] or MHE [START_REF] Chisci | Feasibility in predictive control of constrained linear systems: the output feedback case[END_REF], the proposed method offers many advantages: (i) it does not require any steady-state assumption for the observer nor any further development for compensating the initial uncertainty, since this is automatically handled by the convergence of the pair IO/IP, (ii) the scheme is constructive, since all gains are obtained by the solution of LMIs (it should be noted that a naive selection of observer/controller gains in approaches that uses set approximations for propagating the uncertainties might dramaticallly influence the performance of the MPC, see the discussion in [Chisci and Zappa, 2002, Sec. 5.2] and the comparison below), (iii) thanks to the guar-

is similar to the one of the nominal MPC.

The LPV case

Approaches using zonotopic estimation often require extra procedures to limit their increasing complexity, as well as real-time knowledge on the scheduling parameter [START_REF] Wang | Zonotopic setmembership state estimation for discrete-time descriptor lpv systems[END_REF].

Techniques such as presented in [START_REF] Ding | Dynamic output feedback robust MPC with one free control move for LPV model with bounded disturbance[END_REF]Pan, 2016], that also assume no measurement of θ k , are computationally more expensive. Indeed, the authors fix a prediction horizon N = 1 and impose a min-max optimization problem accounting for all vertexes of the polytopic system, aiming to obtain a robust prediction. This implementation requires estimation in real-time. Clearly, this would imply a reduction in conservatism and better performance.

Also, slightly out of the scope of the MPC, it would be interesting to investigate if the robust constrained control problem could be solved using solely the interval observers. Obviously, an IO would allow the computation of the interval estimates merely one step ahead, instead of predicting it over an entire prediction horizon of arbitrary length. Although this fact possibly reduces the feasible region, it would alleviate the drawback of the design of the interval predictor mentioned above, enlarging the possible applications.

Appendix A

Preliminaries on input-to-state stability

Continuous-time case

Consider a general, continuous-time system given by

where

continuously differentiable function that satisfies f (0, 0) = 0. For each initial state x(0) = ξ and each u ∈ L ∞ , the trajectory of (A.1) is denoted x(t, ξ, u).

Definition A.1. [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF] System (A.1) is input-to-state stable if there exist a KL function β and a K function γ such that, for each input u ∈ L ∞ and each ξ ∈ R n , it holds that

for each t ≥ 0.

Definition A.2. [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF] A smooth function V : R n → R + is called an ISS-Lyapunov function for system (A.1) is there exists K ∞ -functions α 1 , α 2 and K-functions α 3 and α 4 such that

for any ξ ∈ R n and

for any ξ ∈ R n and any u ∈ R m so that |ξ| ≥ α 4 (|u|).

The relation between the above definitions characterizes a Lyapunov-like criterion for ISS [START_REF] Sontag | On characterizations of the input-to-state stability property[END_REF].

Lemma A.1. If system (A.1) admits an ISS-Lyapunov function, then it is ISS.

Discrete-time case

Consider a general, discrete-time system given by

where x : Z + → R n , u : Z + → R m and f (0, 0) = 0 is assumed. Let x(•, ξ, u) denote the trajectory of system (A.2) for a initial condition x 0 = ξ and the input u. Then, the concepts of input-to-state stability (ISS) and ISS-Lyapunov function are stated as follows Definition A.3. [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF]] System (A. 2) is input-to-state stable if there exists a KL-function β and a K function γ such that, for all ξ ∈ R n and u ∈ ℓ ∞ , it holds that

for each k ∈ Z + .

Definition A.4. [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF] A continuous function V : R n → R + is called an ISS-Lyapunov function for (A.2) if the following holds:

1. There exist K ∞ functions α 1 , α 2 such that

2. There exists a K ∞ function α 3 and a K function σ such that

The following result from [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF] establishes the relation between the definitions above:

Lemma A.2. If system (A.2) admits a continuous ISS-Lyapunov function, then it is ISS.

Appendix B Preliminaries on sliding-mode control

Sliding-mode and sub-optimal controller

Consider the following non-linear system:

where x ∈ R n is the state vector, u ∈ R m is the control input vector and d ∈ R n is a vector containing perturbations and possibly uncertainties within the system. First, the first-order conventional sliding mode is discussed. Consider the following surface in the state space:

where σ : R n → R is a continuously differentiable function called sliding variable.

The objective is to design σ(x) in such a way that system (B.1), under feedback control, behaves with a prescribed performance. In this sense, if the control input is selected guaranteeing σ = -ηsign(σ), where η > 0, one can consider a Lyapunov function as

and, computing its time derivative, we have that

which implies V ≤ 0, hence the origin of σ(x) is globally stable. Moreover, it is clear that |σ| = √ 2V 0.5 and, consequently, it implies that

meaning that the solution of (B.3) becomes zero in a finite time [START_REF] Shtessel | Sliding Mode Control and Observation[END_REF].

Then, it is clear that the first sign change in the control law will take place in a second instant of time t" > t ′ , in which we have (due to the form of (B.5)) that σ(t") = λσ(t ′ ).

Consequently, the evaluation of (B.11) in the time interval [t ′ , t"] shows that

and hence, we can finally conclude that

Hence, as contraction and convergence properties of such a control algorithm have been proven (see [START_REF] Bartolini | Output tracking control of uncertain non linear second order systems[END_REF]] and proofs therein), we can then conclude that | σ(t)| ≤ | σ(0)| for all t > 0.

The proof of the remaining estimate follows by analysing the signal change of the control law. Clearly, due to the definition of the functional ∆(σ(t)), the signal of the control law will switch every time that σ(t) reaches zero or when σ(t) = λσ(t * ) (see the definition of t * above). This fact allows us to deduce that, as σ(t") < 0, it is clear that σ(t) < 0 holds for all t > t ′ . Thus |σ(t)| < |σ(t ′ )| holds for all t > 0, as claimed.

Sliding-mode differentiator

As it will be discussed, the observers proposed in this thesis will require the time derivatives of the available measurements. This differentiation can be performed numerically by means, for instance, of the robust exact differentiator, proposed by Levant [Levant, 2003]. The n-th order differentiator is realized by the following observer: where f (t) = f 0 (t) + v(t) is a signal composed by an unknown base signal f 0 (t) with its n-th derivative having a known Lipschitz constant L > 0, and v ∈ L ∞ is the measurement noise.

The constants λ i are tuning parameters. Under (B.12), the estimation errors are bounded according to the following theorem:

Theorem B.1. [Levant, 2003] Let the input noise satisfy |v(t)| ≤ ϵ for almost all t ≥ 0.

Then the following inequalities are established in finite-time T > 0, for some positive constant 1, . . . , n. (B.13) According to Levant, all solutions of (B.12) are stable. Furthermore, the convergence of z 0 → f (t), z 1 → d dt f (t), . . . , z n → d n dt n f (t) is established in a finite-time in the noise-free case.