Olivier Bernard
email: olivier.bernard@inria.fr

Georges Bastin
email: bastin@auto.ucl.ac.

Identification of reaction networks for bioprocesses: determination of a partially unknown pseudo-stoichiometric matrix

Keywords: Modelling, Nonlinear systems, Bioreactors, Validation

In this paper we propose a methodology to determine the structure of the pseudo-stoichiometric coefficient matrix K in a mass balance based model and to identify its coefficients from a set of available data. The first stage consists in estimating the number of reactions that must be taken into account to represent the main mass transfer within the bioreactor. This provides the dimension of K. Then we propose a method to directly determine the structure of the matrix (i.e. mainly its zeros and the signs of its coefficients). These methods are illustrated with simulations of a process of lipase production from olive oil by Candida rugosa.

Introduction and motivation

Macroscopic modelling provides simple dynamical models which have proved of great interest in bioengineering for the design of on-line algorithms for bioreactor monitoring, control and optimisation [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF][START_REF] Bastin | Nonlinear and adaptive control in biotechnology: a tutorial[END_REF]. In such works, the dynamical behaviour of a stirred tank bioreactor is often described by the following general macroscopic massbalance model:

dξ(t) dt = K r(t) + v(t), (1)
In this model, the vector ξ = (ξ 1 , ξ 2 , . . . , ξ n ) T is made-up of the concentrations of the various species inside the liquid medium. The term v(t) represents the net balance between inflows, outflows and dilution effects. The term K r(t) represents the biological and biochemical conversions in the reactor (per unit of time) according to some underlying reaction network. The (n × p) matrix K is a constant pseudo-stoichiometric matrix. r(t) = (r 1 (t), r 2 (t), . . . , r p (t)) T is a vector of reaction rates (or conversion rates). It is supposed to depend on the state ξ and on external environmental factors such as temperature, light or pressure, etc

The pseudo-stoichiometric (PS) matrix K is associated to a macroscopic reaction network that lumps together the many intracellular metabolic reactions of the various involved microbial species. The reaction network summarises then the main mass transfer throughout the bioreactor by a few reactions involving mainly extracellular compounds and biomasses without describing into all details the intracellular behaviour. Each column of the matrix corresponds to a chemical or biological reaction of the underlying macroscopic reaction network. The coefficients k ij j = 1, . . . , p are associated with the j th reaction. A positive k ij means that the i th species ξ i is a product of the j th reaction, while a negative k ij means that ξ i is a substrate of the j th reaction. If k ij = 0 the species ξ i is not involved in the j th reaction.

In this paper, we are concerned with modelling situations where the on-line concentrations ξ i of the involved species are measured but the structure of the reaction network is a priori questionable and therefore the matrix K is partially unknown. The objective, as in [START_REF] Bogaerts | Sytematic generation of identifiable macroscopic reaction schemes[END_REF], is to provide guidelines to the user for the identification of the structure of a macroscopic reaction network and the determination of the PS matrix K from the available data. Note that the method can also be applied to simplify a known detailed intracellular metabolic network and provide a simpler reaction network that represents the main mass transfers throughout the system and directly connect initial substrates to final products. In such a case the concentrations ξ i would result from simulations of a model based on the detailed known reaction network and matrix K is to be found from these "data".

The usual approach dedicated to the determination of reaction networks relies on the linearisation of the dynamics around a reference solution [START_REF] Eiswirth | Mechanistic classification of chemical oscillators and the role of species[END_REF][START_REF] Chevalier | Toward a systematic determination of complex reaction mechanisms[END_REF] and identification of the local jacobian matrix. This approaches are then suitable for data close to steady state. Here, in the spirit of [START_REF] Chen | Structural identifiability of the yield coefficients in bioprocess models when the reaction rates are unknown[END_REF][START_REF] Bernard | Structural identification of nonlinear mathematical models for bioprocesses[END_REF], we use linear algebraic properties to exploit the structure of the bioprocesses (Equation ( 1)) and our arguments do not rely on any linearisation. As a consequence we are not limited to steady state data and we can exploit all the available measurements, even when associated to transient states.

The problem is illustrated with the following example.

Example: Let us consider the example of a competitive growth on two substrates [START_REF] Serra | A model for lipase production by Candida rugosa[END_REF] which could represent, for instance, the production of lipase from olive oil by Candida rugosa. Here the microorganism is supposed to grow on two substrates that are produced by the hydrolysis of a primary complex organic substrate.

The following 3-step reaction network has been assumed in the literature [START_REF] Chen | Structural identifiability of the yield coefficients in bioprocess models when the reaction rates are unknown[END_REF]:

-Hydrolysis:

k 11 S 1 + E -→ S 2 + k 31 S 3 + E -Growth on S 2 : k 22 S 2 + k 62 O -→ X + k 72 P -Growth on S 3 : k 33 S 3 + k 63 O -→ X + k 43 E + k 73 P
where S 1 is the primary substrate (olive oil, made of several compounds, mainly triglycerides), S 2 (glycerol) and S 3 (fatty acids) are the secondary substrates. E is the enzyme (lipase), X the biomass (Candida rugosa), O the dissolved oxygen and P the dissolved carbon dioxide.

The associated PS matrix is:

K =           -k 11 0 0 1 -k 22 0 k 31 0 -k 33 0 0 k 43 0 1 1 0 -k 62 -k 63 0 k 72 k 73           with k ij > 0.
We shall assume that this reaction network is unknown to the user and has to be discovered from data of the species concentrations. Here the data will be simulated by a model but of course in practice the data are obtained from experiments.

Generally, the choice of a reaction network and its associated PS matrix K results from modelling assumptions. Sometimes however, a complete description of the reaction network is a priori not available. This can be a consequence of a lack of phenomenological knowledge on some of the involved mechanisms, letting a part of the reaction network questionable. The problem can also arise when it is desired to reduce a complicated given reaction network to a much simpler model. This situation especially occurs for models describing wastewater treatment processes involving a bacterial consortium made of a broad range of bacterial species degrading a mixture of organic substrates. For example, more than 140 bacterial species have been found [START_REF] Delbès | Bacterial and archaeal 16s rdna and 16s rrna dynamics during an acetate crisis in an anaerobic digestor ecosystem[END_REF] in an anaerobic digestion wastewater treatment plant.

We first propose a method to determine the size of matrix K ı.e. the number of independent reactions that are distinguishable from the available data. Then we show how the structure of matrix K can be estimated, using the a priori available knowledge on the process. By structure we mean the sign and the location of the non-zero entries of matrix K. In addition, the method can also provide an estimate of the parameters k ij if the available knowledge is sufficient.

2 Determination of the number of reactions

Introduction

In this section, we intend to determine the minimum number of reactions which are needed in order to explain the observed behaviour of the process, without any prior knowledge on the underlying reaction network. We assume that the vectors ξ(t) of species concentrations and v(t) of inflow/outflow balances are measured during some time interval and exhibit significant variations with time. We assume also that the number of measured variables is larger than the number of reactions: n > p. The PS matrix K and the vector of reaction/conversion rates r(t) are unknown.

Theoretical determination of dim(Im(K))

The model equation (1) can be viewed as a linear dynamical system with state ξ and inputs r(t) and v(t) (although we know obviously that r and v may be state dependent). If we take the Laplace transform of this equation, we get:

sΞ(s) = KR(s) + V (s) (2)
where Ξ(s), R(s) and V (s) are the Laplace transforms of ξ(t), r(t) and v(t) respectively. A linear filter or smoother with transfer function G(s) can then be used in order to clean the data (noise reduction, decrease of autocorrelations etc ...):

U (s) = KW (s) with U (s) = G(s)[sΞ(s) -V (s)] and W (s) = G(s)R(s).
Or, in the time domain:

u(t) = Kw(t) (3)
with u(t) and w(t) the inverse Laplace transforms of U (s) and W (s) respectively. The vector u(t) can be computed directly from the data by appropriate filtering/smoothing techniques possibly involving delay operators.

For example, the moving average is a very simple filter that can be applied to [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF], and provides an expression of the form (3) with (T denotes the considered moving average window):

u(t) = 1 T ξ(t) -ξ(t -T ) - t t-T v(τ )dτ (4)
and

w(t) = 1 T t t-T r(τ )dτ
Now the question of the dimension of matrix K can be formulated as follows: what is the dimension of the image of K ? In other words, what is the dimension of the space where u(t) lives ? Note that we assume K to be a full rank matrix. Otherwise, it would mean that the same dynamical behaviour could be obtained with a matrix K of lower dimension, by defining other appropriate reaction rates. The determination of the dimension of the u(t) space is a classical problem in statistical analysis. It corresponds to the principal component analysis (see e.g. [START_REF] Johnson | Applied multivariate statistical analysis[END_REF]) that determines the dimension of the vector space spanned by the vectors k i which are the rows of K. To reach this objective, we consider the n × N matrix U obtained from a set of N estimates of u(t):

U = (u(t 1 ), . . . , u(t N ))
We will also consider the associated matrix of reaction rates, which is unknown:

W = (w(t 1 ), . . . , w(t N ))
We assume that matrix W is full rank. It means that the reactions are independent (none of the reaction rates can be written as a linear combination of the others). We consider more time instants t i than state variables: N > n.

Property 1 For a matrix K of rank p, if W has full rank, then the n × n matrix M = U U T = KW W T K T has rank p. Since it is a symmetric matrix, it can be written:

M = P T ΣP
where P is an orthogonal matrix (P T P = I) and

Σ =             σ 1 0 . . . 0 0 σ 2 0 0 . . . . . . σ p 0 . . . . . . 0 . . . 0             with σ i-1 ≥ σ i > 0 for i ∈ {2, ..., p}.
Moreover, the eigenvectors associated with the σ i generate an orthonormal basis of ImK.

This property is a direct application of the singular decomposition theorem [START_REF] Horn | Matrix analysis[END_REF] since rank (M ) = rank (KW ) = rank (K) = rank (Σ) = p. Now from a theoretical point of view, it is clear that the number of reactions can be determined by counting the number of non zero singular values of U U T .

Practical implementation

In practice, the ideal case presented above is perturbed for three main reasons:

-The reaction network that we are looking for is a first approximation of chemical or biochemical reactions which can be very complex. The "true" matrix K is probably much larger. The reactions that are fast or of low magnitude can be considered as perturbations of a dominant low dimensional reaction network that we are actually trying to estimate -The measurements are corrupted by noise. This noise can be very important, especially for the measurement of biological quantities for which reliable sensors are not available. -In order to compute u(t) we need a numerical implementation of the filter G(s). Moreover an interpolation is often required to estimate the values of ξ(t i ) and v(t i ) at the same time instants t i . These processes generate additional perturbations.

Data normalisation

In order to avoid conditioning problems and to give the same weighting to all the variables, the data vectors u(t i ) are normalised as follows:

ũi (t j ) = u i (t j ) -a(u i ) √ N s(u i )
where a(u i ) is the average value of the u i (t k ) for k ∈ {1..N }, and s(u i ) their standard deviation.

Practical determination of the number of reactions

In practice, for the reasons we have mentioned above, it is well known that there are no zero eigenvalues for the matrix M = U U T .

The question is then to determine the number of eigenvectors that must be taken into account in order to produce a reasonable approximation of the data u(t).

To answer that question, let us remark that the eigenvalues σ i of M correspond to the variance associated with the corresponding eigenvector (inertia axis) [START_REF] Johnson | Applied multivariate statistical analysis[END_REF].

The method then consists in selecting the p first principal axis which represent a total variance larger than a fixed confidence threshold.

For instance, in the next example, we will consider a threshold (depending on the information available on noise measurements) at 95% of the variance. This leads to the selection of 3 axis, and therefore p = 3.

Remark: if rank (M ) = n it means that rank(K) ≥ n. In such a case we cannot estimate p and measurements of additional variables are requested in order to apply the method presented here.

Parameter Value Unit c0 0.5 g/l.day -1 c1 3 day -1 c2 1 g/l c3 0.2 g/l c4 20 g.day -1 l -1 c5 1 g/l c6 0.2 g/l c7 2 g 2 /l 2 c8 2 g/l c9 0.2 g/l c10 5 day -1 c11 15 g/l c12 5 day -1 c13 0.5 g/l
Table 1 Parameters values.

Example: competitive growth on two substrates

We come back to the example which has been introduced above. Consistently with Equation (1), the model for the state ξ = (S 1 , S 2 , S 3 , E, X, 0, P ) t involving 3 main reactions can thus be written:

dξ dt = K   r 1 r 2 r 3   + v(t)
where

v(t) = D(ξ in -ξ) -Q(ξ)
with ξ in = (S 1in , S 2in , S 3in , 0, 0, 0, 0) t the vector of influent concentrations and Q(ξ) = (0, 0, 0, 0, 0, qo 2 (O), qco 2 (P )) t the vector of gaseous flow rates.

The matrix K was chosen as follows:

K =           -3 0 0 1 -5 0 0.3 0 -0.5 0 0 0.2 0 1 1 0 -2 -1 0 0.3 1.5          
For the simulation purpose, we assume that the kinetics of the three reactions are given by the following expressions. They have not been selected on a realistic basis, but more in order to illustrate our approach on a broad variety of kinetics:

r 1 (S 1 , E) = c 0 S1 S1+c8 E E+c9 X r 2 (S 2 , O, X) = c 1 S2 S2+c2 O O+c3 X r 3 (S 2 , S 3 , O) = c 4 S3 (S3+c5)(S2+c6) O 2 O 2 +c7 X;
The transfer between liquid and gaseous phase is represented by the classical Henry's law: The values of the coefficients c i can be found in Table 1. A 30 day run of the model has been performed using the initial conditions provided in Table 2. The collected data have been corrupted with a white noise of high magnitude (30% of the standard deviation of each component) and sampled. Finally 380 data points are available.

Initial Value Unit condition S1(0) 10 g/l S2(0) 0 g/l S3(0) 5 g/l E(0) 5 g/l X(0) 15 g/l O(0) 0 g/l P (0) 0 g/l
The data (after sampling) are presented in Figure 2. The state variables S 2 , S 3 , E, X, P , O and of the gaseous flow rates q O2 and q CO2 have been measured. We assume here that the state variable S 1 was not recorded in order to illustrate the fact that our approach is applicable even if the full set of state variables is not available for measurement. Moreover the dilution rate and the substrate inflow rate (see Figure 1) have been selected in order to guarantee that the system is sufficiently excited and therefore that the recorded signals will have a sufficiently informative content to expect good identification results. The vectors u(t i ) are then computed by applying a simple moving average from these data and subsequently normalised as explained before. Finally, the eigenvectors of U U T are computed. Figure 3 represents the cumulated variance associated with the number of considered inertia axis. For instance, we can see that two reactions are sufficient to explain 82% of the observed variance. Since three reactions explain 95% of the total variance, it seems reasonable in this example to use 3 reactions for the model.

The reader is referred to [START_REF] Bernard | On the estimation of the pseudo-stoichiometric matrix for mass balance modeling of biotechnological processes[END_REF] for an application to real data, for growth and vanillin production by cultures of the fungus Pycnoporus cinnabarinus in bioreactors.

3 Estimation of the pseudo-stoichiometric matrix K Since we have a value for the number of involved reactions, we are in a position to start the estimation of the (totally or partially) unknown matrix K.

Determination of ImK

Let us use Property 1 which states that ImK is spanned by the eigenvectors ρ i associated with the non zero eigenvalues of U U T . Now, from the experimental data collected through the matrix U U T we get p eigenvectors ρ i that span K. It means that each column k i of K is a linear combination of the ρ i . In other terms, there exists a p × p matrix G such that

K = ρ G
where the columns of matrix ρ are the eigenvectors ρ j . In other words, the family of possible PS matrices K is parameterised by G.

Remark: In general, since the reaction rates are unknown, matrix G (and therefore matrix K) is not identifiable: this can be easily understood on a very simple example. If r 1 (ξ) and r 2 (ξ) are two reaction rates, the term Kr(ξ) can be written:

Kr(ξ) = k 1 r 1 (ξ) + k 2 r 2 (ξ) = k1+k2 2 (r 1 (ξ) + r 2 (ξ)) + k1-k2 2 (r 1 (ξ) -r 2 (ξ))
And therefore both matrices

K = [k 1 k 2 ] and K = [ k1+k2 2 k1-k2
2 ] can produce the same result. The reaction rates associated with the second matrix are then: r1 (ξ) = r 1 (ξ) + r 2 (ξ) and r2 (ξ) = r 1 (ξ)r 2 (ξ).

Additional hypotheses

In order to make matrix G (and K) uniquely identifiable, we need to introduce additional structural constraints. At this stage, all the a priori knowledge on the reaction network should be considered to improve the estimation process.

Normalisation

First, we shall impose (without loss of generality) that each reaction rate is normalised with respect to one species, and therefore that each column of matrix K contains one +1 or one -1. This induces obviously additional constraints on the possible matrices G. Note that sometimes we may not know the sign of the element: the two possible cases must then be considered.

Physical assumptions

One can impose the conservation of elementary mass balances. For example if one wants the carbon to be conserved in the model, if c i is the carbon content of one unit of the state ξ i , it means that we should have for each of the p reactions (j ∈ {1, . . . , p}):

n i=1 c i k ij = 0 (5)
Note however that for macroscopic mass balance where the state variables represent a collection of compounds, the carbon content of the variable can be undetermined. However it can be bounded:

c - i ≤ c i ≤ c + i .
Then equation ( 5) becomes an inequality.

An inequality can also be obtained if we assume that some of the products are not measured, for example if we have a loss of carbon through unmeasured products, we get:

n i=1 c i k ij ≥ 0 3.2.

Biological and biochemical assumptions

When additional constraints are still necessary, we use biochemical assumptions.

When only a subset of the components are present in the reaction at the initial time, the production or not of new components with consumption or not of substrates is an indicator of the variables that are necessary for the reaction. It is clear for example that the first reaction will involve only the substrates which were present at the beginning of the fermentation.

We can for example deduce from this analysis that a specific component is not involved in one of the p reactions and therefore impose a zero in matrix K.

Other assumptions

One can also try to find a matrix K involving the minimum number of components in each reaction (i.e. containing the maximum number of zeros). If these hypotheses are not sufficient, several matrices K can then be identified, parameterised by some parameter, and their biochemical meaning must then be assessed.

Validation

The main result provided by the previous analysis is the determination of the variables which are substrates or products in the reactions or, in other words, the obtained signs of the entries of K.

Another expected result can be the determination of the variables which are not involved in a reaction, corresponding to zero elements in the matrix K. However it is actually very unlikely that the analysis will provide estimates of the elements of K which are exactly zero. The idea consist then in replacing the very small elements by zeros, and to validate the corresponding reaction network using the techniques presented in [START_REF] Bernard | Structural identification of nonlinear mathematical models for bioprocesses[END_REF][START_REF] Bernard | On the estimation of the pseudo-stoichiometric matrix for mass balance modeling of biotechnological processes[END_REF]. These methods are based on the 5 following steps (see [START_REF] Bernard | On the estimation of the pseudo-stoichiometric matrix for mass balance modeling of biotechnological processes[END_REF] for more details):

1. Determination of the vectors which are in the left Kernel of K, i.e. the vectors λ ∈ R n such that:

λ t K = 0 (6)
2. Determination of the associated regression, which is directly deduced by left multiplying Equation ( 3) by λ t . Using Equation ( 6) it gives:

λ i u i (t) = 0 (7)
3. Verification that the u i (t) involved in [START_REF] Chevalier | Toward a systematic determination of complex reaction mechanisms[END_REF] are not related by any other linear relation associated to another left Kernel vector λ (soundness of λ) 4. Computation of the u i (t) from the available data and test of the significance of regression [START_REF] Chevalier | Toward a systematic determination of complex reaction mechanisms[END_REF]. 5. Test of the positivity of the PS coefficients identified in the previous step.

This validation procedure will be illustrated in the example study.

Example (continued)

3.4.1 Statement of the problem and considered data We shall now illustrate the proposed approach with the simulation study of lipase production from olive oil. From the previous study of the number of reactions, we know that 3 reactions should be considered.

We assume here that the first reaction is known, and therefore we only focus on the two other reactions. We are thus in the process of estimating the submatrix K extracted from K by removing the first line and the first column.

A set of noisy data of the state variables S 2 , S 3 , E, X, P , O and of the gaseous flow rates q O2 and q CO2 is produced by simulation as described in Section 2. The goal is to determine the 6×2 matrix K from this data set. More specifically, a question that we want to address is to determine, from the data, which of the two reactions produces the enzyme E.

K identified K identified K after validation         -5 0 0 -0.5 0 0.2 1 1 -2 -1 0.3 1.5                 -3.
Using Equations ( 8) and ( 9) with the obtained values of ρ 1 and ρ 2 , we can now write matrix K parametrised by α 11 and α 22 as follows: In other words this means that S 2 is the only substrate of one reaction and that S 3 is the only substrate of the other one. Thus we will impose k12 = 0 and k21 = 0.

K =         -1.42α
These additional constraints allows us to compute α 11 (0.621) and α 22 (0.93).

Finally we end up with an estimate of matrix K (see Table 3). It is worth noting that the identified matrix K is close to the true one. The value of the (theoretically zero) coefficient k13 is 0.01 which can be neglected with respect to the other coefficients of K. Hence, the unknown part of the structure of matrix K has been recognised. Moreover the estimates of the non-zero entries of the matrix K are quite accurate.

Validation

Here we will validate the identified structure for K with respect to the available data. As it was shown in the previous step, the following structure for matrix K has been identified:

K =         -k11 0 0 -k22 0 k32 1 1 -k51 -k52 k61 k62        
Now the kernel of K T is spanned by the following 4 vectors:

λ1 =         0 k32 k22      k61-k62 k11 0 0 -k62 0 1        
The associated regressions are the following: R 1 : k32 k22 u 3 (t) + u 4 (t) = 0 R 2 : k52-k51 k11 u 2 (t) + k52 u 5 (t) + u 6 (t) = 0 R 3 : -k32 k11 u 2 (t) + u 4 (t) -k32 u 5 (t) = 0 R 4 : k61-k62 k11 u 2 (t) -k62 u 5 (t) + u 7 (t) = 0 [START_REF] Horn | Matrix analysis[END_REF] Note that these regressions are sound [START_REF] Bernard | On the estimation of the pseudo-stoichiometric matrix for mass balance modeling of biotechnological processes[END_REF] in the sense that they do not involve a set of components that are related together by another linear relationship.

The numerical results obtained from the considered regressions are presented in Table 4. It results that all the regressions (10) are highly significant, showing that the estimated reaction network is validated.

Moreover, the following quantities are estimated in Table 4:

{ k32 k22 , k52 -k51 k11 , k52 , k32 k11 , k32 , k61 -k62 k11 , k62 }
It is easy to compute the values of k11 to k62 from this set, leading to the estimate of matrix K proposed in Table 3. The final step consists in verifying that the estimates of the PS coefficients are all positive. This concludes the validation procedure.

Conclusion

Determining a macroscopic reaction network for a bioprocesses is a difficult issue mainly because of the complexity inherent to biological systems. This problem is fundamentally ill stated since the Pseudo-stoichiometric matrix K is generally not identifiable from a data set. We show in this paper how to identify the space generated by the columns of K and how to add constraints in order to determine a unique (or a set of) matrix K.

Through the studied example we have demonstrated that the proposed method can accurately estimate the values of the PS coefficients in spite of noises due to measurements and low sampling frequency.

It is worth noting that this approach does not necessarily require the availability of all the state variables ξ i measurements. Of course, if the measurement of the i th biochemical component ξ i is not available the i th line of matrix K cannot be determined by the method.

qco 2 (

 2 P ) = c 10 (Pc 11 ) and qo 2 (O) = c 12 (Oc 13 )

Fig. 1

 1 Fig. 1 Values of dilution rate, S1in, S2in and S3in used for the simulation.

Fig. 2

 2 Fig.2Experiment simulated from the kinetic modelling corrupted with an additive white noise.

Fig. 3

 3 Fig.3Total variance explained with respect to the number of reactions for the production of lipase from olive oil by Candida rugosa.

3. 4 . 2

 42 Estimation of K Now, using a moving average, we can compute the quantities U i associated with the 6 state variables. Next we compute the matrix M = U T U . The eigenvectors ρ i associated with the two largest eigenvalues are then the basis of ImK. Since G is a 2×2 matrix, the columns k1 and k2 of K can be written:k1 = α 11 ρ 1 + α 12 ρ 2 and k2 = α 21 ρ 1 + α 22 ρ 2 (8)Now we proceed in two successive steps: i. Normalisation. The PS coefficients associated with the biomass growth are normalised : k41 = 1 and k42 = 1. We get then: k41 = 1 = α 11 ρ 41 + α 12 ρ 42 k42 = 1 = α 21 ρ 41 + α 22 ρ 42

Table 2

 2 Initial conditions used for the simulation.

Table 3

 3 True coefficients of matrix K and identified values.

		54 0			-4.54 0	
	0 -0.51			0 -0.50	
	0.01 0.22	 	 	0	0.19	 
	1	1	 	 	1	1	 
	-1.34 -0.87	 	 	-1.33 -0.72	 
	0.18 1.51			0.34 1.24

  Hypothesis: A reaction still takes place when only S 2 [resp. S 3 ] is present at the initial time, and no S 3 [resp. S 2 ] is produced.

	11 -2.65 0.2α 11 -0.13 -0.08α 11 + 0.062 -0.071α 22 + 0.28 -1.2α 22 + 1.12 0.17α 22 -0.67 1 1	     
	-0.19α 11 -1.2 -0.53α 11 + 0.51 -0.45α 22 + 1.93 -0.16α 22 -0.72	 
	ii. Biological hypotheses. Now to determine uniquely matrix K two additional
	assumptions must be introduced.

Table 4

 4 Significance (threshold 5%) of the regressions[START_REF] Horn | Matrix analysis[END_REF] and parameter values.

		Significance Unknown Value	Interval
	R1	YES	k32 k22	2.36	[2.06 2.66]
	R2	YES	k52 -k51 k11 k52	-0.14 0.72	[-0.18 0.09] [0.66 0.77 ]
	R3	YES	k32 k11 -k32	-0.04 -0.2	[ -0.05 -0.03] [ -0.21 -0.18]
	R4	YES	k61 -k62 k11 -k62	-0.35 -1.25	[ -0.42 -0.29 ] [-1.33 -1.16]

Acknowledgement: This work has been carried out with the support provided by the European commission, Information Society Technologies programme, Key action I Systems & Services for the Citizen, contract TELEMAC number IST-2000-28256. It also presents research results of the Belgian Programme on Inter-University Poles of Attraction initiated by the Belgian State, Prime Minister's office for Science, Technology and Culture. The scientific responsibility rests with its authors.


B Chachuat
email: bchachua@mit.edu

O Bernard

Probabilistic Observers for a Class of Uncertain Biological Processes

Keywords: observers, unknown inputs, uncertain parameters, bioprocesses, wastewater treatment, anaerobic digestion

In this paper, probabilistic observers are considered for a class of continuous biological processes described by mass-balance-based models. It is assumed that the probability density functions (PDFs) of the uncertain parameters and inputs of the model, as well as the PDFs of the missing initial conditions are known. Then, the PDFs of the unmeasured state variables are obtained, at any time, by considering the image of these initial PDFs by the flow of the dynamic model (differential system). In comparison to classical open-loop asymptotic and interval observers, the method provides information on the confidence level of the estimates rather than simple upper and lower bounds. Moreover, unlike Kalman filters, probabilistic observers are not restricted to Gaussian distributions for the uncertain parameters. The design and application of a probabilistic observer to an industrial wastewater treatment plant is presented. Finally, a number of practical considerations is discussed in connection to both implementation and utilization issues.

INTRODUCTION

One of the main difficulties in the monitoring and control of biological reactors lies in the absence, in most applications, of cheap and reliable sensors capable of providing direct, online measurements of the biological state variables. The design of state observers, also called software sensors, for the online monitoring of the state variables which are not measurable in real time has thus received much attention in the literature, and has given rise to many practical applications [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF][START_REF] Bernard | State estimation for bioprocesses[END_REF].

Biological processes typically exhibit large uncertainties, which make it difficult to apply classical observation techniques. This is the case in particular for the class of high gain observers [START_REF] Gauthier | A simple observer for nonlinear systems -Applications to bioreactors[END_REF][START_REF] Gauthier | Observability and observers for nonlinear systems[END_REF] that generally lack of robustness when dealing with noisy and/or low frequency measurements, as well as parameter uncertainties. More reliable estimates can be obtained by combining the aforementioned high gain observers with Kalman filtering techniques [START_REF] Deza | Exponential observers for nonlinear systems[END_REF][START_REF] Busvelle | High-gain and non-high-gain observers for nonlinear systems[END_REF], thus yielding extended Kalman filters. Such observers have demonstrated improved ability to cope with measurement noise, and successful applications have even been reported for the estimation of biological systems [START_REF] Bernard | Nonlinear observers for a class of biological systems. Application to validation of a phytoplanktonic growth model[END_REF]. However, the large uncertainties resulting from model mismatch as well as the inherent variability of biological processes are not really handled by these approaches that are sensitive to the picking phenomenon [START_REF] Sussmann | The peaking phenomenon and the global stabilization of nonlinear systems[END_REF]. This may also explain why practical applications of high gain observers to full-scale biological processes are scarce.

For biological processes, another approach has been popularized by Bastin and Dochain [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]. It relies on the underlying structure of many bioprocess models that consist of two parts: (i) a linear part based on mass-balance considerations; and (ii) a nonlinear term which describes the biological reaction rates (kinetics). It should be noted that this latter term refers to the core of the biological activity and is often very sensitive to medium growth conditions; moreover, the associated kinetic parameters are not known accurately in practice. These considerations have motivated the approach proposed by [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] that takes advantage of a linear change of variable to eliminate the uncertain kinetic terms from the observer design; the resulting observers were termed asymptotic observers, and have been widely used in practice [START_REF] Bastin | Nonlinear and adaptive control in biotechnology: a tutorial[END_REF][START_REF] Bernard | State estimation for bioprocesses[END_REF]. Note that the idea of performing the estimation without knowing the reaction rates has its roots in the notion of reaction invariants [START_REF] Fjeld | On a pitfall in stability analysis of chemical reactions[END_REF][START_REF] Asbjørnsen | Response modes of continuous stirred tank reactors[END_REF][START_REF] Asbjørnsen | Reaction invariants in the control of continuous chemical reactors[END_REF] and constitutes an extension of mass-conservation principles. Since they are mainly based on mass-conservation principles (assuming known and constant biological stoichiometry), asymptotic observers are most efficient when the influent composition is well known. In the counterpart however, the corresponding state predictions may converge very slowly; moreover, the estimates may be corrupted in the presence of large measurement noise. In response to this, some authors have proposed a way of combining the convergence and filtering advantages of the high gain approach together with the robustness of the asymptotic observers. For such observers, the state estimates first reach a neighborhood of the unmeasured states very rapidly, from which neighborhood the estimates then converge more slowly, but also more accurately, towards the true values of the state variables [START_REF] Bogaerts | A hybrid asymptotic-kalman observer for bioprocesses[END_REF][START_REF] Lemesle | An hybrid bounded error observer for uncertain bioreactor models[END_REF].

More recently, interval observers have been proposed as a rigorous way of dealing with uncertain processes [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF][START_REF] Rapaport | Parallelotopic and practical observers for non-linear uncertain systems[END_REF]. The idea is to estimate bounds enclosing the state variables, provided that conservative bounds are known for the uncertainties (influent composition, model parameters, etc.). The derivation of interval observers requires that the observation principle be weakened in the sense that the observation error is no longer expected to converge to zero asymptotically. Several successful applications have been reported in the field of biological systems (e.g., [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF]). But although interval observers are more robust than the aforementioned asymptotic observers, still limited information is obtained with respect to the unmeasured state variables, since "only" an upper bound and a lower bound are derived and these bounds may be weak.

In many situations however, more information than simple upper and lower bounds is known for the uncertainties. In particular, an estimate of their probability density functions (PDFs) is often available, e.g., when the uncertain parameters have been obtained from the application of a parameter identification procedure. In this work, we propose to exploit this additional information, in conjunction to the two-part structure of mass-balance base bioprocess models, for estimating the PDFs of the unmeasured state variables. Our approach is based on the same reaction invariants [START_REF] Asbjørnsen | Reaction invariants in the control of continuous chemical reactors[END_REF] as those used by [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] to formulate the asymptotic observers. We first estimate the PDFs of a set of auxiliary variables, and then deduce the PDFs of the unmeasured part of the state by combining these former PDFs with those of the measured state variables. Loosely speaking, these observers can therefore be seen as a hybrid combination between a (linear) Kalman filter and an asymptotic observer. Moreover, because they provide an estimate of the PDFs of the state variables, the term probabilistic observer has been coined for these observers.

The remainder of the paper is organized as follow. Some results on mass-balance bioprocess modeling and asymptotic observers are first recalled in section 2. Probabilistic observers are then described in section 3 for a class of continuous bioprocesses, and their numerical implementation is closely considered. The method is illustrated in section 4 with the application of probabilistic observer to the estimation of the substrate concentration in an industrial wastewater treatment plant. Finally, a number of practical considerations regarding the use of probabilitic observers are discussed in section 5.

PRELIMINARIES

Continuous bioprocesses that are operated in stirred tank bioreactors are often described by means of general mass-balance-based models of the following form [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]:

ξ = Kr (ξ) + D ξ in -ξ -q (ξ) (1)
where ξ ∈ R n ξ (resp. ξ in ) denotes the concentration vector in the liquid phase (resp. in the influent), K ∈ R n ξ ×nr the yield coefficient matrix, r ∈ R nr the reaction rate vector, q ∈ R n ξ the gaseous exchange vector, and D ∈ R the dilution rate.

Throughout the paper, we assume that (n ξn b ) components of the state variables are measured on-line and we denote by ξ a ∈ R n ξ -n b and ξ b ∈ R n b the measured and unmeasured state variables, respectively. We also assumed that the gaseous flow rates q are measured on-line. It follows that Eq. (1) can be rephrased as:

ξa = K a r (ξ) + D ξ in a -ξ a -q a (ξ) ξb = K b r (ξ) + D ξ in b -ξ b -q b (ξ)
where (K a , K b ), ξ in a , ξ in b and (q a , q b ) are the induced partition of K, ξ in and q, respectively.

We now state the main Assumption.

Assumption 1. Matrix K a ∈ R (n ξ -n b )×nr is a full rank column matrix (i.e. with n ξ -n b ≥ n r ).
If Assumption 1 holds, the following n b -dimensional auxiliary system can be derived from a linear change of variables:

ż = -D z -z in -Pq a -q b ( 2
)
ξ b = z -Pξ a (3)
where:

P := -K b K †
and K † a is a left inverse of K a (such an inverse is guaranteed to exist from Assumption 1). It should be noted that the measured state variables ξ a play the role of inputs in the auxiliary system [START_REF] Bernard | State estimation for bioprocesses[END_REF][START_REF] Gauthier | A simple observer for nonlinear systems -Applications to bioreactors[END_REF]. Note also that the auxiliary system is a consequence of both stoichiometric and mass-conservation arguments, and that it corresponds to the reaction invariants [START_REF] Asbjørnsen | Reaction invariants in the control of continuous chemical reactors[END_REF].

An asymptotic observer for the unmeasured variables is obtained by simply making a copy of Eqs. [START_REF] Bernard | State estimation for bioprocesses[END_REF][START_REF] Gauthier | A simple observer for nonlinear systems -Applications to bioreactors[END_REF], provided that a condition of persistence of excitation holds for the input variable D [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]. Such asymptotic observers have however several known defects due to their open-loop nature. In particular, their derivation relies on the implicit assumption that the mass-balance part of the model is known perfectly and that neither the measurements, nor the values of the feeding inputs, nor the estimates of the yield coefficients (matrix K) are biased. Otherwise the predictions of the mass-balance-based observer would be corrupted and may provide very loose estimates of the true concentration values ξ b .

For large uncertainties in particular, it is no longer possible to build exact observers guaranteeing that the observation errors converges to zero. Therefore, the observation principle must be revisited and the results must be weakened. In the next section, we propose new probabilistic observers that estimate the PDF of the unknown state variables ξ b .

PROBABILISTIC OBSERVERS

As mentioned earlier in the paper, biological processes generally exhibit very large uncertainties in both the influent composition and the model parameters. This is particularly tricky for wastewater treatment processes where sensors that are capable of measuring the organic pollutant concentrations accurately and in real-time are still lacking. But although the influent composition is generally not known precisely, still some information is available for it in most practical applications. For the model parameters on the other hand, estimates are often obtained on application of a parameter identification procedure. PDFs (or at least good estimates of them) can therefore be known for either kinds of uncertainties. However, such PDFs can be of any type (especially for the influent concentrations), i.e., they are not necessarily normally distributed.

Mathematical Background

We consider continuous bioprocesses described by mass-balance-based bioprocess models of the form (1). For such models, major sources of uncertainties arise from:

• the unmeasured state variables at initial time ξ b0 ;

• the unmeasured influent concentrations ξ in ; • the unknown model parameters k (e.g., the coefficients of the yield matrix K);

• the noise ǫ in the measurements of the state variables ξ a and of the gaseous flow rates q a , q b .

In the sequel, we make a distinction between the uncertainties originating from the initial conditions, and those corresponding to either influent concentrations, model parameters and measurement noises. These latter uncertainties are gathered within a unique vector p = ξ in , k, ǫ .

In order to calculate the PDF of the unmeasured state variables ξ b , we consider the auxiliary linear dynamical system (2,3) defined in section 2. It is worth reminding that this system requires that Assumption 1 holds and results from a linear change of variables that utilizes the available state measurements ξ a in order to eliminate the uncertain kinetic rates r (ξ). Here, we rewrite the auxiliary system (2,3) in a somewhat more general form, thus yielding the following linear differential system:

   ż (t) = A (p, t) z (t) + b (p, t) ξ b (t) = z (t) + c (p, t) z (0) = ξ b0 -c (p, 0) (S)
where p ∈ P ⊂ R np denotes the vector of the uncertain parameters; ξ b0 ∈ X 0 ⊂ R n b are the unknown initial conditions; and the coefficients of

A ∈ R n b ×n b , b ∈ R n b and c ∈ R n b are C 1 functions of the parameters p.
Assumption 2. The uncertain parameters p in system (S) are constant with known PDFs (stationary disturbances). In addition, we suppose that the PDFs of p and of the initial conditions ξ 0 are independent.

It should be noted that Assumption 2, although not fundamentally necessary, considerably simplifies the synthesis of a probabilistic observer for (S). In the more general case where p is piecewise constant however, all subsequent developments extend in a straightforward manner.

Notations 1. We use the following notations for the random variables: X 0 = (ξ b0 , p), Z 0 = (z 0 , p), X t = (ξ b (t) , p) and Z t = (z (t) , p). In addition, we denote by f X0 , f Z0 , f Xt and f Zt their respective PDFs. Likewise, the PDFs of the unknown state variable ξ b at time t is denoted by f ξ b t .

From classical linear systems theory, the solution ξ b of (S) is given by:

ξ b (t) = ψ t • ϕ t • ψ -1
0 (ξ b0 ) where ψ t and ϕ t are the linear diffeomorphisms defined as:

ψ t : R n b -→ R n b z -→ z + c (p, t) ,
and

ϕ t : R n b -→ R n b z -→ Φ (p, 0, t) z + t 0 Φ (p, τ, t) b (p, τ ) dτ.
Here, Φ (p, •, •) ∈ R n b ×n b denotes the transition matrix of (S), i.e., the solution of the matrix differential equation:

d dt Φ (p, τ, t) = A (p, t) Φ (p, τ, t) Φ (p, τ, τ ) = I n b ×n b .
For the calculation of the PDF of the unmeasured state variables ξ b (t), the idea is to take advantage of the aforementioned analytical relationship between the initial condition ξ b0 , the parameters p and the available measurements (contained in both terms b (p, t) and c (p, t)). This is formalized in the algorithm described in Proposition 1 below. Proposition 1. Under Assumption 2 and given the PDFs of the initial conditions, parameters and measurements, an estimate of PDF f ξ bt of the unmeasured state variables at a given time t > 0 can be computed from the following three-step algorithm:

1. Estimate the PDF f Z0 of the random variable Z 0 at initial time, from the joint PDF X 0 of both the unmeasured state variables ξ b0 and the parameters p. 2. Compute the image f Zt , at time t, of the PDF f Z0 by the flow of the differential equation ż = Az + b in the auxiliary system (S). 3. Estimate the PDF f Xt of the random variable X t , at time t, from the PDF of Z t .

Steps 1 and 3 rely on the diffeomorphisms ψ -1 0 and ψ t , respectively, with ψ t being defined as:

ψ t : V t -→ U t Z t -→ X t = ξ b (t) = ψ t (z (t) , p) , p .
On the other hand, step 2 is based on the diffeomorphism ϕ t defined as:

ϕ t : V 0 -→ V t Z 0 -→ Z t = (z (t) = ϕ t (z 0 , p) , p) . Moreover, we have U 0 ⊃ X 0 × P, V 0 = ψ -1 t (U 0 ), V t = ϕ t (V 0 ) and U t = ψ t (V t ).
Overall, the PDF f ξ bt corresponding to the image ξ b (t) of the random variables ξ b0 ∈ X 0 , for p ∈ P, is given by:

f ξ bt (x) =      P f X0 ψ 0 • ϕ -1 t • ψ -1 t (x, p) exp t 0 Tr [A (p, τ )] dτ dp if x ∈ ψ t • ϕ t • ψ -1 0 (X 0 , P) 0 otherwise. (4)
Remark 1. In the particular case of a standard bioprocess, the matrix A is simply defined as:

A (p, t) = -D (t) I n b ×n b -see Eq.
(2). Accordingly, the following simplification holds:

Tr [A (p, t)] = -n b × D (t) .
In order to prove Proposition 1, we first recall the following result (see, e.g., Corollary 4.7.4 in [START_REF] Grimmett | Probability and Random Processes[END_REF]), which provides a general framework for calculating the image of a given PDF by a C 1 -diffeomorphism.

Theorem 1. Let U be a subset of R n , and let φ : U → V = φ(U) be a C 1 diffeomorphism. If Z is a random variable with a PDF f Z , then the random variable W = φ (Z) has a PDF f W given by:

f W (ω) = f Z φ -1 (ω) det J φ -1 (ω) if ω ∈ V 0 otherwise
where J φ denotes the Jacobian matrix of φ.

Proof of Proposition 1. It is clear, from the analytical solution of (S), that ϕ t and ψ t are C 1 -diffeomorphisms. Therefore, the overall transformation

ψ 0 • ϕ -1 t • ψ -1 t is itself a C 1 - diffeomorphism (chain
, det J (ψ0•ϕ -1 t •ψ -1 t ) = det J ϕ -1 t . Moreover, det J ϕt = det Φ (p, 0, t) ⋆ 0 np×n b I np×np = det Φ (p, 0, t)
and from the Jacobi-Liouville formula, we obtain

det J ϕt = exp t 0 Tr [A (p, τ )]dτ .
Finally, the PDF of the unmeasured state variables ξ b is obtained, at any time t > 0, by integrating the joint state/parameter PDF with respect to the uncertain parameters p ∈ P, which completes the proof.

Numerical Implementation and Practical Application

The numerical implementation of the previous framework can be decomposed into two consecutive steps:

1. Compute the image f Xt of the joint PDF f X0 by considering the diffeomorphism

ψ 0 • ϕ -1 t • ψ -1 t
in Proposition 1. This can be done, e.g., by defining a grid in the joint state/parameter space. Note that due to the linearity of the auxiliary system (S) with respect to z, one only needs to calculate the trajectories of (S) at the vertices of the variation domain of z; the intermediate values are then easily calculated from linear interpolation. 2. Integrate the joint state/parameter PDF f Xt , at a given time t, with respect to the parameters p ∈ P in order to obtain an estimate of the PDF f ξ bt of the unmeasured state variables ξ b (t).

Mention should be made at this point that step 2 is by far the most complex and computationally demanding step in the numerical procedure. Because the influence of the initial condition ξ 0 vanishes as time increases, the joint PDF tends to become degenerate (it becomes infinitely thin), hence giving rise to serious numerical issues. Particular attention must therefore be payed to the integration of the joint PDF on P (see Eq. ( 4)) in order to guarantee accurate results. These issues are illustrated in Example 1 below. Different approaches are then discussed in subsection 3.3, which allow one to address these problems.

Example 1. An illustration of the practical difficulties that are encountered for the numerical implementation of the probabilistic observers is provided through the following simple differential system:

ż (t) = -z (t) + p ∀t > 0, z (0) = z 0 .
We assume that both PDFs of z 0 and p are Gaussian distributions, denoted by N (z 0 , σ z0 ) and N (p, σ p ), respectively. Then, on application of Proposition 1, the joint PDF f Zt of the random variable Z t = (z (t) , p), at a given time t > 0, is given by:

f Zt (ω) = f Z0 ϕ -1 t (ω) exp(-t)
for all ω = (y, q) ∈ R 2 , where

f Z0 (ω) = 1 2πσ z0 σ p e -1 2 » " y-z 0 σz 0 " 2 + " q-p σp " 2 - ϕ -1 t (ω) = e t y + (1 -e t ) q q .
By combining these expressions and letting t go to infinity, it can be easily checked that:

lim t→+∞ f Zt (ω) = +∞ if y = q, lim t→+∞ f Zt (ω) = 0 otherwise.
This is illustrated in Fig. 1.

Practical Measures to Avoid Degenerate PDFs

We start this subsection by noting that the case of interest in this work corresponds to matrix A being Hurwitz, because (S) is intended to be used as an auxiliary system/observer. In this case, the influence of the uncertain initial conditions ξ b0 progressively vanishes with time, thus eventually resulting in a degenerate PDF (see Example 1 above for an illustration). In fact, the diffeomorphism ϕ t itself becomes degenerate for asymptotically stable differential systems, in the sense that it is no longer invertible asymptotically, i.e., lim t→+∞ Φ (p, 0, t) = 0 with • being a matrix norm. In order to avoid such a situation in practice, we consider a small threshold value η > 0. Then, for a given η, ∃t η > 0 such that Φ (p, 0, t) < η ∀t > t η .

For those problems verifying n p ≥ n b , we give the following Proposition that defines an approximate diffeomorphism of ϕ t . Proposition 2. Let ε be a small positive constant, and assume there exists a time t η such that |ϕ t (Z 0 )ϕ ∞ t (p ′ )| ≤ ε for any initial PDF Z 0 ∈ V 0 , where ϕ ∞ t is a diffeomorphism from P ′ (resulting from an appropriate projection of P onto R n b ) to R n b . Then, the solution ξ b (t) of (S) can be approximated as:

ξ b (t) ≈ ψ t • ϕ ∞ t (p ′ ) ∀t > t η .
Remark 2. This proposition is in fact quite general, and the hypotheses are naturally verified for many practical problems. In the case where they are not satisfied, one may still be able to find an appropriate reduced diffeomorphism either by considering a subset of the components of ϕ t only, or by viewing ϕ t as the projection of another diffeomorphism from a space of higher dimension.

By construction, the reduced diffeomorphism ϕ ∞ t is simpler than ϕ t for it only depends on a reduced number of components of the uncertain parameters p, and provides an appropriate approximation of ϕ t for t > t η . From a more practical point of view, ϕ ∞ t is indeed trivial for many applications, as it is the case for instance in Example 1 (see below).

At this point, it becomes clear that a similar procedure as the one described in Proposition 1 can be utilized to calculate the PDF f ξ b , based on the reduced diffeomorphisms ϕ ∞ t and ψ t . Here, the algorithmic procedure is reduced to steps 1 and 2 only, and one finally needs to integrate the joint PDF of X t = (ξ b (t) , p ′ ) on the reduced parameter space P ′ in order to obtain the desired PDF f ξ bt .

Example 1 (continued). On application of Proposition 2, the following trivial diffeomorphism is obtained (the same notation as in Example 1 are used):

ϕ ∞ t (q) = q.
Accordingly, the following approximation is obtained for the PDF of z (t):

f zt (q) = f p ϕ ∞ t -1 (q) = f p (q) ∀q ∈ R.
Note that the PDF of z is then equal to the PDF of p.

An alternative approach that may be used to prevent the PDF f ξ bt from becoming degenerate consists of adding an artificial term that acts as a diffusion [START_REF] Zeeman | Stability of dynamical systems[END_REF]. A simple way of doing it is to replace the transition matrix Φ (0, t, p) in the expression of ϕ t by [Φ (0, t, p) + η × I n b ×n b ], with η ≪ 1. This alternative approach also presents the advantage for the user to keep the same numerical procedure for all t > 0.

In the coming practical case study, note that we do not have to deal with the aforementioned numerical issues. This is because the PDF of the state estimates does not become degenerate within the considered time frame, as the observer is reinitialized on a weekly basis by using new off-line measurements.

INDUSTRIAL CASE STUDY

In this section, we illustrate our approach with the concrete case study of an industrial anaerobic digestion (AD) wastewater treatment plant. AD consists of several consecutive biological degradation steps that take place in the absence of oxygen, each step being performed by different groups of specialized bacteria. Complex organic compounds are decomposed and bioconverted, on the one hand, into biogas (i.e., a mixture of CO 2 and CH 4 ) and, on the other hand, into microbial biomass and residual organic matter.

The main advantages of AD wastewater treatment plants, as compared to the more classical aerobic activated sludge plants, are their ability to handle substrates that are otherwise difficult to degrade at high concentrations, their high efficiency, as well as their low sludge production. Moreover, the produced methane can be utilized as a valuable energy resource. However, many industrials are still reluctant to use AD processes, mainly because AD processes have two stable equilibrium points [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF], one of which corresponds to biomass washout and must be avoided for obvious reasons. Accordingly, AD processes therefore require a strong monitoring procedure in order to operate safely. Moreover, the observer's robustness towards modeling uncertainties (kinetic rates in particular) is of the utmost importance for AD, since still no reliable model has been proposed to date.

Mass-Balance Modeling of the Digester

We use a very simple model that describes the key features of AD [START_REF] Mailleret | Robust nonlinear adaptive control for bioreactors with unknown kinetics[END_REF]. Note that this model only takes the main degradation processes into account and is therefore a rough summary of the overall biological reaction network. It accounts for a single substrate (index of the total organic matter concentration, called Chemical Oxygen Demand (COD)) and a single bacterial population, denoted by s and x, respectively. By assuming perfect mixing in the digester, the following 2-dimensional dynamical model is obtained:

ẋ (t) = -αD (t) x (t) + µ (•) x (t) (5)
ṡ (t) = -D (t) s (t) -s in (t) -k 1 µ (•) x (t) (6)
where s in (t) stands for the concentration of substrate in the feeding stream; k 1 is the yield coefficient associated to biomass production; and α represents the fraction of bacteria being present in the liquid phase remaining fraction (1α) is attached onto a support within the digester). In addition, the methane flow rate q CH 4 is defined as follows:

q CH 4 (t) = k 2 µ (•) x (t) (7)
where k 2 stands for the yield coefficient associated to methane production. It should be noted that the bacterial growth rate, µ (•), is often described as a non-monotonic function of Haldanetype [START_REF] Mailleret | Robust nonlinear adaptive control for bioreactors with unknown kinetics[END_REF] that is responsible for the process instability.

Probabilistic Observer for the COD Concentration

The objective in this application is to design an observer capable of estimating the COD concentration s, based on the methane measurements, and without making any assumption on the growth rate function µ (•). Also, we want to take advantage of the off-line COD measurements that are generally performed at a very low sampling frequency, e.g., once a week. In practice, the observer will therefore run in open-loop mode between two consecutive COD measurements, and will then be reinitialized after a new COD measurement has become available (see subsection 3.3).

Based on the model equations ( 5)- [START_REF] Bernard | Nonlinear observers for a class of biological systems. Application to validation of a phytoplanktonic growth model[END_REF], the following auxiliary system can be derived for the biomass concentration, provided that the methane gaseous flow rate as well as the influent COD are measured on-line (see section 2):

ṡ (t) = -D (t) s (t) + D (t) s in - k 1 k 2 q CH 4 (t) . (8)
A set of experimental data obtained from a 2, 000 m 3 industrial digester located at Agralco Ltd. in Estella, Spain, is considered subsequently. This digester is fed with wine distillery vinasses, and measurements are available for the dilution rate D and for the methane outflow rate q CH 4 , both at a high sampling frequency and for more than 140 days.

In order to design a probabilistic observer that estimates the COD concentration based on Eq. ( 8), the following uncertainties are taken into account:

i) Parameter γ = k1
k2 . The value of γ was estimated off-line based on mass-balance considerations along with its standard deviation: γ = 5.30 10 -2 ± 0.98%. It is assumed that the probability density for this parameter is Gaussian.

ii) Influent concentration s in . The uncertainty on s in is handled by considering a multiplicative noise ε defined as s in (t) = s in (1 + ε (t)) , ∀t. It was found from experimental data that s in = 35.9 g.L -1 and that the PDF of ε can be modeled by a normal distribution N (0, σ ε ) with σ ε = 0.124. iii) measured COD concentration s mes . The observer is reinitialized on a weekly basis by using the available off-line COD measurements. The uncertainty associated to the COD measurements is modeled by a uniform PDF (within the range s mes ± 20%).

On the other hand, the measurements of D and q CH 4 , which are obtained at a very high sampling frequency, have been low-pass filtered and, for the sake of simplicity, no error is taken into account for these measurements.

A probabilistic observer for the COD concentration is derived from Eq. ( 8) on application of Proposition 1. In this case, note that ψ t = ψ -1 0 = I and Φ (τ, t) = exp -t τ D (s) ds . The resulting PDFs of the estimate ŝ are represented in Figure 2 (upper plot) at different time instants. For the sake of observer validation, experimental COD concentrations that were measured on a daily basis are also displayed on this figure; it is worth reminding to the reader that these additional points were not used for designing the probabilistic observer. Note that the variation range of ŝ (corresponding to the vertical segments) is rather conservative, but a significant part of the range has zero probability. Also note that most of the experimental points are enclosed within the non-zero probability part, and that the maximum likelihood estimate is close to the experimental points.

A quite straightforward utilization of the resulting PDFs of ŝ consists of extracting their maximum likelihood values (or their average values) as well as their standard deviations. In practice, it is very likely that only these two informations will be made available to the user. Both are shown in figure 2 (lower plot). It can be seen from this plot that the estimates are close to the experimental measurements and that the trajectories corresponding to the standard deviation of ŝ mainly enclose these measurements. However, the fact that the data are not in the appropriate range between days 43 to 49.5 warrants some further explanations. In fact the off-line measurements of the COD concentration turn out to be corrupted due to chemicals expiration. Interestingly enough, it appears that the probabilistic observer could therefore be used in such situations for efficiently detecting measurement biases or failures.

The interested reader is also referred to [START_REF] Chachuat | Design of probabilitic observers for mass-balance based bioprocess models[END_REF] for discussion of different application of the probabilistic observers in the field of bioprocesses.

PRACTICAL CONSIDERATIONS

The proposed probabilistic observers present some peculiarities in comparison to usual observers. The most important, maybe, lies in the fact that they compute, at any time instant, the exact PDF of the estimated state variables. Accordingly, there is no required convergence period regarding the use of these observers. What can be said concerning convergence, however, is that the estimated state PDF converges to a limit PDF that is independent of the initial state PDF, provided that matrix A in (S) is Hurwitz; moreover, the corresponding rate of convergence is dictated by the smallest eigenvalue of matrix A. This "natural" convergence is closely related to the problem of degenerate PDFs discussed earlier in subsection 3.3. Therefore, probabilistic observers could be equally well seen as detectors, since they are fully based on the (reduced) dynamics of the system. Moreover, these observers are much more efficient when an accurate estimate of the initial state PDF is provided, since convergence to the aforementioned limit PDF is then much faster. This was illustrated by the industrial case study presented previously.

In practice, the information contained in the calculated state PDFs is often too complex to be displayed for an industrial operator. Instead, the end-user can be given the average value together with an estimate of the standard deviation value. An interval enclosing those state values having a cumulated probability equal to a given confidence level (e.g. 95%) can be also be easily calculated.

On the other hand, the information contained in the state PDF can be fed directly into a supervision procedure [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF], hence allowing the supervision system to account for the level of uncertainty when making a decision regarding the state of the process. This is of crucial importance in certain practical applications having multiple equilibrium points, such as the anaerobic digestion process. The ability to detect, early enough during the operation, potential state trajectories outside of the attraction basin of the desired operating point could indeed help prevent the process from dramatic problems.

Finally, since many applications in biotechnology are characterized by a high level of noise, probabilistic observers may produce very widespread PDFs that are useless from a practical point of view. However, determining which source of uncertainty must be reduced in order for the observer to produce a more accurate prediction is quite straightforward. Therefore, the proposed observers could be used as tools to identify which measurements should be improved in priority for producing more accurate state estimates. In such situation, specific effort could then be made on the designated sensors.

CONCLUSIONS

In this paper, the design of probabilistic observers has been investigated for a class of continuous biological processes described by mass-balance-based models. It has been shown how the PDF of the unmeasured state variables can be obtained from the knowledge of the PDF of the uncertain parameters, inputs and measurement noises. The practical implementation of probabilistic observers has also been closely considered. These observers share some common features with both asymptotic observers and classical Kalman filtering techniques, but also differ in several points; the main features of probabilistic observers are summarized hereafter.

• in order to combat the inherent uncertainties characterizing biological systems, probabilistic observers are based on a reduced system obtained by eliminating a number of functions (e.g., microbial kinetics) based on linear algebra manipulations; • in line with asymptotic observers [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF], probabilistic observers are open-loop by nature, in the sense that the estimation is not corrected by the difference between the predicted outputs and the measurements; • unlike classical Kalman filters finally, probabilistic observers can handle any kind of PDF, including uniform and normal PDFs as well as compact support based PDFs characterizing bounded uncertainties.

The design and application of a probabilistic observer for the estimation of the COD concentration in an industrial anaerobic wastewater treatment plant has been presented. This case study demonstrates both the tractability and the good performances of the proposed observers for a practical problem. Finally, a number of general considerations have been discussed regarding potential applications of probabilistic observers in connection to supervision and fault diagnosis of bioprocesses.

  rule of differentiation), and Theorem 1 applies. Because, det J ψ0 = det J ψt = I n b +np Copyright c 2005 John Wiley & Sons, Ltd. Int. J. Robust Nonlinear Control 2005; ##:1-16 Prepared using rncauth.cls

  then

Copyright c 2005 John Wiley & Sons, Ltd. Int. J. Robust Nonlinear Control 2005; ##:1-16 Prepared using rncauth.cls

ACKNOWLEDGEMENTS

This work has been carried out with the support provided by the European Commission, Information Society Technologies program, Key action I Systems & Services for the Citizen, contract TELEMAC number IST-2000-28256. This information is provided under the sole responsibility of the author and does not necessarily represent the opinion of the European Commission, which is not responsible for any use that might be made of it.