DESIGN OF A NONLINEAR OBSERVER USING THE FINITE ELEMENT METHOD WITH APPLICATION TO A BIOLOGICAL SYSTEM

An observer for a nonlinear biological system — biomass production in a bioreactor —is proposed. The speciﬁc growth rate is estimated. The key point of the observer design is ﬁnding a solution of a certain partial differential equation. Conditions guaranteeing existence of its solution are presented. The solution is approximated using ﬁnite element method. The results are illustrated by a numerical example.


Introduction
Observer problem for biological systems form a large part of practical applications of observer theory. One reason for this is enormous importance of biotechnology in nowadays, the other one being a consequence of the fact that many systems in biotechnology contain quantities that cannot me directly measured or, in some cases, measurement of these quantities is difficult or prohibitively expensive.
Let us mention few examples of observers used for control of processes in biotechnology. Due to a sheer number of results we can pick a small number of them. The first possibility is to use the so-called high-gain observers. These are in fact linear observers with gain designed sufficiently large to deal with nonlinear systems. High-gain adaptive observer for a biological systemwaste water treatment -was presented inČelikovský et al. (2018) with estimation of the specific growth rate.
Unfortunately, the common drawback of high gain observers for nonlinear systems is a strong sensitivity to noise. A discontinuous (sliding-mode) observer was derived in Daaou and Dochain (2017); De Battista et al. (2010); Vargas et al. (2014) etc. Finally, let us note Rehák and Papáček (2013) where the nonlinearity was treated as robustness which, however, allows only a rough approximation of the nonlinearity. This observer has a particular feature that it uses a combination of continuous and discrete measurements. These examples cover applications of three main directions of the observer theory: robust observers, sliding-mode observers and high-gain observers. For a thorough overview of observers for chemical and biological processes, see e.g. Mohd Ali et al. (2015).
A different approach to state observation of nonlinear systems was proposed by Kazantsis and Kravaris in Kazantzis and Kravaris (1998), further modified to systems with delayed measurement in Kazantzis and Wright (2005). Here, an equation -a counterpart of the Sylvester equation known from the observer problem of linear systems -was derived. This equation is a partial differential equation (PDE) of first order. Originally, expansion of all expressions into Taylor polynomials and seeking an approximation of the aforementioned PDE also in the form of Taylor polynomials was proposed as a viable method for solution. This method has rather restrictive assumptions: the approximation is guaranteed to exist if the linearization of the observed system has all eigenvalues in the left complex half-plane or all eigenvalues in the right half-plane.
This assumption was removed in Sakamoto et al. (2014) by proposing alternative method for approxima-tion of the PDE. This method is based on successive computation of trajectories of certain system of ordinary differential equations followed by interpolation of the results. Observe that this method is based on analogous method for the solution of the regulator equation arising from the nonlinear output regulation problem as described in Sakamoto and Rehák (2011). This observer construction was successfully applied to practical control problems, e.g. in Horibe and Sakamoto (2018).
The regulator equation known from the nonlinear output regulation problem is also a PDE of first order. Solution of this equation using finite-element method (FEM) was studied in Rehák andČelikovský (2008), Rehák et al. (2009) with proof of existence of FEM approximation of the solution given in Rehák (2011). As this equation is closely related to the PDE derived in Kazantzis and Kravaris (1998), the idea of application of FEM to the observer problem is straightforward. First attempt was presented in Rehák, B. (2018), in this case for systems with delayed measurements.
Let us mention that a related problem -estimation of parameters of a system -was successfully solved via the approach using similar PDEs for the case of a biological system (algae growth) in Papáček et al. (2010); . Hence also motivation for our research.
Motivated by the previous positive experience with applications of the observer design based on Kazantzis and Kravaris (1998) in connection with a successful application of FEM to the solution of PDEs arising from this or related problems, we propose to adopt this approach to the observer design to a biological process -to be specific, to the problem of estimation of the specific growth rate in a bioreactor.

Biomass Growth Equations
The biomass growth in a bioreactor can be described by the following differential equations (see e.g. De Battista et al. (2010) ( Equation (1) describes changes in the biomass concentration x dependent on the dilution rate D. The second equation, Eq.
(2) shows how the specific growth rate changes. The function ρ attains a specific form according to the kinetic model used. In most cases, the kinetics is reformulated using the Monod or Haldane curves: µ is a function of a variable s given by the equatioṅ where the specific growth rate is expressed as in case of the Haldane curve used in this paper or as (called the Monod curve). The parameters µ m , k s and k i are assumed to be known. The biomass concentration is supposed to be continuously measured, however the variables µ or s (after using the Monod or Haldane curves) are not measurable. Hence the need for an observer. Also, the dilution rate it is continuous and bounded, moreover, this function is known. Naturally, the biomass concentration is supposed to be strictly positive.
3 Nonlinear Observer Consider a nonlinear system (the plant) where F : (7) describes the output of the plant, the other quantities are supposed not to be measurable. For system (6,7) we construct an observer adopting the approach proposed in Kazantzis and Kravaris (1998). First, define matrixÃ ∈ R n×n , a vector b ∈ R n chosen so that the pair (Ã, b) is observable. Then, one constructs mapping Φ : R n → R n such that it obeys the following equation Then, the observer for system (6,7) is defined aṡ with the observer gain L(x) defined as Proposition 3.1. (Kazantzis and Kravaris (1998)) If a solution of Eq. (8) exists, then the observation error e = x −x satisfies lim t→∞ e(t) = 0.
Remark 3.2. Kazantsis and Kravaris prove existence of a solution of equation (8) under rather restrictive assumptions: the observed system must be either exponentially stable around the origin or all eigenvalues of the Jacobi matrix ∂F ∂x (0) have positive real parts. This is since the proof of existence of a solution of Eq. (8) relies upon the so-called Lyapunov auxiliary theorem which was used for construction approximations of Φ by Taylor polynomials. This restrictive assumption was removed in Sakamoto et al. (2014) where the solution was found using an iterative method. FEM was applied to solution of a related equation in Rehák, B. (2018), under milder assumptions than those required in Kazantzis and Kravaris (1998).
To improve numerical properties of the proposed method, the linear case is treated separately. If there exist matrices A, C so that system (6,7) is linear, that iṡ Eq. (8) attains the form of the Sylvester equation: Here,Φ ∈ R n×n is an unknown matrix. Note that for its solution suffices to assume max Re eigÃ < min Re eigA.
Let us turn our attention to the original nonlinear problem. Define matrices A, C as Moreover, denote also the remaining higher-order terms in F and h by ϕ, κ: Lemma 3.3. Let A, C, ϕ κ be defined as above, let ma-trixÃ is chosen so that (14) holds and the pair (Ã, b) is observable. LetΦ solve Eq. (13) and assume there exists a smooth function φ : R n → R n satisfying . (17) ThenΦ(x) + ϕ(x) solves Eq. (8).

Solution of the Nonlinear Equation
Eq. (17) is a partial differential equation of first order and, moreover, it is of rather unusual form. It must be solved numerically as, in general, one cannot expect to be able to find an exact solution analytically.
First issue is that one needs to search the approximation of the solution of Eq. (17) on a bounded open set (a domain) Ω ⊂ R 2 rather than on the whole Euclidean space. It is required that this domain contains the origin: 0 ∈ Ω. Moreover, it is assumed the domain Ω has Lipschitz boundary: the boundary of the domain Ω (denoted by ∂Ω) can be locally described as a graph of a Lipschitz-continuous function (for every point x ∈ ∂Ω there exists a neighborhood of x so that an open subset of the boundary containing the point x can be described as a graph of a Lipschitz-continuous function).
Existence conditions of a solution of Eq. (17) on a domain Ω are formulated in Lemma 1.6 in Roos et al. (1996). For the reader's convenience, this lemma is repeated here as theorems guaranteeing existence of a solution of this type of equations are not easy to find in literature. As an alternative, the reader can find details in Rehák (2011)  Let Ω ⊂ R n be a bounded domain with Lipschitz boundary, 0 ∈ Ω and n(x) be the outward normal vector defined at the point x ∈ ∂Ω. Assume a > 0, β ∈ (C 1 (Ω)) n and g ∈ L 2 (Ω).
As one can see, Lemma 4.1 guarantees existence of a solution for scalar equations. How to apply this lemma to the two-dimensional equation (17) shows the subsequent result: Lemma 4.2. LetÃ = diag(−a 1 , −a 2 ). Assume the following holds for i = 1, 2: Then there exists a neighborhood of the origin U ⊂ R 2 so that condition (18) is satisfied in U .

Proof. First, observe that divβ(x) = TraceA+divϕ(x).
Denote ω i = a i − 1 2 TraceA. Moreover, note that all derivatives of function ϕ are continuous and ∂ϕ ∂x (0) = 0, there exists a neighborhood of the origin denoted by U so that for all x ∈ U holds f (x) < ω i . This implies validity of (18).
The main result can be summarized as follows Theorem 4.3. Let Ineq.
5 Remarks About the Numerical Approximation of the Partial Differential Equation The theorems in the previous section form a basis for a numerical approximation of the function φ. As a method of this approximation, FEM was chosen thanks to its good approximation properties as well as availability of software capable of handling such problems.
The boundary conditions cause some undesired error. Fortunately, as numerical experiments show, if the mesh covering the domain Ω is sufficiently fine, this error is concentrated on a small area around the boundary. Another issue crucial for the precision of the results is sufficiently fine discretization in the area around the origin.
Note that since matrix T is invertible, as follows from properties of the Sylvester equation (see Birkhoff and Lane (2017)), inversion (T + ∂φ ∂x (x)) −1 at the origin always exists as ∂φ ∂x (0) = 0. Nevertheless, invertibility of the expression (T + ∂φ ∂x (x)) −1 forx = 0 is not guaranteed. Even if this matrix is nonsingular but the condition number of this matrix is small, computation of the observer gain is difficult.
As a remedy, a modified term was used in place of the inversion of the matrix (T + ∂φ ∂x (x)). Let us define Then, matrix T is symmetric positive definite, hence its eigenvalues are positive. Denote by u(x) the eigenvector corresponding to the smallest eigenvalue of T (x).
Choose also a parameter δ > 0. Then, the new observer gain is given by Note that if δ = 0 then L = L. Let us comment how to choose the matrixÃ. If its eigenvalues are too close to eigenvalues of matrix A, the observer will not be fast enough. On the other hand, eigenvalues with too large absolute values of their real parts will cause other problems. As known, the observer will be acting too fast, resulting in high sensitivity to noise. Moreover, the matrix T might have a small norm, thus small changes in the term ∂φ ∂x lead to large changes of the expression (T + ∂φ ∂x (x)) −1 so that even the aforementioned method might not be sufficient. Moreover, this property results in the requirement of a precise approximation of the function φ, hence the mesh must be very fine. This in turn poses high demands on the computational resources.

Example
The example is as in the paper De Battista et al. (2010). The system (1,2) is rewritten using the Haldane curve and dilution rate in the forṁ The values of the parameters were chosen as in De Battista et al. (2010). For the reader's convenience, we repeat them here in the following table.
µ m = 0.22, k s = 0.14, The dilution rate is chosen so that system (22,23) has equilibria x e = 13.84, s e = 0.1996. The output equation is y = x, hence C = (1, 0). Linearization around the equilibrium point (x e , s e ) yields matrix A as Observer is designed using parameters Solution of the linear equation (13) can be evaluated as The function φ was approximated using finite elements on the domain of elliptical shape with center at the equilibrium (13.84, 0.1996) and semiaxes parallel to coordinate axes. The length of the semiaxis parallel to the x-axis is 1, semiaxis parallel to the y-axis is 0.15. Trajectories of the simulated system must then remain in this domain. The approximation of PDE (17) was computed using the software Comsol Multiphysics. Advantage of this software is providing not only the values of the solution φ(x, s) but also its first derivatives. As these are needed for the observer construction, this is a facilitation of the observer implementation. The resulting functions φ 1 and φ 2 are depicted in Figs. 1 and 2, respectively.
Results of simulations can be seen in the following figures. Fig. 3 illustrates the state s (solid line) and its estimate obtained by the aforementioned nonlinear observer (dashed line). The observer gain L was replaced by the gain L with a parameter δ = 0.0002. For a comparison, dotted line represents the estimate obtained by a linear observer with observer gain L =Φ −1 b. One can see that the linear observer has a much larger initial undershoot. This is also visible in Fig. 4 where the norm of the observation error of both observers is depicted. The solid line stands for the nonlinear observer while the dashed line represents the linear observer. We can see faster response of the nonlinear observer. The state s in this case is shown in Fig. 5, the norm of estimation error is in Fig. 6. We can see a smaller initial undershoot but the price for this is a slower convergence.

Conclusion and Outlook
An observer for a biological system was presented. The goal was to estimate the specific growth rate. The main problem was a solution of a partial differential equation, this was carried out using the finite element method. Results are illustrated by a numerical example.
For the future work, observers with delay will be investigated as well as observers depending on a parameter. Also, the results of the present research will be used for synchronization in complex networks ), multi-agent systems (Rehák and Lynnyk (2019)) and robotics (Anderle andČelikovský (2017)).