Abstract
A state estimator and various model-based control systems have been designed for a real anaerobic digestion (AD) pilot reactor fed with dairy manure. The model used is a modified Hill model which is a relatively simple dynamical AD process model. The state estimator is an Unscented Kalman Filter (UKF) which uses only methane gas flow measurement to update its states. The model and the state estimates are used in different control systems. One of the control systems aims at controlling the methane gas flow to a setpoint. Simulations indicate that the setpoint tracking performance of a predictive control system is considerably better comparing with PI control, while disturbance compensation is not much better. Consequently, assuming the setpoint is constant, the PI controller competes well with the predictive controller. A successful application of predictive control of the real reactor is presented. Also, three different control systems aiming at retaining the reactor at an operating point where the volatile fatty acids (VFA) concentration has a maximum, safe value are designed. A simulation study indicates that the best control solution among the three alternatives is PI control based on feedback from estimated VFA.
1. Introduction
Anaerobic digestion (AD) of organic substrates can produce biogas which consists mainly of methane and carbon dioxide, [1–3]. In a well-operated AD reactor, the methane content is sufficiently large to make the biogas combustible; that is, the AD process produces applicable energy. Moreover, the reactor digestate is often high in nutrients and can be used in fertilization. Animal waste, in many cases combined with, for example, food waste, is a typical feedstock of AD reactors. A presentation of AD of animal wastes, from dairy, beef, poultry, and swine, is provided, for example, in [4].
UASB (upflow anaerobic sludge blanket) type reactors are effective AD reactors as they allow for relatively high load rates (feed rates) and/or small reactor volumes, [1, 5]. The effectiveness is due to relatively large solids retention time (SRT), which is the retention time of the microorganisms which degrades the substrate and generates, for example, methane, compared with the hydraulic retention time (HRT) of the reactor. The AD reactor studied in the present paper is a UASB reactor.
Anaerobic digestion is a complex and nonlinear dynamic process and most plants suffer from a lack of robust online-measurement systems for online process monitoring [3]. Therefore, automatic plant control is a challenging task. The present paper presents an attempt to use a mathematical dynamic model to estimate online, nonmeasured AD state variables and to use these estimates in a model-based control system. Results of the application of state estimation and model-based control to a real pilot AD reactor using dairy waste as feedstock are shown. The reactor is situated at Foss Farm, Skien, Norway. The results from the pilot reactor are assumed to be transferable to a planned full-scale reactor at the farm.
In this paper, state estimates are used both in industry-standard PI controllers and in predictive controllers. The only online measurement used by the estimator, and thus by the controllers, is the methane gas flow. The reactor temperature is retained at a constant setpoint by means of a temperature control system [6].
Several control systems are designed and applied to the reactor. One aims at retaining the produced methane flow at a setpoint which can stem from a specified power production. Another control system aims at retaining the reactor at a safe operating point, where the concentration of the VFA (volatile fatty acids) is not above a certain value.
The model-based design and the simulations are based on the modified Hill model adapted to the pilot reactor [7]. This model is summarized in Section 3.3.
This paper is organized as follows. Section 2 gives a literature review. Section 3 provides a system description, including the mathematical reactor model used as the basis for state estimation and model-based control. In Section 4, safe reactor operation conditions are defined in terms of an acceptable range of VFA (volatile fatty acids). Section 3.2 presents a general structure of a model-based optimization and control system, applicable to the reactor. Application of the Unscented Kalman Filter (UKF) to estimate the state variables of the reactor and its main disturbance, namely, , is described in Section 5. These estimates are used for control of , which is described in Section 6, which includes both simulated and real results. The estimates are also used for control of , which is described in Section 7, which is simulation study. Conclusions are given in Section 8.
Matlab and Simulink (MathWorks, Inc.) are used for numerical computations and simulations. The real control system is implemented in LabVIEW (National Instruments, Inc.) running on a laptop PC. In the LabVIEW program, the algorithms of the UKF and the predictive controller are implemented in a Matlab Script Node.
2. Literature Review
2.1. State Estimators for AD Reactors
Literature about state estimators applied to AD reactors fed specifically with dairy manure has not been found. Below are references to state estimators applied to reactors fed with other types of substrates, assumed to be also relevant for the present application.
In a simulation study, Jones et al. [8] apply an Extended Kalman Filter (EKF) to estimate four states of a simplified version of the AD model by Hill and Barth [9], using five online measurements.
Bernard et al. [10] estimate the six states of a real AD reactor fed with effluents from a wood processing plant using an asymptotic observer [11]. Available online measurements were CH4 gas flow and CO2 gas flow. Influent concentrations are assumed to be known. The estimator is based on a state variable transformation leading to a model having auxiliary state variables where the reaction rates are eliminated. These rates are then estimated from the state estimates. The estimator is designed so that the estimation errors converge towards zero with dynamics of the mass balances of the model, determined by, for example, the feed rate. The asymptotic observer is an open-loop estimator and has no tuning parameters, contrary to a Luenberger observer and a Kalman Filter which are closed loop, or feedback estimators with parameters which can readily be used for performance adjustment.
Alcaraz-González et al. [12] estimate four out of six states of a real AD reactor fed with industrial wine distillery vinasses, namely, the methanogens and acidogens concentrations, COD (chemical oxygen demand), and alkalinity, by using online measurements of CO2 gas flow, VFA, and TIC (total inorganic carbon). The AD process model is as in [10]. The estimator is an interval observer based on the structure of an asymptotic observer. An important property of an interval observer is that the estimates are guaranteed to be within bounds given by uncertainty bounds of model parameters and AD process inputs.
In a study based on real data, Theilliol et al. [13] estimate the six state variables and three unknown inflow concentrations, namely, COD, VFA, and TIC, of an AD reactor fed with industrial wine distillery vinasses, using five online measurements: COD, VFA, alkalinity, CH4 gas flow, and CO2 gas flow. The estimator is based on manipulating the original state space model using SVD (singular value decomposition) to find an observable subsystem insensitive to unmeasured inputs. Then, a Luenberger observer based on this subsystem is used to estimate the state and the unmeasured inputs.
In a simulation study based on a full-scale agricultural biogas plant, Gaida et al. [14] use discriminant analysis and classification-based pattern recognition methods to find the static mapping function between the measurement data, which are biogas flow, CH4 and CO2 gas concentrations, pH in the reactor, the amount of each substrate, and the state of the AD process. The state variables are those of the ADM1 model (Anaerobic Digestion Model No. 1) [15]. The various substrates considered are maize silage, grass, manure, and manure solids.
Dochain [16] and Bogaerts and Vande Wouwer [17] give an overview of various state estimators suitable for bioprocesses, including the estimators applied in the references above.
In the applications referred to above, the estimators use two or more online measurements. In the present paper, only one measurement is used, namely, (CH4 gas flow). Furthermore, in the present paper the Unscented Kalman Filter (UKF) is used. The UKF can be used without any linearization or model manipulation; that is, it uses the nonlinear state space model directly in the algorithm. We have not found literature on application of the UKF to AD reactors.
2.2. Model-Based Control of AD Reactors
We have not found literature on model-based control systems of AD reactors fed specifically with dairy waste. Below are references to model-based control systems of reactors fed with other types of substrates, assumed to be also relevant for the present application.
Bernard et al. [10] have implemented a model-based adaptive linearizing controller and a fuzzy controller designed to maintain the intermediate alkalinity (VFA, volatile fatty acids) and the total alkalinity within specified limits to ensure stable process conditions and to avoid VFA accumulation despite organic load disturbances. The so-called AM2 model, [18], is used for design and simulation.
Puñal et al. [19] have designed an automatic fuzzy logic-based control system to maintain the online-measured VFA concentration at a proper setpoint.
Méndez-Acosta et al. [20] have designed a model-based controller for maintaining the COD (chemical oxygen demand) of the reactor effluent at its setpoint, using the AM2 model, [18].
Méndez-Acosta et al. [21] have designed a multivariable control system for controlling the concentration of VFA in the reactor to its setpoint using the feed rate and controlling the total alkalinity to its setpoint using the addition of an alkali solution.
Strömberg et al. [22] have identified, using simulations, three controllers for AD processes to be the most suitable ones for maximizing gas production, while being able to react properly to process disturbances due to variations in pH, ammonia, and concentration in the reactor feed. The simulations use the ADM1 model [15]. All of the controllers have the feed rate as control variable (controller output). The controllers resemble an expert system, with logics (if-clauses) in the control function. The three controllers are the extremum-seeking variable gain controller by Liu et al. [23], the disturbance monitoring controller by Steyer et al. [24], and the hydrogen-based variable gain controller by Rodríguez et al. [25]. Strömberg et al. [22] note that no uniform tuning method could be derived to tune the three controllers. Instead, trial-and-error procedures are used.
In a simulation study, Gaida et al. [26] have implemented a nonlinear predictive controller to control a simulated ADM1, assuming all states are available, and, therefore, a state estimator is not used. The controller allows alternative optimization criteria, for example, economical optimization and minimum methane concentration of the biogas. The plant is the same as in [14], compared with the above section about state estimation.
In a simulation study, Ordace et al. [27] have implemented a predictive controller based on transfer functions adapted to the ADM1 model to control the ADM1. The optimization criterion of the controller contains the square of the control error, while the control signal usage is not included; that is, it has no cost in the criterion.
3. System Description
3.1. AD Reactor with Control System
Figure 1 depicts the AD reactor with its control system. The reactor type is UASB (upflow anaerobic sludge blanket). The reactor is part of a pilot biological plant for nutrient and energy recovery named Foss Biolab, situated at Foss Farm, Skien, Norway. Input to the plant is dairy manure diluted with 25% water and filtered with a sieve, and outputs are fertilizer and biogas consisting of approximately 70% methane. The reactor’s temperature is kept fixed at its setpoint with an automatic temperature control system.
In Figure 1, the block denoted ‘‘Model-based controller” may comprise a state estimator and alternative controller functions (predictive controller and PI controller with feedback from state estimates). The model-based controller uses an online measurement of which is provided by sensor FT. This measurement is obtained by multiplying the online biogas flow measurement from a thermal gas flow sensor and the online methane concentration measurement from an IR-based sensor. The raw measurement signals are smoothed using software filters.
is used as control variable. The demanded flow is obtained with a peristaltic feed pump operated with PWM (Pulse width modulation) with a cycle time of 700 sec.
In principle, is also a candidate as control variable since it has a clear impact on , but in [28] we argue why is not considered a usable control variable.
An online measurement of is used by the controller, since is an important model variable. is retained at its (fixed) setpoint with a separate temperature control system, where the controller is a PI (proportional plus integral) controller [6].
In this paper, is kept at 35°C because this is a typical temperature at which AD reactors are operated (mesophilic conditions). However, this temperature is not necessarily optimal. In [29] we show how the temperature can be specified using model-based optimization.
3.2. Control System Structure
Figure 2 shows the structure of the control system.
In the block diagram: , and . comprises here the four state variables of the modified Hill model, compared with Section 3.3: . Depending on the applications in this paper, , compared with Section 6, or , compared with Section 7. Furthermore, the Process is the reactor. The Controller implements predictive control, PI control, or manual control. The Estimator is an Unscented Kalman Filter (UKF). The Control Designer is the algorithm or strategy used to transform the specifications of the optimal operation into (optimal) setpoints and/or control signals. The Control Designer may also set parameters for controller tuning, for example, cost factors in the optimization criterion of a predictive controller, or it may be an optimization algorithm to calculate optimal setpoints.
The symbol in various blocks in Figure 2 represents the assumed mathematical model used in the block. The symbol in the Process block is the model representing the real system (process). Only if model errors are assumed to be zero, and will be identical.
The connections from and/or to the Control Designer are due to being an input to the process, and the value of or is included in the model-based optimization. For example, the value of in the feed of the reactor has an impact on the specific value of needed to produce a specified which in turn is closely related to the power production in the reactor.
In general, the operational objectives, which are the inputs to the Control Designer in Figure 2, may be adjusted based on results of an evaluation of the factual process operation, but this possible adjustment is not depicted in Figure 2.
A large number of model-based controllers exist [30]. In this paper, a predictive controller [31, 32] is selected (a predictive controller is also denoted as model-based predictive controller (MPC)). The selection of a predictive controller is due to its popularity (as model-based controller) in the process industry [33] and due to our view that it implements most of the important controller features which would otherwise require a number of special solutions, that is, feedback, feed forward, integrator antiwindup, constraints handling, and time-delay compensation. When nonlinear predictive control is used, as in this paper, process nonlinearities are taken into account naturally and without approximations. Furthermore, a predictive controller is relatively easy to tune, if the process model is accurate.
3.3. AD Process Model
The mathematical model of the AD processes in the reactor is a modification of the Hill model [34] adapted to the pilot reactor [7]. The model is based on material balances of biodegradable volatile solids, volatile fatty acids, acidogens and methanogens, and a calculation of the produced methane gas flow. The model is summarized below:
material balances:
methane gas production:
reaction rates:
Table 1 shows model parameter values as adapted to AD reactor at Foss Farm, [7].
One example of a set of steady-state values of the AD process variables is given in Table 2.
4. Safe Operation Condition
The various control systems proposed in this paper are designed to retain the reactor at a safe reactor operation condition, defined below. Hill et al. [35] have found, from a comprehensive study of literature reporting operational data for reactors fed with swine and beef manure and confirmed by their own laboratory experiments, that g/L indicates an impending reactor failure, causing a reduction of methane production. Hence, it is here stated that defines safe operation conditions for the reactor. For practical reasons, we have not been able to conduct our own experiments to verify inequality (4) or to identify a different . However, a new value of will not change the principal results of this paper.
Hill et al. [35] found that also the propionic to acetic acid (P/A) ratio is a good indicator of health. However, this ratio cannot be calculated from the mathematical model used in this paper, and, therefore, the analysis here is not based on this ratio.
Hill et al. [35] did not use dairy manure in their analysis since reliable data for such manure were not available. Nevertheless, it is here assumed that the aforementioned safe range of also applies approximately for reactors fed dairy manure. A support for this assumption is that the validated AD reactor model by Hill [34] has the same parameters describing the AD process for dairy, swine, poultry, and beef manure, except for parameters expressing the fraction of the organic feed that is degradable, but the AD process dynamics are independent of the latter parameters.
Figure 3 shows simulated static (steady-state) responses in a number of variables to a range of constant feed rates (). The cyan horizontal line in the plot represents g/L. The green intervals on the abscissas indicate safe reactor operation, and, conversely, the red interval indicates unsafe operation.
(a)
(b)
(c)
(d)
(e)
(f)
Table 2 shows the values of several variables at the ultimate safe steady-state operating point. The set of three corresponding values (, , and ) constitutes the ultimate safe steady state operating point of the reactor. Table 2 also shows, for completeness, values of other model parameters and variables than those discussed here.
One question arises about the applicability of the modified Hill model to predict safe/unsafe operation of the reactor. Is it necessary to include g/L explicitly to find the ultimate (maximum) safe operating point? Assuming the reactor model is accurate, safe operating points should be implicit in the model; that is, they can be calculated from the model, for example, by simulations. The modified Hill model used in the present paper is relatively simple. It is not clear to what extent the model is able to predict unsafe operation of the real reactor due 10 to high concentration of VFA. Therefore, as long as this simple model is chosen, it will be safer to define explicitly instead of relying on the model alone to predict a possible failure.
Defining explicit limits on model variables for safe operation is consistent with the approaches in, for example, [10, 21], where limits on VFA and TA (total alkalinity) are set explicitly.
5. State Estimation
State estimation is used in the control systems described in Sections 6 and 7. State estimators can also be useful solely for monitoring purposes, that is, for estimation of state variables in the lack of sensors. The state estimator used in the present paper is a Kalman Filter [36] algorithm based on the modified Hill model presented in Section 3.3. While there exist several state estimation algorithms (cf. Section 2), we select here the Kalman Filter because it has a relatively simple and straightforward structure and because it can be easily fine-tuned.
The modified Hill model is a nonlinear model. The Extended Kalman Filter (EKF) is a commonly used extension of the basic Kalman Filter for nonlinear models. The EKF involves linearization of the process model. An alternative to the EKF is the Unscented Kalman Filter (UKF) [36]. Two benefits of the UKF, compared to the EKF, are that no linearization is necessary and that the estimates are more accurate as the propagation of the estimation covariances, needed to calculate the optimal state estimates, are calculated more accurately. Because of these two benefits, the UKF is selected as state estimator in this paper.
5.1. Variables and Parameters of the Model
The state variables of the modified Hill model are (cf. Section 3.3) , , , and . They are estimated with the UKF. It is decided to also estimate with the UKF since it is assumed that its value may vary, though slowly. As is common, is modeled as a ‘‘random walk”: , where is a random disturbance. Thus, the augmented state vector to be estimated by the UKF is
is regarded as an input variable to the UKF. is the control variable, which is always known.
The model parameters are known from model adaptation [7]. may vary but is always known as it is measured continuously.
The process measurement, , used by the UKF is available from sensor FT in Figure 1. Hence, in the UKF.
5.2. Observability
The linearized reactor model, augmented with , is found observable at a number of typical operating points using the obsv function of the Matlab Control System Toolbox (further details are not shown here).
5.3. Tuning of the UKF
The tuning parameters of the UKF are as follows: (initial estimated state; the initial a posteriori estimate), (initial state estimation error covariance), (measurement noise covariance), and (process noise covariance). Ideally, these parameters are set equal to their known values, but some of them may not be available. Good tuning guidelines are actually hard to find. Even an otherwise thorough book as [36] gives little advice. In this paper the tuning is done as follows.(i) is set equal to the values from laboratory analysis at the start of the pertinent time interval. This applies ideally to , , and . However, for , we impose, for the purpose of demonstration, a large initial estimation error, by setting the initial estimate of equal to 20% of the value known from laboratory analysis. and are not known, but their initial values are calculated from the model assuming steady state (details of the calculation can be found in [7]).(ii) is set as a diagonal matrix as follows: with .(iii) is a diagonal matrix, which, since the number of measurements is one (), is reduced to a scalar—the measurement variance. From a representative real time series, (iv) is typically set as a constant matrix (diagonal). Assuming that , , and are set, can be used as final tuning parameter.(a)Increasing makes the estimate for state variable converge faster to the assumed true value, but with the drawback that the estimate for becomes more noisy (caused by the increased propagation of the measurement noise, via the Kalman Filter gain(s)).(b)Reducing has the opposite effects. It is proposed to relate the diagonal element (i.e., the process noise variance) to the magnitude of the pertinent state variable: With the initial setting of , it is found that is a proper value. Then the ultimate tuning is made by adjusting . By trial-and-error, .
5.4. Results and Discussion
Figure 4 shows estimates with the UKF together with real data from online sensors and laboratory analysis over a time interval of 85 days. (This time interval includes the interval where the UKF is applied to the real reactor as part of the predictive controller, cf. Section 6.5.) The process measurement used by the UKF is .
(a)
(b)
(c)
(d)
(e)
(f)
Overall, the UKF gives reasonably good estimates (real values of and are not known).
The large initial estimation error imposed on purpose is effectively reduced during approximately 15 days.
From d, there is a noticeable difference between the estimate and the laboratory analysis of . It is not clear what the cause of this difference is. If the model is trusted, the difference may indicate an inaccuracy of the laboratory analysis.
6. Control of Methane Gas Production
6.1. The Effect of Feedback Control
To demonstrate the effect of feedback (or automatic or closed-loop) control of , Figure 5 shows, for the real pilot reactor, experimental time-series of and (and ) with feedback control and without control. It is clear that varies less with control than without control. remains close to even after the setpoint is changed. The variations are due to inevitable disturbances. In the case of feedback control, is of course varying, while it is constant in the case of no control (i.e., open-loop control). is actually different in the two cases, but it is assumed that the difference between the two cases is independent of the temperature difference.
(a)
(b)
(c)
(d)
Whether the variation in in open-loop control is acceptable or not must be decided in each specific application. A comparison of the performance of closed-loop control and open-loop control when disturbances are assumed can be made using simulations with the AD model presented in Section 3.3.
6.2. Operational Objective and Control Strategy
It is here assumed that a sufficient rationale for feedback control of exists. The operational objective is stated as producing a demanded methane gas flow. A specific value of is related to the power, (kW), as the energy content of methane gas is 9.95 kWh/m3 at NTP.
The methane gas flow setpoint must be feasible. The feasibility can be checked with steady-state simulations. More specifically, it can be checked using the upper-left plot in Figure 3.
Furthermore, safe reactor operation must be ensured, which here means that inequality (4) is satisfied.
Relating to Figure 2, the above specifications concerning , the limitation of variations of , and the condition inequality (4) are inputs to the Control Designer. Outputs from the Control Designer are and . The latter is the cost factor of the control signal variations of a predictive controller.
6.3. Control Functions
In control system design, the PI(D) controller should normally be taken into account when different controllers are evaluated. If oscillations can be tolerated, even the on-off controller should be considered. Using on-off controllers and PI controllers for control of the pilot reactor is discussed in detail in [28].
In many cases, advanced controllers can give improved control compared with the simple PI(D) controller and the on-off controller, but typically the implementation is considerably more demanding. As argued in Section 3.2, a predictive controller is used as advanced controller in this paper. A predictive controller to retain at its setpoint is implemented both on a simulator of the reactor and on the real reactor. The model is the modified Hill model (cf. Section 3.3). A time-delay of d is included at the control input of the model: where is the feed rate of the modified Hill model and is the control signal. This time-delay accounts approximately for the dynamics not included in the modified Hill model presented in Section 3.3. The optimization objective of the predictive controller is where with constraint which is included in the optimization problem formulation; that is, it is an input argument in the fmincon function call in Matlab. is the present time instance. is the control error, . The time derivative, , represents the control signal changes. The larger the , the smoother the control actions.
In implementations, the discretized version of is minimized, giving an optimal control sequence, , over the prediction horizon. The first element of this sequence, that is, , is applied as control signal at the present time point. The prediction horizon is receding and the procedure of obtaining and is repeated as time evolves.
The prediction made by the controller is based on the modified Hill model discretized with the Euler explicit (forward) method. is calculated with the nonlinear optimization function fmincon in the Optimization toolbox of Matlab. The present state, , needed for the prediction, is calculated with the augmented Unscented Kalman Filter presented in Section 5.
6.4. Simulations
6.4.1. Controller Settings
The settings of the predictive controller in the simulations are as follows.
A time-step of d is used in the discrete-time version of the modified Hill model used for prediction. This is also the time-step of the discretization of . corresponds to time-steps, which is then the prediction horizon in number of time-steps.
in (11) is found by trial-and-error on a simulator. A proper value of the prediction horizon is found as d (with d, a change in performance can be observed).
In the simulations, the predictive controller is compared with the PI controller. The PI controller is tuned at the operating point shown in Table 2 using the Skogestad method [37], with the modification of the setting as proposed in [38]. The PI settings are [(L CH4/d)/(L feed/d)] and d.
6.4.2. Performance and Robustness Measures
The control system performance and robustness measures applied in the simulations are described in the following.
(1) IAE (Performance). The IAE index (Integral of Absolute Error) is a commonly used measure of control system performance. measures the setpoint tracking:
The measures the disturbance compensation:
(2) Control Signal Variations (Performance). As measures of the variation of the control signal, both the standard deviation, , and the mean of the absolute value of the rate of changes, , are calculated.
(3) Stability Margins (Robustness). The traditional measures for robustness of linear control systems are the gain margin (GM) and the phase margin (PM). The predictive controller is a nonlinear controller, and the (reactor) is a nonlinear process. Thus, the predictive control system and the PI control system are nonlinear systems. We propose here to expand the use of GM and PM as stability margins also for these nonlinear systems, as explained in the following.
An adjustable gain, , is inserted into the loop (between the controller and the process); see Figure 6. Normally, . The (ultimate) value that brings the (simulated) control system to the stability limit, with sustained oscillations, is found by trials. Then,
To calculate the PM, an adjustable time-delay, , is inserted into the loop; see Figure 6. Normally, . The value that brings the control system to the stability limit, that is, causing a sustained oscillation, is found experimentally on the simulator. Denote the period of the oscillation as [s]. As shown in [39] (Appendix 1), Seborg et al. [40] propose the following ranges for appropriate values of the stability margins: dB dB and .
Relating to Figure 2, and are included before the Process block, after the branch from to the Estimator.
6.4.3. Simulations
Figure 7 shows simulated time-series with predictive control and, for comparison, PI control. The initial operating point of the reactor is as shown in Table 2, which is the ultimate (maximum) safe steady-state operating point. The setpoint is varied as a sequence of two ramps of slope ±2 (L CH4/d)/d each lasting for 1 d. The disturbance is varied as a ramp of slope 2 (g/L)/d during 1 d which is a realistic variation [28].
(a)
(b)
(c)
(d)
The simulations shown in Figure 7 are without measurement noise. To measure the control signal variations, simulations have been run with measurement noise in the form of a normally distributed random signal with zero mean and standard deviations L CH4/d, which is realistic for the present reactor. The simulations are run over 10 d with a constant setpoint and a constant disturbance (simulations are not shown here).
6.4.4. Results and Discussion
Table 3 shows performance and robustness measures with predictive control and with PI control. The IAE indexes, (12) and (13), are calculated with d, d, d, and d.
Comments on the results shown in Table 3 are the following.(i) with predictive control is 13% of with PI control. Hence, predictive control is clearly the best.(ii) with predictive control is 66% of with PI control. Again, predictive control is the best, but the improvement compared with PI control is not large.(iii) with predictive control is approximately 67% of the value with PI control, while with predictive control is approximately 28% of the value with PI control. These numbers vary with the realization of the random processes generated in the simulation, but they are representative. By detuning the PI controller for more relaxed control (reducing and increasing according to Skogestad’s formulas), both and are reduced. By a proper retuning, either of them can become approximately equal to the value with predictive control. The consequence of such a retuning is that the IAE measures with PI control will increase. In one simulated example, the PI controller was retuned so that with predictive control and PI control was approximately equal. The with PI control then increased 4.5 times; that is, the control performance became radically worse. The smoother control action with predictive control compared with PI control has been observed from experiments on the real reactor.(iv)GM is acceptable with PI control. With predictive control the notion of GM is questionable, since the simulated control system does not actually become unstable for any gain increase at the process input. Rather, the gain increase is seen by the UKF as a change in the disturbance, or, more specifically, as an increase in . Consequently, the estimate of is increased, which in turn is used in the prediction by the predictive controller, causing a large overshoot or undershoot in before it eventually reaches (plots of simulations not shown). From simulations it is found that is back at its setpoint during 1-2 d for .(v)PM is larger with predictive control (63.9°) compared with PI control (47.6°).
6.4.5. Concluding Remarks
Above, the predictive controller has been compared with the PI controller tuned using a standard method, namely, the Skogestad method [37]. Simulations indicate that predictive control has better performance and better robustness than the PI controller. It can also be claimed that the predictive controller, here including the state estimator, is more intuitive to adjust since its parameters have a direct relation to practical factors such as measurement noise and control signal variation. The drawbacks with predictive control are that a mathematical model of the reactor is required and that it is more complicated to implement.
The setpoint tracking performance of the predictive controller is considerably better than that of the PI controller, while the improvement in disturbance compensation is not large. Taking into account that the PI controller is much easier to implement, it may be claimed that the PI controller is the preferred controller if the setpoint is constant.
6.5. Experiments on the Real Reactor
Predictive control has been applied to the real reactor. Some of the settings in the practical experiment differ from those used in the simulation study presented in Section 6.4, which has been accomplished approximately one year after the practical experiment. (However, simulations were used to test the control system before the practical implementation.) The differences in settings are shown below.(i)In the practical experiments, d and d. In the simulations in Section 6.4, d and d. d has been tested in simulations, giving a slight change in performance, probably due to less accurate numerical integration (explicit Euler is used). is the same both in the practical experiments and in the simulations.(ii) is limited to 40 L/d, which is also used in the simulations in Section 6.4. This limit is reached in the practical experiment but is not reached in the simulations since the perturbations are relatively small there.(iii)No time-delay term is included in the model used by the predictive controller in the practical experiment, while it is found appropriate to include a time-delay in the simulation study as the model analysis in [7] indicates that a time-delay is present.(iv)The cost factor in (11) was set to 0.8 in the practical experiment, while 0.01 was found appropriate in the simulation study (cf. Section 6.4). The smaller in the simulations may be due to dynamic phenomena of the real reactor not encapsulated by the model. In any case, is typically a tuning parameter.
6.5.1. Results and Discussion
Figure 8 shows the time-series of the practical experiments. Below are comments on this figure.(i)At d, was reduced instantly from 190 to 150 L/d. Since the reduction was instant, the predictive controller could not take any control action in advance. The response in the gas flow is stable and shows acceptable stability, but the stability is reduced compared with the simulated response. The control error is less than 3 L/d after approximately 1 d. A possible explanation of the damping of the real response being less than in the simulated response is that the predictive controller does not include any process model time-delay while, as pointed out above, there is actually a time-delay in the real process.(ii)At d, a preset-ramped setpoint profile started. The predictive controller adjusts before starts increasing. The tracking is accurate. The upper bound of of 40 L/d is eventually reached, causing the rate of change of to become less than the rate of change of .(iii)At d, the rate of change of is instantly adjusted from (L CH4/d)/d to (L CH4/d)/d. The observed lag in can be explained with the instant change of which prevents predictive control action.(iv)At d, a preset step change of from 150 to 155 L CH4/d is applied. The predictive adjustment of is obvious. shows a clear overshoot, but it is expected that the response will stabilize.(v)At d, the predictive control experiment had to be stopped as other experiments were scheduled to start at this point of time. The controller was actually set to manual mode. The saved future control signal sequence generated by the predictive control shows a declining behavior, indicating that eventually would have been brought back to its setpoint.
(a)
(b)
As pointed out earlier, the methane gas flow setpoint must be feasible. For the above experiments, the feasibility can be checked using the upper-left plot in Figure 3. According to this plot, the setpoint values used in the experiments (cf. Figure 8) are actually feasible.
7. Control for Safe Reactor Operation
7.1. Objective and Control Strategies
Here, the operational objective of the reactor is defined as retaining the reactor at the ultimate safe steady-state operating point given in Table 2 (this is the input to the Control Designer in Figure 2). To this end, the following three alternative control strategies are tested (they comprise the “output” from the Control Designer in Figure 2).(1) is controlled to a setpoint of , which is L/d, assuming the operating point shown in Table 2. This control strategy is described in Section 7.2.(2) is controlled to a setpoint of , which is g/L according to Table 2. In principle, this control requires feedback from the measurement of . Such sensors do exist [41, 42], but they are not in use on the present reactor. Instead, the estimate of calculated continuously with a state estimator (Kalman Filter) is used (cf. Section 5). This control strategy is described in Section 7.3.(3) is controlled to a setpoint of , which is L CH4/d according to Table 2. This control requires feedback from the measurement of . This control strategy is described in Section 7.4, where also PI control is applied for comparison.
In each of the control strategies, the feed rate is used as control variable, (cf. Section 3.2).
The applicability of the three control strategies described above is demonstrated with simulations in the following subsections. In each of the simulations, a disturbance in is applied.
7.2. Control of
is held constant at 35.3 L/d (cf. Table 2). On the real reactor, this can be implemented easily since the feed pump is a peristaltic pump which gives the demanded flow without feedback (flow) control.
Figure 9 shows the simulated response with constant . Table 4 shows performance measures.
(a)
(b)
(c)
(d)
7.3. Control of
is controlled to its setpoint, , using feedback from from the Kalman Filter (cf. Section 5). Both predictive control and PI control are tested.
7.3.1. Predictive Control
The optimization criterion of the predictive controller is selected as where with constraint . The control error is . Comparing with the criterion of predictive control of , (11), the term , which is at the end of the prediction horizon, is now included. The term brings approximately to zero. Without this term, is 0.1 g/L, and the control signal is actually constant. It is found that and are proper settings.
It is found that the predictive control is considerably smoother with d than with d which is used in Section 6. Increasing from 0.025 d, which is used in Section 6, to d, here, has very little impact on the control system performance over the simulation time interval used here, while the computational burden is noticeably less.
7.3.2. PI Control
PI controller is also applied. The PI settings are (L/d)/(g VFA/L) and d found using the Relaxed Ziegler-Nichols closed-loop method based on relay oscillations [38] which is a quick method to use on a simulator.
7.3.3. Simulations
The initial operating point is as shown in Table 2. The setpoint is g/L. At d, the disturbance is changed as a ramp of slope 2 (g/L)/d during 1 d, which is the same variation as in control (cf. Section 6). This is a reasonable variation for the real reactor. Measurement noise is not included in simulations.
Figure 9 shows simulated responses in , , , and with predictive control and PI control and with constant . Table 4 shows performance measures.
7.3.4. Results and Discussion
Table 4 shows performance and robustness measures with the three control strategies above. , defined by (13), is calculated over the simulated time interval. is the maximum control error. GM and PM are found as explained in Section 6.4.
Comments(i)In Figure 9 it is seen that the setpoint tracking works for both predictive control and PI control. However, the PI controller gives a more smooth response in .(ii)The lower-right plot in Figure 9 shows the manipulated which is adjusted by the controllers. The control action is smoother with PI control than with predictive control.(iii)The performance measures shown in Table 4 indicate that PI control of , based on feedback from UKF, is the best control strategy here.(iv)Also using a constant can be regarded as acceptable with the disturbance change simulated.(v)The upper-right plot in Figure 9 illustrates that is not under control. Although not shown here, settles at steady state at approximately d.(vi)GM is large with PI control. With predictive control, the notion of GM is questionable, since the simulated control system does not actually become unstable for any gain increase at the process input. Rather, the gain increase is seen by the UKF as an increase in . The relatively large estimate of is used in the prediction by the predictive controller, causing a large overshoot in before it eventually reaches (plots are not shown here). This behavior is the same as with predictive control of (cf. Section 6.4).(vii)PM is large with PI control. With predictive control, no limit was found; that is, the controller handles unmodeled time-delays in the controlled process even as large as 10 d.
7.4. Control of
7.4.1. Controllers
The third control strategy proposed in Section 7.1 is controlling to a setpoint, , set equal to the value of at the ultimate operating point (cf. Table 2). Both predictive control based on feedback from UKF estimates and PI control based on measurement of are simulated.
7.4.2. Simulations
The simulation scenario differs from the scenario of the simulations in Section 7.3 as is now decreased instead of increased. Decreasing is selected here because, in the corresponding response, increases (in steady state), and an increase of is more critical than a decrease.
In the predictive controller, is set as d, and is 1 d.
Figure 10 shows simulated responses in , , , and with predictive control based on feedback from UKF estimates and PI control based on measurement of .
(a)
(b)
(c)
(d)
7.4.3. Results
(i)As seen in Figure 10, is much closer to with predictive control than with PI control.(ii)With both predictive control and PI control, increases. Simulations over 400 d show that goes toward approximately 1.05 g/L, which is 0.25 larger than the critical value of 0.8 g/L. This makes this control strategy questionable.
7.5. What Is the Best Control Strategy for Safe Reactor Operation?
From the results in Sections 7.3, 7.2, and 7.4, it can be concluded that the best control strategy for safe reactor operation is controlling to a (fixed) setpoint using feedback from the state estimator (UKF). In the aforementioned control strategy, PI control is evaluated as better than predictive control. These two controllers give similar disturbance compensation, but the control signal is smoother with PI control than with predictive control.
8. Conclusions
The original four states of the modified Hill model, , , , , and the assumed unknown organic content, , of the feedstock of a real pilot AD reactor have been mainly successfully estimated with an Unscented Kalman Filter (UKF), but with an estimation error for in a part of the time interval.
These estimates, together with the model, have been applied in two different model-based control systems. The first system aims at retaining at a possibly time-varying setpoint, which may originate from a demanded power production by the reactor. Simulations indicate that the setpoint tracking performance of the predictive controller is considerably better while disturbance compensation, assuming that the disturbance has an unknown value, is not much better compared with PI control, confirming a well-known fact, compared to, for example, [33]. Consequently, assuming the setpoint is constant, the PI controller competes well with the predictive controller. A successful application of predictive control of the real reactor is reported.
The second control system aims at retaining the reactor at an ultimate safe operating point, where has a critical maximum value. This operating point is characterized by three corresponding values of , , and , as found from steady-state simulations of the reactor model. These operating point values can be used as setpoints in pertinent control systems. Simulations indicate that the best control solution among the three alternatives is PI control based on feedback of estimated by Kalman Filter.
The results of this paper indicate that a model-based control system, using a relatively simple mechanistic dynamical reactor model, can be designed and implemented on real AD reactors.
Abbreviations
| AD: | Anaerobic digestion |
| BVS: | Biodegradable volatile solids |
| COD: | Chemical oxygen demand |
| EKF: | Extended Kalman Filter |
| FC: | Flow controller |
| FT: | Flow transmitter (sensor) |
| HRT: | Hydraulic retention time |
| IAE: | Integral of absolute error |
| MPC: | Model-based predictive control |
| NTP: | Normal temperature and pressure—0°C, 1 atm |
| PI: | Proportional plus integral (control) |
| PWM: | Pulse width modulation |
| TIC: | Total inorganic carbon |
| UKF: | Unscented Kalman Filter |
| VFA: | Volatile fatty acids |
| VS: | Volatile solids. |
| : | Amplitude of the control error and the process output (measurement) |
| ((g VFA/L)/(g BVS/L)): | Acidity constant |
| : | Amplitude of the on-off control signal |
| ((g BVS/L)/(g VS/L)): | Biodegradability constant |
| : | Cost (weight) factor of in predictive control |
| : | Cost (weight) factor of in predictive control |
| (d−1): | Dilution rate |
| : | Control error |
| : | Objective function |
| (L/d): | Influent or feed flow or load rate, assumed equal to effluent flow (constant volume) |
| (L CH4/d): | Methane gas flow |
| (L CH4/d): | Setpoint of |
| GM: | Gain margin |
| : | Discrete-time index |
| (g BVS/(g acidogens/L)): | Yield constant |
| (g VFA/(g acidogens/L)): | Yield constant |
| (g VFA/(g methanogens/L)): | Yield constant |
| (L/g methanogens): | Yield constant |
| (g BVS/L): | Monod half-velocity constant for acidogens |
| (g VFA/L): | Monod half-velocity constant for methanogens |
| (d−1): | Specific death rate of acidogens |
| (d−1): | Specific death rate of methanogens |
| (d−1): | Reaction (growth) rate of acidogens |
| (d−1): | Reaction (growth) rate of methanogens |
| (d−1): | Maximum reaction rate for acidogens |
| (d−1): | Maximum reaction rate for methanogens |
| (kW) | Power |
| (d): | Period of oscillation |
| : | State estimation error covariance |
| PM (degrees): | Phase margin |
| : | Process noise covariance |
| : | Measurement noise covariance |
| (g VFA/L): | Concentration of VFA acids in reactor |
| (g VFA/L): | Estimate of |
| (g VFA/L): | Concentration of VFA in biodegradable part of influent |
| (g VFA/L): | Upper limit of safe range of concentration of VFA in reactor |
| (g VFA/L): | Setpoint of |
| (g BVS/L): | Concentration of BVS in reactor |
| (g BVS/L): | Concentration of BVS in influent |
| (g VS/L): | Concentration of volatile solids in influent |
| (L/d): | Standard deviation of control signal |
| (°C): | Reactor temperature |
| : | Prediction horizon |
| (d): | Controller integral time |
| (d): | Integration variable in predictive control criterion |
| (L): | Effective reactor volume (assumed filled with liquid) |
| : | Estimated state vector |
| (g acidogens/L): | Concentration of acidogens |
| (g methanogens/L): | Concentration of methanogens. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Funding of this project is supplied by Innovasjon Norge, Statens Landbruksforvaltning, Utdannings-og Forskningsdep., The Research Council of Norway, and Telemark University College (TUC). Thanks are due to E. Fjelddalen, W. Bergland, M. Torabzadegan, and students in the Environmental Engineering M.S. degree study at TUC for practical support.