Measurement processing for state estimation and fault identification in batch fermentations

Abstract

This work describes an application of maximum likelihood identification and statistical detection techniques for determining the presence and nature of abnormal behaviors in batch fermentations. By appropriately organizing these established techniques, a novel algorithm that is able to detect and isolate faults in nonlinear and uncertain processes was developed. The technique processes residuals from a nonlinear filter based on the assumed model of fermentation. This information is combined with mass balances to conduct statistical tests that are used as the core of the detection procedure. The approach uses a sliding window to capture the present statistical properties of filtering and mass-balance residuals. In order to avoid divergence of the nonlinear monitor filter, the maximum likelihood states and parameters are periodically estimated. The maximum likelihood parameters are used to update the kinetic parameter values of the monitor filter. If the occurrence of a fault is detected, alternative faulty model structures are evaluated statistically through the use of log-likelihood function values and chi2 detection tests. Simulation obtained for xanthan gum batch fermentations are encouraging.

fermentation process; stochastic model; maximum likelihood state

BIOTECHNOLOGY

Measurement processing for state estimation and fault identification in batch fermentations

R. Dondo

Instituto de Desarrollo Tecnológico para la Industria Química, INTEC (CONICET-UNL), Güemes 3450 (3000) Santa Fe, República Argentina. E-mail, [e-mail]

ABSTRACT

This work describes an application of maximum likelihood identification and statistical detection techniques for determining the presence and nature of abnormal behaviors in batch fermentations. By appropriately organizing these established techniques, a novel algorithm that is able to detect and isolate faults in nonlinear and uncertain processes was developed. The technique processes residuals from a nonlinear filter based on the assumed model of fermentation. This information is combined with mass balances to conduct statistical tests that are used as the core of the detection procedure. The approach uses a sliding window to capture the present statistical properties of filtering and mass-balance residuals. In order to avoid divergence of the nonlinear monitor filter, the maximum likelihood states and parameters are periodically estimated. The maximum likelihood parameters are used to update the kinetic parameter values of the monitor filter. If the occurrence of a fault is detected, alternative faulty model structures are evaluated statistically through the use of log-likelihood function values and c² detection tests. Simulation obtained for xanthan gum batch fermentations are encouraging.

Keywords: fermentation process, stochastic model, maximum likelihood state.

INTRODUCTION

The importance of on-line monitoring of biotechnological processes has increased during the last twenty years. Advantages include gaining knowledge about the state of the process and the possibility of detecting and isolating abnormal process developments at early stages. This reduces process costs, contributes to process safety and helps in trouble-shooting and process accommodation. The main problem in fermentation monitoring and control is the fact that process variables usually cannot be measured on-line (e.g., biomass, substrate and product concentrations). Monitoring and controlling these processes can therefore be difficult because only indirect measurements are available on-line, while calculated values may be rather uncertain. This can be due to uncertainty with respect to the equations used, measurement errors or both. For automatic control this may have serious consequences, especially as the actual variables of interest often cannot be directly controlled and related variables are controlled instead.

In fermentation processes, on-line and off-line measurements are the main source of information about the state of the process. In combination with model-based calculations, they are used to produce estimations for monitoring purposes as well as for automatic and manual process control. Model parameters are also established by means of measurement. It is therefore important to have an accurate and consistent set of measurements. In practice, measurement errors will always occur. The most common type is a relatively small random error due to minor disturbances in the measurement equipment. The magnitude of these errors, commonly referred as measurement noise, defines the accuracy of the measurement. They are usually regarded as zero-mean with Gaussian distribution. This kind of noise can be eliminated by the use of state estimators such as Kalman filters. On the other hand, multi rate estimators [Halme, 1987] are observers that are well suited for state estimation in fermentations. In these estimators, the measurement vector is expanded to include infrequent off-line measured variables when these measurements are available. This expansion is only made functional at the time of measurement. To overcome problems with the time delay caused by laboratory analysis the technique uses stored data. The estimates are then recalculated from the time of measurement to the present as soon as the measurement value becomes available. In many cases, there is a certain amount of “overlap” between off-line and on-line measurements. This overlap together with conservation equations provides constraints to improve the accuracy of the measurements and to detect significant errors in the measurements or in the model used by the fermentation observer. Faulty sensors and omitted components can be detected in this way. This results in enhanced reliability and accuracy of on-line state and parameter estimates.

Much research on state estimation in bioprocesses can be found in the literature. Some of the most relevant are by Stephanopoulos and San [1984a and b], Bastin and Dochain [1986, 1990] and Gudi et al. [1994]. Two different detection methods can be used for fault detection in batch fermentations. The first is herein referred to as the “residual-based detection method”. It focuses on the analysis of estimation residuals of a Kalman-filter-type observer. The second is herein referred to as the “balance-based detection method” and it uses conservation principles for testing the consistency of the variables measured. Isermann [1984] and Frank [1990] offer survey of the residual-based detection method. Alcorta García and Frank [1997] reviewed observer based approaches to several classes of deterministic nonlinear systems. Significant works related to the balances-based detection method were published by Wang and Stephanopoulos [1983] and Van der Heijden [1994a and b]. Dondo [2003] proposed the simultaneous use of both methods. The idea behind the use of both methods is that the limitations of one be compensated by the use of the other. In the present work an evolution of the idea developed in Dondo [2003], which is designed for obtaining robust and accurate state estimation and fault diagnostics under parametric uncertainty, is presented.

This work is organized as follows: section 2 discusses the specifics of estimation and detection in batch fermentations. In section 3, a methodology of state estimation and fault detection for batch fermentations is presented. Numerical results are shown in section 4 and the conclusions are outlined in section 5.

PROBLEM DISCUSSION

Instrumentation failures and abrupt kinetic changes can be understood as a deviation of a process variable that is not permitted and that leads to an inability to maintain control of the running fermentation. In the present work, these deviations are generically referred to as faults. Fig. 1 shows a block diagram for fault detection and isolation in fermentations. Checking whether measured and/or unmeasured estimated variables are within a given tolerance of their normal values means detection. If the check is not passed, this leads to a fault message. Tasks related to detection and isolation can be divided into the following stages:

• Residual generation: computation of functions that are sensitive to the occurrence of a fault.

• Fault decision: checking residuals if there is a fault.

• Fault isolation: identification of fault occurrence time, type, size and source.

Computational requirements are a practical problem regarding fault detection because algorithms for detection and diagnosis are often computation intensive. Nevertheless, this is not a problem for the reason that batch-fermentations are generally carried out over many hours or even days. Furthermore, a detection algorithm must have two important capabilities:

• The ability to quickly detect the occurrence of an abnormal event within a short period following its occurrence.

• The ability to correctly identify the event, its occurrence time and its magnitude.

• One of the fundamental aims of supervision of a biotechnological process is to promptly detect and identify abnormal behaviors (faults) in order to take corrective action for maintaining the fermentation running. This capability is crucial for enhancing the reliability of the operating equipment and to ensuring a profitable operation. Examples of sources of faults in batch fermentations are

• A measured variable has a significant error.

• The system description is incorrect because a component has a composition different from that specified or a component is not included in the description of the fermentation.

• Abrupt kinetic changes are produced during the course of the fermentation.

• The assumed measurement variances are incorrect resulting in a poorly tuned estimation algorithm.

• Since detection methods must be sensitive to the occurrence of faults but robust to noises, modeling errors and signal variations, the following trade-off exist [Isermann, 1984]:

• Size of faults vs. detection time.

• Parameter estimation rate vs. false alarm rate.

• Detection time vs. false alarm rate.

Methods that are designed for detection of abrupt changes are usually not suitable for state and parameter estimation and vice versa. These considerations call for developing an innovative approach and have motivated the methodology presented below.

THE ESTIMATION AND DETECTION METHODOLOGY

An aerobic fermentation with production of a single metabolite can be seen as three parallel “chemical reactions” denoted partial metabolisms [Minkevich, 1983]. These reactions are biomass production, metabolite production and main substrate oxidation. Thus, the aerobic growth of biomass (X) consuming a carbon and energy source (S) and an independent nitrogen source that can also contain carbon (S_N) while generating a metabolite P, CO₂, and H₂O can be written as

Biomass production:

Metabolite production:

Main substrate oxidation:

Compositions of components X, S, P and S_N are expressed by their atomic formulae CH_b1O_c1N_d1, C_a2H_b2O_c2N_d2,_, C_a3H_b3O_c3 and C_a4H_b4O_c4N_d4, respectively (the metabolite is assumed to be a nitrogen-free component). The kinetics of each reaction are characterized by the evolution of each one of the relevant reaction components: X, P and S^R (quantity of oxidized main substrate). Expressions (1) can be expressed as

or

where C_{Y I/J} is a matrix of stoichiometric yields Y_I/J and DI is a vector of net production of the system components. Since element balances are constraints that must always be satisfied, they are constraints to be met by the fermentation “reactions”. These balances mean four constraints (one for each element considered: C, H, O and N) to be met by the relation between seven components (X, S, P, S_N, O₂, CO₂ and H₂O). Thus, an aerobic fermentation with formation of a single metabolic product has ( 7-4 ) = 3 degrees of freedom and unknown component evolutions may be obtained from the knowledge of the stoichiometric yields Y_I/J and three component evolutions. Thus, if there are more than three measurements of component evolutions, an overlap of measurements is produced. This overlap and conservation equations (2) provide constraints to improve the accuracy of measurements and to detect significant errors in measurements or in the model used by a fermentation observer. In this way, constraints (2) can be lumped into a conventional detection methodology for building an efficient estimation and detection procedure. To do this, let us assume that the dynamic model of a batch fermentation is represented by the usual nonlinear state-space formulation:

In this formulation, kinetic parameters p appear in the dynamic function f(•), while if some of the fermentation components are measured, balances constraints are in the state-measurement relations c(•) [Dondo and Marqués, 2002]. In order to explain this, let us assume that the fermentation states are the biomass concentration [X], the metabolic product concentration [P] and the amount of main substrate oxidized (DS^R). The amount of oxygen consumed (DO₂) and the amount of carbon dioxide produced (DCO₂) are on-line measurements, and the biomass concentration [X], the metabolic product concentration [P] and the main substrate concentration [S] are off-line measurements. This is probably the most common measurement arrangement in batch fermenters. From expression (2.a) it is clear that the relation between these main components is linear. Thus, the relation between states and measurements is also linear [Dondo and Marqués, 2002], and it can be expressed as the product of a matrix of stoichiometric coefficients C by the vector of states x_t(k):

Since most kinetic models of batch fermentations are nonlinear and have parametric uncertainty, adaptive nonlinear observers are used for monitoring this kind of processes [Gudi et al., 1994; Dondo, 2003]. Although the use of off-line information within the adaptive estimation procedure makes the estimates more robust, the uncertainty of estimated variables (given by the state covariance matrix) is relatively high [Dondo, 2003]. This uncertainty is also manifested as high values of the measurement covariance matrix and reduces the sensitivity of any detection test [Wilsky, 1986]. Hence, hypothesis distinguishibility is rather difficult. Consequently, to promptly detect a fault, a sensitive logic that takes into account the past history of the system as well as parametric uncertainty must be designed. The approach to the estimation and detection problem presented here relies on an intensive use of statistical criterions as indicators of a faulty process. These criterions are based on “signatures” of the fermentation which are monitored and compared with a priori estimations based on the unfaulty model of the system. Statistical indicators are also used to determine the occurrence time of faults and their identity. The operative logic of the estimation-detection procedure is detailed in the following subsections.

Normal Process Operation

Let us assume that a multi rate extended Kalman filter (EKF) is used as fermentation states and parameters estimator (Table 1). The innovation sequence of the filter is defined by

where y_t(k) is the measurement vector value at time t(k) and c(x_t(k/k-1)) is the prediction of the measurement vector value based on prediction of the state vector value x_t(k/k-1). If there are no convergence problems and under the no-fault hypothesis (the model corresponds to the reality and the measurements are unbiased and corrupted by independent sequences of “white noises”), sequence g_t(k), should be a zero-mean V_t(k) covariance sequence, where V_t(k) is defined by

In eq. (6), c_X is the Jacobian matrix of c(x_t(k/k-1)), P_t(k/k-1) is the predicted state-covariance matrix at time t(k) and S is the matrix of noise variances of the measurement vector. Faults and abrupt dynamic changes are usually manifested as unexpected values of y_t(k). Therefore, a non-zero-mean sequence indicates a possible fault. For testing the zero-mean hypothesis, the following statistical indicators are proposed:

In eq. (7.a), l_t(k) represents the sum of normalized squared innovations on time horizon w and m_t(j) is the dimension of the measurement vector at time t(j). In eq. (7.b), lⁱ_t(k) is the sum on time horizon w of the squared innovations of measurement i normalized by their variances v_i². Finally, l^-i_t(k) represents the sum of normalized squared innovations of all but i measurements on time horizon w. The measurement covariance V_t(j) is computed by the EKF. Matrix V^-i_t(j) and the variance nⁱ_t(j)² are extracted from V_t(j). For a selected window size, w, the effect on residuals at times t(j) < t(w) is neglected. Under the no-fault-hypothesis, variables l_t(k), lⁱ_t(k) and l^-i_t(k) should be c² distributed variables with S^w_j=1 m_t(j), w and S_j=1^w (m_t(j)-1) degrees of freedom, respectively. Therefore, by defining thresholds j, jⁱ and j^-i with confidence levels selected a-priori, it is possible to carry out the following tests:

Statistical indicators, l_t(k), lⁱ_t(k) and l^-i_t(k), computed on the w-lag sliding window, provide simple and efficient detection tools. However, since the sliding window involves residuals from nonlinear filters and a limited sample size, actual indicator values will not be exactly c²-distributed. Detection thresholds cited above, which are based on asymptotic properties, should therefore be approximate. Thus, persistence tests (the indicators must exceed their thresholds over a given time period) should be used to cut down false alarms due to spurious and unmodeled events.

In order to avoid convergence problems due to the effect of nonlinearities and to keep variables l_t(k), lⁱ_t(k) and l^-i_t(k) sensitive to occurrence of a fault, the value of P_t(k/k-1) must be keep as small as possible [Dondo, 2003]. For this purpose, each time that there is an off-line measurement available, a maximum likelihood optimization is used in a time window W = t(0),...,t(k):

subject to

This maximization fulfills two tasks: (i) it keeps estimated parameters as near as possible to the true parameter values to avoid divergence of the EKF monitor and (ii) it keeps the covariance matrix V_t(k) as small as possible in order to maintain l_t(k), lⁱ_t(k) and l^-i_t(k) sensitive to minor variations in the innovation sequence. This permits use of small sliding-window lags, w, and then variables l_t(k), lⁱ_t(k) and l^-i_t(k) will be able to quickly react to an unexpected event. The use of small sliding windows is critically necessary because of the use of adaptive observers. This is because effects of unexpected measurement values on innovations are manifested as correction of the estimated parameter values and will promptly disappear. Maximization (9) had been a very difficult task, particularly in the case of nonlinear systems [Young, 1981]. Main difficulties reported in the literature are the need for considerable computational power and the computation of analytical Jacobian and Hessian for the maximization algorithm. Nevertheless, these difficulties have been practically overcome because of the tremendous advances in computational power and the development of efficient minimization methods that do not use Jacobian and Hessian matrixes (Downhill simplex method due to Nelder and Mead and Powell’s method [Press et al., 1992]). Since the maximum likelihood estimation gives the minimum-variance estimates [Mendel, 1995] it is utilized for on-line state and parameter estimation in a specified time-window W = (t₍₀₎ t_(k)) in order to reinitialize the monitor Kalman filter with minimum variance states and maximum likelihood parameters. Time-window lag, W, must be large enough to allow a significant collection of information, but small enough to avoid lumping parameter variations.

On the other hand, if there are redundant state-measurement relations when an off-line measurement is available, the following nonlinear least-squares estimation of states and measurement can be obtained:

If the difference between and x_t(k) is not acceptable, it is possible to re-estimate the state and measurement vectors by recalculating eqs. (11) and (12), replacing x_t(k) by and c(x_t(k)) by , respectively. The procedure can be repeated until no significant modifications of estimates are obtained [Mendel, 1995]. Thus, the following residuals vector and covariance matrix can be defined as follows [Dondo, 2003]:

Now the following statistical indicators can be computed:

where variance s_i² is extracted from variance matrix P_e. Variable h_t(k)^-i is computed using all but the i measurement and P^-i_e is calculated by eq. (14), eliminating columns and rows related to measurement i from matrixes S and c_X. If the measurement arrangement is given by eq. (4.b), the Jacobian c_X is to be replaced by C, and, for tests h_t(k), hⁱ_t(k), and h^-i_t(k), the nonlinear least-squares eqs. (11) and (12) are simplified to the linear case:

If elements of e_t(k) are assumed to be zero-mean and Gaussian-distributed, under the no-fault hypothesis, h_t(k), hⁱ_t(k) and h^-i_t(k) are approximately c²-distributed variables with (n-m), 1 and (n-m-1) degrees of freedom, respectively. Thus, the following tests can be conducted:

In expression (4.b) there are n = 3 states and m = 5 measurements and therefore the degree of redundancy is 2. Furthermore, as there are two on-line measurements (DO₂ and DCO₂), it follows that l_t(k)^O2 = l_t(k)^-CO2 and viceversa.

Variables defined by eqs. (15), when available, and by eqs. (7) form a set of statistical indicators that provide strong indications of the occurrence of a fault. For example, if on-line measurement i is suddenly biased, l_t(k) should indicate the occurrence of an unexpected event, lⁱ_t(k) should show a sharp increase in its value and l^-i_t(k) should be quasi-invariant to this bias. When an off-line measurement is available, indicators h_t(k), hⁱ_t(k) and h^-i_t(k) should also have a similar behavior.

Tests (8) and (18) give an intuitive justification not only for use in fault detection, but also for formulating a detection/diagnosis scheme. This approach consists of including in the extended Kalman filter various possible faulty models, estimating their parameters by a maximum-likelihood approach while testing indicators h_t(k), hⁱ_t(k) and h^-i_t(k) for these models and then choosing the model with the maximum log-likelihood function. This diagnosis scheme will be detailed in the next subsection.

Faulty Process Operation

If a fault is detected, its cause should be identified. Once the information needed to detect and diagnose faults (residuals and measurement history) has been accumulated, it is necessary to interpret the information in various ways: whether or not there is a failure, the probability of occurrence of a failure and the failure most likely to have occurred. Each hypothesis (i.e., sensor drift, formation of a by-product, etc) will demonstrate a specific time-dependent pattern in measurement evolution and tests l_t(k), lⁱ_t(k), l^-i_t(k), h_t(k), hⁱ_t(k) and h^-i_t(k). The idea behind this is that the signature of the measurement evolution contains information on the kind and magnitude of the fault. Thus, every suspected fault characterized by a given type (I), identity (J), magnitude (u) and occurrence time (t) is simulated, and data from these simulations are used to estimate the faulty model states and parameters and to define hypothesis log-likelihood functions L_I,J(u, t). The technique can be viewed as an extension of the generalized likelihood ratio method (GLR) [Wilsky, 1986] to the nonlinear case. The general form of the resulting log-likelihood maximization problem can be written as

subject to

where p denotes the previously estimated parameters of the unfaulty process model and (y_t(j) – c_IJ(x_t(j/j-1),u,t)) denotes residuals from the (I, J) faulty-model filter. Log-likelihood function values L_I,J(u,t) are computed for each alternative fault location and structure and are ranked from largest to smallest to assess the appropriateness of a particular hypothesis about the unknown event. In addition, the evolution of indicators h_t(k), hⁱ_t(k) and h^-i_t(k) provides further information for discriminating between different hypotheses.

a) Detection and Isolation of Significant Measurement Errors (I = 1)

A measurement bias can frequently be found. Then the mean of its measurement noise is different from zero. Sensor drift or inaccurate calibration may cause the bias. This type of error can be disastrous when the measured variable is used to determine another process variable for control purposes and it must be promptly detected. But if on-line measured variable i has a significant error, the i element of the innovation sequence is biased, and therefore the value of l_t(k) should increase, a sudden and large change in the value of lⁱ_t(k) is expected and l^-i_t(k) should remain below its threshold. Furthermore, when an off-line sample is analyzed, h_t(k), hⁱ_t(k), and h^-i_t(k) should show a similar pattern. Thus, the error can be detected by analyzing the behavior of these statistical indicators. Failure parameters are estimated by maximizing the log-likelihood function (19). In the same way, if off-line measured variable i has a significant error, since indicators l_t(k), lⁱ_t(k) and l^-i_t(k) cannot be calculated, the value of h_t(k) should increase, a sudden and large change in the value of hⁱ_t(k) is expected and h^-i_t(k) should remain below its threshold.

If the state-measurement relations of the fermentation observer are given by eq. (4.b), variables l_t(k), l^O2_t(k) and l^CO2_t(k) allow detection of significant errors in on-line measurements DO₂ and DCO₂. Analogously, variables h_t(k), hⁱ_t(k), and h^-i_t(k) allow detection of significant errors in off-line measurements DS, DX and DP. In principle two simultaneous but related measurement errors can be detected. For example, if there is a fault in the gas-flow measurement device, DO₂ and DCO₂ will have a correlated bias. Thus l_t(k), l^O2_t(k) and l^CO2_t(k) indicators will not trigger the alarm but indicators h_t(k), hⁱ_t(k) and h^-i_t(k) will do it. Since residual vectors are usually rather inaccurate, a search of two simultaneous error sources will not give trustworthy results. Nevertheless, simulation of a hypothesis with a simultaneous and correlated drift in DO₂ and DCO₂ should give the maximum log-likelihood-function value.

b) Detection and Isolation of Incorrect System Descriptions (I = 2)

If the evolution of some system components is measured and introduced into eq. (2.b), due to measurement noises this equality will never be exactly satisfied. It can be better written as

where n is a residuals vector. Its expected value under the non-fault hypothesis is 0. But in addition to measurement noises, error in the specified constraint (2) may occasionally be encountered. This may be due to time-varying or ill-defined component composition (i.e., biomass), a component omitted from the balance equation or an alternative metabolic pathway (e.g., partial substrate oxidation). This kind of error will result in incorrect balance constraints and must be distinguished from measurement errors. As adaptive observers are used for monitoring batch fermentations, tests l_t(k), lⁱ_t(k) and l^-i_t(k) are not sensitive to this type of error. This is because the effect of these errors is manifested as a correction of some parameter values, resulting in residuals close to the zero-mean known-covariance hypothesis. Nevertheless, since the specified matrix C_{Y I/J} is incorrect, indicators h_t(k), hⁱ_t(k) and h^-i_t(k) will trigger alarms. To summarize, alarms triggered by variables h_t(k), hⁱ_t(k) and h^-i_t(k), while variables l_t(k), lⁱ_t(k) and l^-i_t(k) indicate normal process operation, can be interpreted as an incorrect system description. Other evidence that permit to detect a composition error are (i) no measurement error can cause a residual vector of the given form and (ii) measurement errors that may cause a residual vector of the given form are checked by tests l_t(k), lⁱ_t(k) and l^-i_t(k) and found to be correct.

Incorrect Component Composition

Composition of X and P (if P is a complex product) may not be exactly known or may be time-varying, resulting in incorrect or time-varying element constraints, e.g., the biomass N content can vary with time. Therefore, the detection and diagnosis of other errors will be rather difficult. Nevertheless, there may be heuristic information that can be used within expressions (2) and (17) for reducing uncertainty about the biomass composition. For example, it is known that the biomass C content is f_C» 0.48 with a relative variance of 5% and that the degree of reduction in the biomass "mol" ^* is g_X» 4.25 with a relative variance of 4% [Erikson et al., 1978]. These data together with the elemental composition of S, S_N and element balances can be used to derive linear relationships between Y*_X/N, Y*_X/S and Y*_X/CO2. With this information and measurements y_t(k), the isolation problem can be written as the maximization of the log-likelihood function (19) subject to the EKF equations that include a modified yields matrix C^*:

Missing Components

An analogous procedure is useful to identify the formation of a suspect by-product. Thus, its isolation can be written as the maximization problem given by eqs. (19) and (20), but now C* includes a column of stoichiometric yields of the suspected by-product q in the measured components:

The mol of biomass is defined by the formula CH_b1O_c1N_d1. The degree of reduction of a mol of biomass is defined as 4+b1-2c1-3d1 [Erikson et al., 1978].

No distinction can be made between an omitted component and a component composition error. However, many components are precisely defined (O₂, CO₂ and many times S), a limited number of components (X and sometimes P) lack a defined elemental composition and a limited number of substances are suspected to be omitted by-products. Thus, if probable composition errors for X and P do not cause residuals of the given form, it can be inferred that a component is omitted from constraints (2). It must be noted that, since the matrix C* defined by eq. (23) has one more column than the matrix C* defined by eq. (22), the degree of freedom is reduced to 1 and therefore the threshold q has to be corrected. Finally, if a component composition error or a missing component is detected, the states and parameters of the process can be re-estimated by maximization of the log-likelihood function (19) from the identified occurrence time to the present.

c) Detection and Isolation of Abrupt Kinetic Changes (I = 3)

Abrupt kinetic changes are caused by physiological disturbances that are manifested as parametric or structural changes in the process dynamic f(•).Therefore they do not affect the state measurement relations c(•). Thus, if this kind of disturbance is produced, indicators h_t(k), hⁱ_t(k) and h^-i_t(k) will not trigger alarms because constraints (2) are still valid. Nevertheless, abrupt changes in kinetic parameter values will be reflected as a nonwhite-noise innovation sequence that will increase the value of indicators l_t(k), lⁱ_t(k) and l^-i_t(k). Summarizing, alarms triggered by indicators l_t(k), lⁱ_t(k) and l^-i_t(k), while variables h_t(k), hⁱ_t(k), and hⁱ_t(k) indicate normal process operation, mean abrupt kinetic changes. As for the previous hypotheses, they are to be isolated by maximization of the log-likelihood function (19) for various possible kinetic models.

d) Detection of Poorly Specified Variances (I = 4)

A practical detection algorithm must be able to detect small errors and must be reliable. These two properties are in conflict with one another. In order to achieve high reliability, strong indications of error are necessary. But high sensitivity to error means that indicators may respond to minor disturbances. This, however, may also be caused by measurement noises or modeling errors instead of real faults. If specified measurement variances are too small, detection tests l_t(k), lⁱ_t(k), l^-i_t(k), h_t(k), hⁱ_t(k) and hⁱ_t(k) will be very sensitive and they will produce too many false alarms. On the other hand, large specified variances will cause excessive smoothing and detection tests will be insensitive to real faults. Thus, well-tuned filters are crucial for detection reliability and sensitivity. Incorrect measurement variances can only be detected if sufficient samples have been taken. Thus, measurement variances can be estimated on-line and compared with specified variances for performing the following tests:

are larger than the specified ones, .

Estimated variances, , are smaller than the specified ones, .

are similar to the specified ones, .

Measurement variances can be evaluated in a w’ lag sliding window by

Furthermore, effects of poorly specified variances can be evaluated by on-line computing residual means:

Values for these means significantly different from zero indicate that the observer had “lost” information due to specified variances that were too large. If specified variances are incorrect, the estimated ones can replace them and the estimation procedure can continue with the new variances.

In order to avoid frequent retuning of the filters and interference with detection tests, time windows w’ used for monitoring residual means and variances must be considerably larger than the ones used for detection purposes. Noise correlation tests are also advisable, especially in the case of on-line measurement.

COMPUTATIONAL IMPLEMENTATION AND TESTS RESULTS

The joint estimation-detection algorithm performs three tasks: (i) It propagates the states of the system and the error covariance matrix from one observation time to the next one. (ii) It conducts tests to determine whether or not a fault has occurred. This is done by generating l_t(k), lⁱ_t(k) and l^-i_t(k) indicators in a sliding window of selected size and, if an off-line sample is available, by computing h_t(k), hⁱ_t(k) and h^-i_t(k) indicators. These variables are compared with their critical thresholds as described previously and a fault condition is identified when a critical value is exceeded. (iii) It updates the state and parameter estimates and the error covariance matrix after a measurement is processed. This update is done using the EKF correction equations if no off-line measurement is available or by maximization of the log-likelihood function (9) if an off-line measurement is available. A flow chart that summarizes the estimation-detection procedure is presented in Fig. 2.

The Case Study

Shu and Yang [1991] studied the effects of the temperature on xanthan gum batch fermentations. They proposed a kinetic model in which growth is modeled by the logistic equation. The equation of Luedeking-Piret was used to model the product formation rate and the glucose consumption rate. The parameters of these equations were expressed as functions of temperature. These equations and the structured model proposed by Pons et al. [1989] were used to construct a stochastic model (Table 2) for utilization as a fermentation simulator [Dondo, 2000]. On the other hand, Cacik et al. [2001] used the model of Shu and Yang for calculating an optimal temperature profile that permits a given quantity of product to be obtained in a minimum time. The resulting profile shows an abrupt temperature change at a time that is a function of the initial biomass concentration and of the desired final product concentration. When this profile is applied to a real fermentation it may produce unexpected abnormal behaviors [Dondo, 2000]. Thus, it will be applied to a modified fermenter simulator. The modifications were introduced for simulating the following abnormal behaviors:

• Measurement biases;

• Errors in some stoichiometric yields Y

  _I/J;

• Generation of an unidentified by-product (acetate salts);

• Intracellular accumulation of metabolic products as a consequence of abrupt changes in temperature;

• Abrupt changes in some kinetic parameters as a consequence of abrupt changes in temperature.

Equations for simulating these abnormal behaviors are presented in the ^appendix. Models of the observers used by the estimation-detection procedure are also summarized in the ^appendix.

Test Results

Test cases are presented below to demonstrate the ability of the estimation-detection program to detect different faults. For each test, the system was simulated with measurement noises included and the faults were imposed at a specified time. The estimation-detection algorithm processed the measurements on-line. The aim was to assess the ability of detectors to accurately determine the presence of an abnormal behavior within a short time following the actual occurrence of the event as well as the reliability of the detector in not giving false alarms. For each observer, continuous-time nonlinear differential equations were used to predict the states of the system from one observation to the next by using a Runge-Kutta 4th-5th order algorithm [Forsythe et al., 1977]. The Jacobian matrix computed in the linearization process was used to compute the differential Ricatti equation and the covariance matrix for each time interval. By using the EKF monitor, the updated estimates and the updated covariance matrix were computed and used as initial conditions for the next interval. Afterwards, if there was an off-line measurement available, maximum likelihood estimates were obtained by using the downhill-simplex minimization method [Press et al., 1992]. The monitor filter provided initial states and parameter values and after the optimization it was reinitialized with maximum-likelihood states, parameters and covariances. The routines were written in a commercial programming language and optimizations carried out for 20 seconds in a Pentium 166 Mhz 64 MB RAM PC. The monitor filter also computes the innovation covariance matrix that is used to calculate decision variables l_t(k), lⁱ_t(k) and l^-i_t(k) on a 12 sample lag sliding window. Alarm thresholds have been fixed in the values described on (Table 3).

When an alarm was triggered, maximum likelihood suspected faults were diagnosed. Computation times for testing each hypothesis were on the order of 15 seconds and ten different hypotheses were evaluated. The reliability of the resulting diagnosis was also studied.

a) Unfaulty Fermentation

To observe the behavior of the estimation-detection methodology, an unfaulty case is presented. Initial conditions and stochastic parameters used to simulate this run are presented in Table 4. Fig. 3 shows the evolution of indicators l_t(k), l^O2_t(k) and l^CO2_t(k) and Table 5 shows the evolution of indicators h_t(k), hⁱ_t(k) and h^-i_t(k) for this case. These indicators were compared to their thresholds for detection tests and no faults were detected. Thus, estimated states are assumed to be correct. They can be compared with ‘real’ states in Figs. 4. Note that, in spite of errors in some stoichiometric yields and uncertainty in the kinetic parameters, an excellent agreement between estimated and real state variables was obtained.

b) Faulty Scenarios

Test cases are presented below to demonstrate the ability of the estimation-detection program to identify an unexpected event and its characteristic parameters. In these tests, the stochastic model was started with the parameters in Table 4 under normal operating conditions and the states were propagated before a given bias was imposed at a given occurrence time. Fault parameters are summarized in Table 6. The tests always detected the occurrence of a fault, reflecting the increased detection power arising from multiple tests as can be seen in Table 7 and Figs. 5 to 7. Fault isolation is achieved by hypothesizing alternative sources and by using the maximum likelihood identification procedure in a ten-hour time window (the six hours previous to the alarm launch plus four hours from the alarm launch for collection of dynamic information). The identification is to estimate faulty parameters and log-likelihood functions for each alternative. The capability of the technique to distinguish between alternative hypotheses is reflected by the results in Table 8. The first table compares the term S_t(j)=t0^t(j)=t(k) [(y_t(j) – c_IJ(x,u,t)) V^-1_t(j) (y_t(j) – c_IJ(x,u,t))^T] of the likelihood function (19) for each hypothesis (I,J). This is the critical term in eq. (19), since S_t(j)=t0^t(j)=t(k) ln½V_t(j)½ was almost constant. The second table compares S_t(j)=t0^t(j)=t(k) (h_t(k)) for all hypotheses. Note that this summation provides useful information for further discriminating between different hypotheses. Log-likelihood values for many hypotheses are rather similar, indicating similarity of dynamics for different error sources. Furthermore, as can be expected, measurement faults are clearly easier to distinguish than incorrect system descriptions. Nevertheless, the true hypothesis was correctly selected in most cases, even with high noise and parametric uncertainty. Estimated faulty parameters are shown in Table 9. They can be compared with true fault parameters presented in Table 6. The results show that faulty parameters can be estimated with fairly good precision.

c) Estimation of Measurement Variances

In Figures 8 and 9 estimated means and variances of measurement residuals for the run in section 4.2.1 are plotted. They were calculated in a w' = 75-sample lag-time window for on-line measurements and w' = 5-sample lag-time window for off-line measurements. Note that a relatively accurate prediction of on-line measurement variances can be obtained. For off-line measurements, the estimated variances are not very accurate. Residual means for on-line measurements are approximately zero-mean but clearly correlated and nonwhite-noise distributed. Residual means for off-line measurements are clearly biased. This is due to the effect of nonlinearities and compensation for yield errors. Thus, it would be wrong to blindly trust detection tests. Rather, one must expect sudden and large changes in the values of statistical indicators. With this relaxation, failure detectors can give very good results as was shown previously.

CONCLUSIONS

An approach to estimate state and parameters and to isolate unexpected events in batch fermentations with nonlinear and uncertain dynamics was developed. It is based on the application of several statistical detection tests and maximum likelihood state and parameter estimation techniques. The approach is designed for faulty structure discrimination. A maximum likelihood filter is used to identify faults. For computational efficiency, the fault parameters are estimated in a fixed-size sliding window. Under null hypothesis, the outputs of the algorithm are the fermentation state and parameter values. Under fault hypothesis, the outputs are states, maximum-likelihood fault parameters and log-likelihood function values. These values are used for statistical comparison with the alternative faulty hypothesis. The original contributions of the method are

• The application of multiple tests, including measurement-dedicated detection tests of the residuals of the monitor filter and balance equations.

• The on-line implementation of maximum-likelihood state and parameter estimation within the detection procedure for both the unfaulty process model and faulty models using a robust (Jacobian-free and Hessian-free) optimization method.

The technique was illustrated for simulated xanthan gum batch fermentations and ten different faulty scenarios were simulated. In spite of nonlinearities, parametric uncertainty and kinetic variations, hypothesis discriminability was very good.

Areas of continuing work include the application and development of data-fusion techniques to fuse data from tests with maximum-likelihood estimates for gaining computational efficiency and hypothesis discriminability in small sample data. Research on the use of the algorithm in multiple-faults tests and in real fermentations is also advisable.

NOMENCLATURE

j Critical threshold of statistical indicator l j^I Critical threshold of statistical indicator lⁱ j^--i Critical threshold of statistical indicator l^-i q Critical threshold of statistical indicator h q ⁱ Critical threshold of statistical indicator hⁱ q^--I Critical threshold of statistical indicator h^-i Dm Disturbance of the specific biomass growth (h^-1) Da Disturbance of the growth-associated specific metabolite production (g P/g X) Db Disturbance of the steady specific metabolite production (g P /g X h) DCO₂ Cumulative carbon dioxide production (g CO₂/l) DK^R Disturbance of the specific substrate consumption for cellular maintenance (g S/g X h) DO₂ Cumulative oxygen consumption (g O₂/l) DS Cumulative main substrate consumption (g S/l) DX_S Disturbance of the stationary biomass concentration (g/l) g Residuals from the state and parameter estimation filter e Residuals from the balance-based estimation of measurements n Residuals vector from the element-balances constraints s_I Measurement i residuals-variance m_I Measurement i residuals-mean m_M Specific biomass growth on the logistic equation (h^-1) W Time window for estimating the maximum likelihood states and parameters a Growth-associated specific metabolite production (g P /g X) b Steady specific metabolite production (g P /g X h) CPR Carbon dioxide production rate (g CO₂ /l h) f(•) State-space dynamic of the fermentation

c(•)

State-measurements relationship (nonlinear case) c_X Jacobian matrix of the vector of state-measurements relationships C State-measurement relationships (linear case) C_{Y I/J} Matrix of stoichiometric yields h Global balance-based fault indicator hⁱ Balance-based fault indicator dedicated to measurement i. h^-i Balance-based fault indicator dedicated to all but measurement i. K^R Specific substrate consumption for cells maintenance (g S/g X h) l Global residual-based fault indicator lⁱ Residual-based fault indicator dedicated to measurement i. l^-i Residual-based fault indicator dedicated to all but measurement i. L Log-likelihood function OUR Oxygen uptake rate (g O₂ /l h) p Vector of kinetic parameters P Metabolic product concentration (g P /l) P States covariance matrix P _e Covariance matrix of measurements estimated by using balances constraints Q System noise covariance matrix S Measurement-noise variance matrix S Main substrate concentration (g S/l.) S_N Nitrogen source concentration (g S/l.) S^R Cumulative amount of oxidized main substrate (g.S /lt.) u Control variable V Covariance matrix of the estimation residuals w w' Sliding-window lags X Biomass concentration (g X /l) x State-variable vector X_S Stationary biomass concentration (g X /l) y Measurement vector Y_I/J Stoichiometric yield of component I to component J (g I/g J)

Received: January 6, 2003;

Accepted: March 9, 2004

A.1. FERMENTER SIMULATOR

A stochastic model [Dondo, 2000] is used as fermenter simulator. In this model (Table 2), initial conditions, kinetic parameters and physiological parameters are stochastically generated.

A.2. EQUATIONS FOR SIMULATION OF FAULTS:

Measurement biases (I = 1)

Faulty measurements are simulated as the sum of non-unfaulty measurements y_t(k), a measurement bias D_J(u, t) and uncorrelated zero-mean Gaussian noises n_t(k).

y_{I=1, J}(u, t)_t(k) = y_t(k) + D_J(u, t) + v_t(k)

J = O₂, CO₂, S, X, P

Incorrect system description (I =2)

Errors in the stoichiometric yields Y_i/J: Faulty measurements are obtained simulating the fermentation with imposed unexpected values of some physiological parameters.

Generation of unidentified by-products: Faulty measurements are the product of an expanded matrix, which incorporates a column of yield coefficients that relate the production of the by-product q with the net production or consumption of measured component J, by an expanded state vector, which incorporates the by-product as a new state.

J = O₂, CO₂, S, X, P

Abrupt changes in some kinetic parameters (I = 3)

The system has the same structure as the nominal dynamic but with a different and deterministic initialization of parameters Xs and b after the temperature jump.

A.3. ALGORITHM MODELS

The state variables are the biomass concentration [X], the product concentration [P] and the quantity of oxidized main substrate (S^R). Equations m(u), X_S(u), a(u) and b(u) and K^R(u) represent the kinetic parameters as functions of the temperature [Shu and Yang, 1991]. Dm, DX_S, Da, Db and DK^R are disturbance parameters used to compensate for the uncertainty of these functions in respect to the complex fermentation dynamic. The structures of observers are

Monitor filter:

Maximum-likelihood filter for estimation of kinetic parameters:

Maximum-likelihood filter for estimation of fault parameters:

where

Y_X/CO2=2.96+1.19[8.0+ d(2,1)D(t,u)]

Y_X/S =-0.78+0.32[8.0+ d(2,1)D(t,u)]

Y_P/CO2=36.872-5.24[0.38+ d(2,2)D(t,u)]

Y_P/S =0.97875-0.1625[0.38+ d(2,2)D(t,u)]

Y_P/O2=20.9124-13.312[0.38+ d(2,2)D(t,u)]

In this model, d(I,J) is a binary variable denoting that hypothesis (I,J) is tracked:

while D(t,u) is a linear function of the fault occurrence time and magnitude:

A.4. INITIALIZATION OF ESTIMATORS

Monitor filter:

Trace P₀^T = [0.020, 0.005, 0.005, 0.050, 0.100, 0.200, 0.050, 0.020]²

Trace Q ^T = [0.000, 0.000, 0.000, 0.005, 0.020, 0.015, 0.003, 0.003] ²

Trace S ^T = [0.05, 0.05, 0.10, 0.15, 0.40]²

Maximum-likelihood filter for estimation of kinetic parameters:

P₀ = P_t(j)=t0 (From the monitor filter)

Trace Q ^T = [0.001, 0.001, 0.001] ²

Trace S ^T = [0.05, 0.05, 0.10, 0.15, 0.40]²

Maximum likelihood filter for estimation of faulty parameters:

P₀ = P_t(j)=t0 (From the monitor filter)

Trace Q ^T = [0.003, 0.003, 0.003] ²

Trace S ^T = [0.05, 0.05, 0.10, 0.15, 0.40]²

Other parameters:

On-line sampling and estimation frequency: 25 h^-1

Average off-line sampling frequency (monitor mode): 0.33 h^–1

Average off-line sampling frequency (diagnostic mode): 1.00 h^–1

Average off-line measurements delay: 0.25 h

Sliding window lag for estimation of kinetic parameters (W)6.00 h

• Alcorta García, E. and Frank, P. M., Deterministic Nonlinear Observer-based Approaches to Fault Diagnosis: A Survey, Control Engineering Practice, Vol. 5, 5, 663-670 (1997).
• Bastin, G. and Dochain, E., On Line Estimation and Adaptive Control of Bioreactors, Elsevier Amsterdam (1990).
• Bastin, G. and Dochain, E., On-line Estimation of Specific Growth Rates, Automatica, Vol. 22, No. 6, 705-709 (1986).
• Cacik, F., Dondo, R. and Marqués, D., Optimal Control of a Batch Bioreactor for the Production of Xanthan Gum, Computers & Chemical Engineering Journal, Vol. 25, 409-418 (2001).
• Dondo, R. and Marqués, D., Mass and Energy Balances as State-space Models for Aerobic Batch Fermentations, Latin American Applied Research, Vol. 32, 195-204 (2002).
• Dondo, R., A Method for Detection and Diagnosis on Batch Fermentations, ISA Transactions, 42, 135-147 (2003).
• Dondo, R., Estimación y Control Óptimo en Biorreactores Batch, Doctoral Thesis, Universidad Nacional del Litoral - Argentina (2000).
• Erikson, L., Minkevich, I. and Eroshin, V., Application of Mass and Energy Balance Regularities in Fermentation. Biotechnol. Bioeng., Vol. 20, 1595-1621 (1978).
• Forsythe, G., Malcolm, M. and Moler, C., Computer Methods for Mathematical Computations, Englewood Cliffs, N.J., Prentice Hall (1977).
• Frank, P., Fault Diagnosis in Dynamic Systems Using Analytical and Knowledge Based Redundancy A Survey and Some New Results, Automatica, Vol. 26, No. 3, 459-470 (1990).
• Gudi, R.D., Shah, S.L. and Gray, M.R., Multirate State and Parameter Estimation in Antibiotic Fermentation with Delayed Measurements, Biotechnol. Bioeng., Vol. 44, 1271-1278 (1994).
• Halme, A., Measurement and Estimation in Bioreactors, Proceedings One Day International Workshop on Control of Biotechnical Processes. Newcastle, U.K. (1987).
• Isermann, R., Process Fault Detection Based on Modeling and Estimation Methods A Survey, Automatica, Vol. 20, No. 4, 387-404 (1984).
• Mendel, J., Lessons in Estimation Theory for Signal Processing, Communications and Control, Prentice Hall (1995).
• Minkevich, I., Mass-Energy Balances for Microbial Synthesis, Biochemical and Cultural Aspects, Biotechnol. Bioeng., Vol. 25, 1267-1293 (1983).
• Pons, A., Dussap C. and Bros J., Modelling Xanthomonas Campestris Batch Fermentation in a Bubble Column, Biotechnol. Bioeng. Vol. 33, 394-405 (1989).
• Press, W., Teukoslky, S., Veterling, W. and Flannery B., Numerical Recipes. The Art of Scientific Computation, 2^nd ed., Cambridge University Press (1992).
• Shu, C. and Yang, S., Kinetics and Modeling of Temperature Effects on Batch Xanthan Gum Fermentation, Biotechnol. Bioeng. Vol. 37, 567-574 (1991).
• Stephanopoulos, G. and San, K.Y., Studies on On-line Bioreactor Identification. I Theory, Biotechnol. Bioeng., Vol. 26, 1176-1188 (1984a).
• Stephanopoulos, G. and San, K.Y., Studies on On-line Bioreactor Identification. II Numerical and Experimental Results, Biotechnol. Bioeng., Vol. 26, 1189-1197 (1984b).
• Van der Heijden, R. T., Romein, B., Heijnen, J.J., Hellinga C. and Luyben, K. H., Linear Constraint Relations in Biochemical Reaction Systems: II Diagnosis and Estimation of Gross Errors, Biotechnol. Bioeng., Vol. 43, 1, 11-20 (1994b).
• Van der Heijden, R.T., Heijnen, J.J., Hellinga, C., Romein, B. and Luyben, K. H., Linear Constraint Relations in Biochemical Reaction Systems: I Classification of the Caluculability and the Balanceability of Conversion Rates, Biotechnol. Bioeng., Vol. 43, 1, 3-10 (1994a).
• Wang, N. and Stephanopoulos, G., Application of Macroscopic Balances to the Identification of Gross Measurement Errors, Biotechnol. Bioeng., Vol. 25, 2177-2208 (1983).
• Wilsky, A., Detection of Abrupt Changes in Dynamic Systems, Lecture Notes in Control and Information Sciences, 27-49, Springer-Verlag (1986).
• Young, P., Parameter Estimation for Continuous Time Models A Survey, Automatica, Vol. 17, No. 1, 23-39 (1981).

appendix

Publication Dates

• Publication in this collection
  14 July 2004
• Date of issue
  Sept 2004

History

• Received
  06 Jan 2003
• Accepted
  09 Mar 2004