# State Estimation and Model-Based Control of a Pilot Anaerobic Digestion Reactor

Authors: Finn Haugen; Rune Bakke; Bernt Lie
Journal: Journal of Control Science and Engineering
 (2014)

---

## 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.


---

## Body

## 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. Section2 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, Svsin, is described in Section 5. These estimates are used for control of Fmeth, which is described in Section 6, which includes both simulated and real results. The estimates are also used for control of Svfa, 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,Fmeth (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 (1) the extremum-seeking variable gain controller by Liu et al. [23], (2) the disturbance monitoring controller by Steyer et al. [24], and (3) 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.


## 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,Fmeth (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 (1) the extremum-seeking variable gain controller by Liu et al. [23], (2) the disturbance monitoring controller by Steyer et al. [24], and (3) 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

Figure1 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.Figure 1
Control system of the AD reactor.In Figure1, 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 Fmeth 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.F
feed 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,Treac is also a candidate as control variable since it has a clear impact on Fmeth, but in [28] we argue why Treac is not considered a usable control variable.An online measurement ofTreac is used by the controller, since Treac is an important model variable. Treac 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,Treac 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

Figure2 shows the structure of the control system.Figure 2
Block diagram of model-based optimization and control system. (Terms and variables are defined in the text.)In the block diagram:u=Ffeed, and d=Svsin. x comprises here the four state variables of the modified Hill model, compared with Section 3.3: x=[Sbvs,Svfa,Xacid,Xmeth]T. Depending on the applications in this paper, y=Fmeth, compared with Section 6, or y=Svfa, 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 symbolMa in various blocks in Figure 2 represents the assumed mathematical model used in the block. The Mr symbol in the Process block is the model representing the real system (process). Only if model errors are assumed to be zero, Ma and Mr will be identical.The connections fromd and/or dest to the Control Designer are due to d being an input to the process, and the value of d or dest is included in the model-based optimization. For example, the value of Svsin in the feed of the reactor has an impact on the specific value of Ffeed needed to produce a specified Fmeth 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 Figure2, 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:(1)S˙bvs=(B0Svsin-Sbvs)FfeedV-μk1XacidS˙vfa=(AfB0Svsin-Svfa)FfeedV+μk2Xacid-μck3XmethX˙acid=(μ-Kd-Ffeed/bV)XacidX˙meth=(μc-Kdc-Ffeed/bV)Xmeth;methane gas production:(2)Fmeth=Vμck5Xmeth;reaction rates:(3)μ=μmSbvsKs+Sbvsμc=μmcSvfaKsc+Svfaμm=μmc=0.013Treac-0.129for20°C<Treac<60°C.Table1 shows model parameter values as adapted to AD reactor at Foss Farm, [7].Table 1
Parameters in the modified Hill model adapted to the AD reactor at Foss Farm.
Parameter
Value
Unit
A
f
0.69
(g VFA/L)/(g BVS/L)
b
2.90
d/d
B
0
0.25
(g BVS/L)/(g VS/L)
k
1
3.89
g BVS/(g acidogens/L)
k
2
1.76
g VFA/(g acidogens/L)
k
3
31.7
g VFA/(g methanogens/L)
k
5
26.3
L/g methanogens
K
d
0.02
d−1
K
d
c
0.02
d−1
K
s
15.5
g BVS/L
K
s
c
3
g VFA/L
V
250
LOne example of a set of steady-state values of the AD process variables is given in Table2.Table 2
The ultimate safe steady-state operating point.
Variable
Value
Unit
S
vfa
0.8
g/L
F
feed
35.3
L/d
F
feed
/
V
=
D
0.14
(L/d)/L
V
/
F
feed
= HRT
7.1
d
F
meth
174.2
L CH4/d
F
meth
/
V
0.70
(L CH4/d)/L
S
bvs
4.14
g/L
X
acid
1.80
g/L
X
meth
0.39
g/L
S
v
s
in
30.2
g/L
T
reac
35
°C


## 3.1. AD Reactor with Control System

Figure1 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.Figure 1
Control system of the AD reactor.In Figure1, 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 Fmeth 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.F
feed 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,Treac is also a candidate as control variable since it has a clear impact on Fmeth, but in [28] we argue why Treac is not considered a usable control variable.An online measurement ofTreac is used by the controller, since Treac is an important model variable. Treac 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,Treac 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

Figure2 shows the structure of the control system.Figure 2
Block diagram of model-based optimization and control system. (Terms and variables are defined in the text.)In the block diagram:u=Ffeed, and d=Svsin. x comprises here the four state variables of the modified Hill model, compared with Section 3.3: x=[Sbvs,Svfa,Xacid,Xmeth]T. Depending on the applications in this paper, y=Fmeth, compared with Section 6, or y=Svfa, 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 symbolMa in various blocks in Figure 2 represents the assumed mathematical model used in the block. The Mr symbol in the Process block is the model representing the real system (process). Only if model errors are assumed to be zero, Ma and Mr will be identical.The connections fromd and/or dest to the Control Designer are due to d being an input to the process, and the value of d or dest is included in the model-based optimization. For example, the value of Svsin in the feed of the reactor has an impact on the specific value of Ffeed needed to produce a specified Fmeth 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 Figure2, 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:(1)S˙bvs=(B0Svsin-Sbvs)FfeedV-μk1XacidS˙vfa=(AfB0Svsin-Svfa)FfeedV+μk2Xacid-μck3XmethX˙acid=(μ-Kd-Ffeed/bV)XacidX˙meth=(μc-Kdc-Ffeed/bV)Xmeth;methane gas production:(2)Fmeth=Vμck5Xmeth;reaction rates:(3)μ=μmSbvsKs+Sbvsμc=μmcSvfaKsc+Svfaμm=μmc=0.013Treac-0.129for20°C<Treac<60°C.Table1 shows model parameter values as adapted to AD reactor at Foss Farm, [7].Table 1
Parameters in the modified Hill model adapted to the AD reactor at Foss Farm.
Parameter
Value
Unit
A
f
0.69
(g VFA/L)/(g BVS/L)
b
2.90
d/d
B
0
0.25
(g BVS/L)/(g VS/L)
k
1
3.89
g BVS/(g acidogens/L)
k
2
1.76
g VFA/(g acidogens/L)
k
3
31.7
g VFA/(g methanogens/L)
k
5
26.3
L/g methanogens
K
d
0.02
d−1
K
d
c
0.02
d−1
K
s
15.5
g BVS/L
K
s
c
3
g VFA/L
V
250
LOne example of a set of steady-state values of the AD process variables is given in Table2.Table 2
The ultimate safe steady-state operating point.
Variable
Value
Unit
S
vfa
0.8
g/L
F
feed
35.3
L/d
F
feed
/
V
=
D
0.14
(L/d)/L
V
/
F
feed
= HRT
7.1
d
F
meth
174.2
L CH4/d
F
meth
/
V
0.70
(L CH4/d)/L
S
bvs
4.14
g/L
X
acid
1.80
g/L
X
meth
0.39
g/L
S
v
s
in
30.2
g/L
T
reac
35
°C


## 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 Svfa>0.8 g/L indicates an impending reactor failure, causing a reduction of methane production. Hence, it is here stated that
(4)Svfa≤0.8g/L=Svfamax⁡
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 Svfamax⁡. However, a new value of Svfamax⁡ 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 Svfa 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.Figure3 shows simulated static (steady-state) responses in a number of variables to a range of constant feed rates (Ffeed). The cyan horizontal line in the Svfa plot represents Svfa=0.8 g/L. The green intervals on the abscissas indicate safe reactor operation, and, conversely, the red interval indicates unsafe operation.Simulated static (steady-state) values of a number of variables versusFfeed (constant) at Treac=35°C. The green intervals on the abscissas indicate safe reactor operation as defined inequality (4). Conversely, the red intervals represent unsafe reactor operation.
(a)
(b)
(c)
(d)
(e)
(f)Table2 shows the values of several variables at the ultimate safe steady-state operating point. The set of three corresponding values (Svfa, Ffeed, and Fmeth) 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 includeSvfamax⁡=0.8 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 Svfamax⁡ 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 Sections6 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. Section3.3) Sbvs, Svfa, Xacid, and Xmeth. They are estimated with the UKF. It is decided to also estimate Svsin with the UKF since it is assumed that its value may vary, though slowly. As is common, Svsin is modeled as a ‘‘random walk”: S˙vsin=w, where w is a random disturbance. Thus, the augmented state vector to be estimated by the UKF is
(5)x=[Sbvs,Svfa,Xacid,Xmeth,Svsin]T.F
feed is regarded as an input variable to the UKF. Ffeed is the control variable, which is always known.The model parameters are known from model adaptation [7]. Treac may vary but is always known as it is measured continuously.The process measurement,y, used by the UKF is Fmeth available from sensor FT in Figure 1. Hence, y=Fmeth in the UKF.


### 5.2. Observability

The linearized reactor model, augmented withSvsin, is found observable at a number of typical operating points using theobsv 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:x^(t0∣t0) (initial estimated state; the initial a posteriori estimate), P^(t0∣t0) (initial state estimation error covariance), R (measurement noise covariance), and Q (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)
x^(t0∣t0) is set equal to the values from laboratory analysis at the start of the pertinent time interval. This applies ideally to Sbvs, Svfa, and Svsin. However, for Svsin, we impose, for the purpose of demonstration, a large initial estimation error, by setting the initial estimate of Svsin equal to 20% of the value known from laboratory analysis.Xacid and Xmethare not known, but their initial values are calculated from the model assuming steady state (details of the calculation can be found in [7]).(ii)
P^(t0∣t0) is set as a diagonal matrix as follows:
(6)P^ii(t0t0)=[kPx^i(t0t0)]2,
with kP=0.01.(iii)
R is a diagonal matrix, which, since the number of measurements is one (Fmeth), is reduced to a scalar—the measurement variance. From a representative real time series,
(7)var(Fmeth)=1.44=R.(iv)
Q is typically set as a constant matrix (diagonal). Assuming that x^(t0∣t0), P^(t0∣t0), and R are set, Q can be used as final tuning parameter.
(a)
IncreasingQi,i makes the estimate for state variable xi converge faster to the assumed true value, but with the drawback that the estimate for xi becomes more noisy (caused by the increased propagation of the measurement noise, via the Kalman Filter gain(s)).(b)
ReducingQi,i 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:(8)Qi,i=[kQmixi(t0∣t0)]2.With the initial setting ofmi=1, it is found that kQ=0.0005 is a proper value. Then the ultimate tuning is made by adjusting mi. By trial-and-error, {mi}={10,1,1,1,10}.


### 5.4. Results and Discussion

Figure4 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 Fmeth.Estimates and real data from laboratory analysis with UKF. Standard deviations of estimation errors:Sbvs: 0.62 g/L; Svfa: 0.32 g/L; Svsin: 1.02 g/L. (For Svsin, the standard deviation is calculated from t=120 d, since a relatively large initial estimation error is imposed on purpose.)
(a)
(b)
(c)
(d)
(e)
(f)Overall, the UKF gives reasonably good estimates (real values ofXacid and Xmeth are not known).The large initial estimation errorSvsin imposed on purpose is effectively reduced during approximately 15 days.Fromt=150 d, there is a noticeable difference between the estimate and the laboratory analysis of Svfa. 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.


## 5.1. Variables and Parameters of the Model

The state variables of the modified Hill model are (cf. Section3.3) Sbvs, Svfa, Xacid, and Xmeth. They are estimated with the UKF. It is decided to also estimate Svsin with the UKF since it is assumed that its value may vary, though slowly. As is common, Svsin is modeled as a ‘‘random walk”: S˙vsin=w, where w is a random disturbance. Thus, the augmented state vector to be estimated by the UKF is
(5)x=[Sbvs,Svfa,Xacid,Xmeth,Svsin]T.F
feed is regarded as an input variable to the UKF. Ffeed is the control variable, which is always known.The model parameters are known from model adaptation [7]. Treac may vary but is always known as it is measured continuously.The process measurement,y, used by the UKF is Fmeth available from sensor FT in Figure 1. Hence, y=Fmeth in the UKF.


## 5.2. Observability

The linearized reactor model, augmented withSvsin, is found observable at a number of typical operating points using theobsv 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:x^(t0∣t0) (initial estimated state; the initial a posteriori estimate), P^(t0∣t0) (initial state estimation error covariance), R (measurement noise covariance), and Q (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)
x^(t0∣t0) is set equal to the values from laboratory analysis at the start of the pertinent time interval. This applies ideally to Sbvs, Svfa, and Svsin. However, for Svsin, we impose, for the purpose of demonstration, a large initial estimation error, by setting the initial estimate of Svsin equal to 20% of the value known from laboratory analysis.Xacid and Xmethare not known, but their initial values are calculated from the model assuming steady state (details of the calculation can be found in [7]).(ii)
P^(t0∣t0) is set as a diagonal matrix as follows:
(6)P^ii(t0t0)=[kPx^i(t0t0)]2,
with kP=0.01.(iii)
R is a diagonal matrix, which, since the number of measurements is one (Fmeth), is reduced to a scalar—the measurement variance. From a representative real time series,
(7)var(Fmeth)=1.44=R.(iv)
Q is typically set as a constant matrix (diagonal). Assuming that x^(t0∣t0), P^(t0∣t0), and R are set, Q can be used as final tuning parameter.
(a)
IncreasingQi,i makes the estimate for state variable xi converge faster to the assumed true value, but with the drawback that the estimate for xi becomes more noisy (caused by the increased propagation of the measurement noise, via the Kalman Filter gain(s)).(b)
ReducingQi,i 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:(8)Qi,i=[kQmixi(t0∣t0)]2.With the initial setting ofmi=1, it is found that kQ=0.0005 is a proper value. Then the ultimate tuning is made by adjusting mi. By trial-and-error, {mi}={10,1,1,1,10}.


## 5.4. Results and Discussion

Figure4 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 Fmeth.Estimates and real data from laboratory analysis with UKF. Standard deviations of estimation errors:Sbvs: 0.62 g/L; Svfa: 0.32 g/L; Svsin: 1.02 g/L. (For Svsin, the standard deviation is calculated from t=120 d, since a relatively large initial estimation error is imposed on purpose.)
(a)
(b)
(c)
(d)
(e)
(f)Overall, the UKF gives reasonably good estimates (real values ofXacid and Xmeth are not known).The large initial estimation errorSvsin imposed on purpose is effectively reduced during approximately 15 days.Fromt=150 d, there is a noticeable difference between the estimate and the laboratory analysis of Svfa. 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 ofFmeth, Figure 5 shows, for the real pilot reactor, experimental time-series of Fmeth and Ffeed (and Treac) with feedback control and without control. It is clear that Fmeth varies less with control than without control. Fmeth remains close to Fmethsp even after the setpoint is changed. The variations are due to inevitable disturbances. In the case of feedback control, Ffeed is of course varying, while it is constant in the case of no control (i.e., open-loop control). Treac is actually different in the two cases, but it is assumed that the difference between the two cases is independent of the temperature difference.F
meth and Ffeed and Treac, for the real reactor, with (automatic) control and without control. (The length of each of the time intervals for these two cases is different.) (Reprinted from [28] by permission.)
(a)
(b)
(c)
(d)Whether the variation inFmeth 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 ofFmeth exists. The operational objective is stated as producing a demanded methane gas flow. A specific value of Fmeth is related to the power, P (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 Figure3.Furthermore, safe reactor operation must be ensured, which here means that inequality (4) is satisfied.Relating to Figure2, the above specifications concerning Fmeth, the limitation of variations of Ffeed, and the condition inequality (4) are inputs to the Control Designer. Outputs from the Control Designer are Fmethsp and Cdu. 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 forFmeth 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 Section3.2, a predictive controller is used as advanced controller in this paper. A predictive controller to retain Fmeth 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=0.2 d is included at the control input of the model:
(9)Ffeed(t)=u(t-τd),
where Ffeed is the feed rate of the modified Hill model and u 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
(10)min⁡ufobj,
where
(11)fobj=∫tt+τh‍[e2(θ)+Cduu˙2(θ)]dθ
with constraint umin⁡≤u(t)≤umax⁡ which is included in the optimization problem formulation; that is, it is an input argument in the fmincon function call in Matlab. t is the present time instance. e is the control error, e=Fmethsp-Fmeth. The time derivative, u˙, represents the control signal changes. The larger the Cdu, the smoother the control actions.In implementations, the discretized version offobj is minimized, giving an optimal control sequence, {u}opt, over the prediction horizon. The first element of this sequence, that is, (t0)opt, is applied as control signal at the present time point. The prediction horizon is receding and the procedure of obtaining {u}opt and u(t0)opt 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.{u}opt is calculated with the nonlinear optimization function fmincon in the Optimization toolbox of Matlab. The present state, x(t), 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τs=0.025 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 fobj. τh corresponds to 1/0.025=40 time-steps, which is then the prediction horizon in number of time-steps.C
du
=
0.01 in (11) is found by trial-and-error on a simulator. A proper value of the prediction horizon is found as τh=1 d (with τh<0.5 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 Table2 using the Skogestad method [37], with the modification of the τi setting as proposed in [38]. The PI settings are Kc=0.89 [(L CH4/d)/(L feed/d)] and τi=0.8 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. IAEs measures the setpoint tracking:
(12)IAEs=∫tistfs‍|e|dt.TheIAEd measures the disturbance compensation:
(13)IAEd=∫tidtfd‍|e|dt.(
2) Control Signal Variations (Performance). As measures of the variation of the control signal, both the standard deviation, σu, and the mean of the absolute value of the rate of changes, μ|u˙|, 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,ΔK, is inserted into the loop (between the controller and the process); see Figure 6. Normally, ΔK=1. The (ultimate) value ΔKu that brings the (simulated) control system to the stability limit, with sustained oscillations, is found by trials. Then,
(14)GM=ΔKu.Figure 6
An adjustable gain and time-delay are inserted into the loop to find the stability margins (gain margin and phase margin) experimentally.To calculate the PM, an adjustable time-delay,Δτdelay, is inserted into the loop; see Figure 6. Normally, Δτdelay=0. The value Δτdelayu 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 Pu [s]. As shown in [39] (Appendix  1),
(15)PM[deg⁡]=360°ΔτdelayuPu.
Seborg et al. [40] propose the following ranges for appropriate values of the stability margins: 1.7=4.6 dB ≤GM≤4.0=12.0 dB and 30°≤PM≤45°.Relating to Figure2, ΔK and Δτdelay are included before the Process block, after the branch from u to the Estimator.


#### 6.4.3. Simulations

Figure7 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 Fmethsp is varied as a sequence of two ramps of slope ±2 (L CH4/d)/d each lasting for 1 d. The disturbance Svsin is varied as a ramp of slope 2 (g/L)/d during 1 d which is a realistic variation [28].Simulated responses with predictive control and PI control.
(a)
(b)
(c)
(d)The simulations shown in Figure7 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 σn=1.2 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

Table3 shows performance and robustness measures with predictive control and with PI control. The IAE indexes, (12) and (13), are calculated with tis=1 d, tfs=7 d, tid=7 d, and tfd=14 d.Table 3
Performance and robustness measures for predictive control and PI control.
Predictive
PI
Ratio Pred./PI
IAE
s
0.093
0.72
0.13
IAE
d
1.46
2.20
0.66
σ
u
0.69
1.29
0.67
μ
|
u
˙
|
7.07
25.2
0.28
GM
≈

∞
2.6
N/A
PM
63.9°
47.6°
1.34Comments on the results shown in Table3 are the following.(i)
IAEs with predictive control is 13% of IAEs with PI control. Hence, predictive control is clearly the best.(ii)
IAEd with predictive control is 66% of IAEd with PI control. Again, predictive control is the best, but the improvement compared with PI control is not large.(iii)
σu with predictive control is approximately 67% of the value with PI control, while μ|u˙| 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 (reducingKc and increasing τi according to Skogestad’s formulas), both σu and μ|u˙| 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 μ|u˙| with predictive control and PI control was approximately equal. The IAEd 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 inSvsin. Consequently, the estimate of Svsin is increased, which in turn is used in the prediction by the predictive controller, causing a large overshoot or undershoot in Fmeth before it eventually reaches Fmethsp (plots of simulations not shown). From simulations it is found that Fmethsp is back at its setpoint during 1-2 d for 0.5≤ΔKu≤4.(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 Section6.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,τs=0.05d and τp=2 d. In the simulations in Section 6.4, τs=0.025 d and τp=1 d. τs=0.05 d has been tested in simulations, giving a slight change in performance, probably due to less accurate numerical integration (explicit Euler is used). Np=τp/τs=40 is the same both in the practical experiments and in the simulations.(ii)
Ffeed 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 factorCdu 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 Cdu in the simulations may be due to dynamic phenomena of the real reactor not encapsulated by the model. In any case, Cdu is typically a tuning parameter.


#### 6.5.1. Results and Discussion

Figure8 shows the time-series of the practical experiments. Below are comments on this figure.(i)
Att=99.8 d, Fmethsp 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)
Att=102.3 d, a preset-ramped setpoint profile started. The predictive controller adjusts Ffeed before Fmethsp starts increasing. The tracking is accurate. The upper bound of Ffeed of 40 L/d is eventually reached, causing the rate of change of Fmeth to become less than the rate of change of Fmethsp.(iii)
Att=102.8 d, the rate of change of Fmethsp is instantly adjusted from +20  (L CH4/d)/d to -20  (L CH4/d)/d. The observed lag in Fmeth can be explained with the instant change of Fmethsp which prevents predictive control action.(iv)
Att=104.1 d, a preset step change of Fmethsp from 150 to 155 L CH4/d is applied. The predictive adjustment of Ffeed is obvious. Fmeth shows a clear overshoot, but it is expected that the response will stabilize.(v)
Att=104.4 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 Fmeth eventually would have been brought back to its setpoint.Time series from application of predictive control ofFmeth on the real reactor.
(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 Figure3. According to this plot, the setpoint values used in the experiments (cf. Figure 8) are actually feasible.


## 6.1. The Effect of Feedback Control

To demonstrate the effect of feedback (or automatic or closed-loop) control ofFmeth, Figure 5 shows, for the real pilot reactor, experimental time-series of Fmeth and Ffeed (and Treac) with feedback control and without control. It is clear that Fmeth varies less with control than without control. Fmeth remains close to Fmethsp even after the setpoint is changed. The variations are due to inevitable disturbances. In the case of feedback control, Ffeed is of course varying, while it is constant in the case of no control (i.e., open-loop control). Treac is actually different in the two cases, but it is assumed that the difference between the two cases is independent of the temperature difference.F
meth and Ffeed and Treac, for the real reactor, with (automatic) control and without control. (The length of each of the time intervals for these two cases is different.) (Reprinted from [28] by permission.)
(a)
(b)
(c)
(d)Whether the variation inFmeth 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 ofFmeth exists. The operational objective is stated as producing a demanded methane gas flow. A specific value of Fmeth is related to the power, P (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 Figure3.Furthermore, safe reactor operation must be ensured, which here means that inequality (4) is satisfied.Relating to Figure2, the above specifications concerning Fmeth, the limitation of variations of Ffeed, and the condition inequality (4) are inputs to the Control Designer. Outputs from the Control Designer are Fmethsp and Cdu. 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 forFmeth 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 Section3.2, a predictive controller is used as advanced controller in this paper. A predictive controller to retain Fmeth 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=0.2 d is included at the control input of the model:
(9)Ffeed(t)=u(t-τd),
where Ffeed is the feed rate of the modified Hill model and u 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
(10)min⁡ufobj,
where
(11)fobj=∫tt+τh‍[e2(θ)+Cduu˙2(θ)]dθ
with constraint umin⁡≤u(t)≤umax⁡ which is included in the optimization problem formulation; that is, it is an input argument in the fmincon function call in Matlab. t is the present time instance. e is the control error, e=Fmethsp-Fmeth. The time derivative, u˙, represents the control signal changes. The larger the Cdu, the smoother the control actions.In implementations, the discretized version offobj is minimized, giving an optimal control sequence, {u}opt, over the prediction horizon. The first element of this sequence, that is, (t0)opt, is applied as control signal at the present time point. The prediction horizon is receding and the procedure of obtaining {u}opt and u(t0)opt 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.{u}opt is calculated with the nonlinear optimization function fmincon in the Optimization toolbox of Matlab. The present state, x(t), 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τs=0.025 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 fobj. τh corresponds to 1/0.025=40 time-steps, which is then the prediction horizon in number of time-steps.C
du
=
0.01 in (11) is found by trial-and-error on a simulator. A proper value of the prediction horizon is found as τh=1 d (with τh<0.5 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 Table2 using the Skogestad method [37], with the modification of the τi setting as proposed in [38]. The PI settings are Kc=0.89 [(L CH4/d)/(L feed/d)] and τi=0.8 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. IAEs measures the setpoint tracking:
(12)IAEs=∫tistfs‍|e|dt.TheIAEd measures the disturbance compensation:
(13)IAEd=∫tidtfd‍|e|dt.(
2) Control Signal Variations (Performance). As measures of the variation of the control signal, both the standard deviation, σu, and the mean of the absolute value of the rate of changes, μ|u˙|, 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,ΔK, is inserted into the loop (between the controller and the process); see Figure 6. Normally, ΔK=1. The (ultimate) value ΔKu that brings the (simulated) control system to the stability limit, with sustained oscillations, is found by trials. Then,
(14)GM=ΔKu.Figure 6
An adjustable gain and time-delay are inserted into the loop to find the stability margins (gain margin and phase margin) experimentally.To calculate the PM, an adjustable time-delay,Δτdelay, is inserted into the loop; see Figure 6. Normally, Δτdelay=0. The value Δτdelayu 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 Pu [s]. As shown in [39] (Appendix  1),
(15)PM[deg⁡]=360°ΔτdelayuPu.
Seborg et al. [40] propose the following ranges for appropriate values of the stability margins: 1.7=4.6 dB ≤GM≤4.0=12.0 dB and 30°≤PM≤45°.Relating to Figure2, ΔK and Δτdelay are included before the Process block, after the branch from u to the Estimator.


### 6.4.3. Simulations

Figure7 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 Fmethsp is varied as a sequence of two ramps of slope ±2 (L CH4/d)/d each lasting for 1 d. The disturbance Svsin is varied as a ramp of slope 2 (g/L)/d during 1 d which is a realistic variation [28].Simulated responses with predictive control and PI control.
(a)
(b)
(c)
(d)The simulations shown in Figure7 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 σn=1.2 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

Table3 shows performance and robustness measures with predictive control and with PI control. The IAE indexes, (12) and (13), are calculated with tis=1 d, tfs=7 d, tid=7 d, and tfd=14 d.Table 3
Performance and robustness measures for predictive control and PI control.
Predictive
PI
Ratio Pred./PI
IAE
s
0.093
0.72
0.13
IAE
d
1.46
2.20
0.66
σ
u
0.69
1.29
0.67
μ
|
u
˙
|
7.07
25.2
0.28
GM
≈

∞
2.6
N/A
PM
63.9°
47.6°
1.34Comments on the results shown in Table3 are the following.(i)
IAEs with predictive control is 13% of IAEs with PI control. Hence, predictive control is clearly the best.(ii)
IAEd with predictive control is 66% of IAEd with PI control. Again, predictive control is the best, but the improvement compared with PI control is not large.(iii)
σu with predictive control is approximately 67% of the value with PI control, while μ|u˙| 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 (reducingKc and increasing τi according to Skogestad’s formulas), both σu and μ|u˙| 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 μ|u˙| with predictive control and PI control was approximately equal. The IAEd 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 inSvsin. Consequently, the estimate of Svsin is increased, which in turn is used in the prediction by the predictive controller, causing a large overshoot or undershoot in Fmeth before it eventually reaches Fmethsp (plots of simulations not shown). From simulations it is found that Fmethsp is back at its setpoint during 1-2 d for 0.5≤ΔKu≤4.(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.4.1. Controller Settings

The settings of the predictive controller in the simulations are as follows.A time-step ofτs=0.025 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 fobj. τh corresponds to 1/0.025=40 time-steps, which is then the prediction horizon in number of time-steps.C
du
=
0.01 in (11) is found by trial-and-error on a simulator. A proper value of the prediction horizon is found as τh=1 d (with τh<0.5 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 Table2 using the Skogestad method [37], with the modification of the τi setting as proposed in [38]. The PI settings are Kc=0.89 [(L CH4/d)/(L feed/d)] and τi=0.8 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. IAEs measures the setpoint tracking:
(12)IAEs=∫tistfs‍|e|dt.TheIAEd measures the disturbance compensation:
(13)IAEd=∫tidtfd‍|e|dt.(
2) Control Signal Variations (Performance). As measures of the variation of the control signal, both the standard deviation, σu, and the mean of the absolute value of the rate of changes, μ|u˙|, 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,ΔK, is inserted into the loop (between the controller and the process); see Figure 6. Normally, ΔK=1. The (ultimate) value ΔKu that brings the (simulated) control system to the stability limit, with sustained oscillations, is found by trials. Then,
(14)GM=ΔKu.Figure 6
An adjustable gain and time-delay are inserted into the loop to find the stability margins (gain margin and phase margin) experimentally.To calculate the PM, an adjustable time-delay,Δτdelay, is inserted into the loop; see Figure 6. Normally, Δτdelay=0. The value Δτdelayu 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 Pu [s]. As shown in [39] (Appendix  1),
(15)PM[deg⁡]=360°ΔτdelayuPu.
Seborg et al. [40] propose the following ranges for appropriate values of the stability margins: 1.7=4.6 dB ≤GM≤4.0=12.0 dB and 30°≤PM≤45°.Relating to Figure2, ΔK and Δτdelay are included before the Process block, after the branch from u to the Estimator.


## 6.4.3. Simulations

Figure7 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 Fmethsp is varied as a sequence of two ramps of slope ±2 (L CH4/d)/d each lasting for 1 d. The disturbance Svsin is varied as a ramp of slope 2 (g/L)/d during 1 d which is a realistic variation [28].Simulated responses with predictive control and PI control.
(a)
(b)
(c)
(d)The simulations shown in Figure7 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 σn=1.2 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

Table3 shows performance and robustness measures with predictive control and with PI control. The IAE indexes, (12) and (13), are calculated with tis=1 d, tfs=7 d, tid=7 d, and tfd=14 d.Table 3
Performance and robustness measures for predictive control and PI control.
Predictive
PI
Ratio Pred./PI
IAE
s
0.093
0.72
0.13
IAE
d
1.46
2.20
0.66
σ
u
0.69
1.29
0.67
μ
|
u
˙
|
7.07
25.2
0.28
GM
≈

∞
2.6
N/A
PM
63.9°
47.6°
1.34Comments on the results shown in Table3 are the following.(i)
IAEs with predictive control is 13% of IAEs with PI control. Hence, predictive control is clearly the best.(ii)
IAEd with predictive control is 66% of IAEd with PI control. Again, predictive control is the best, but the improvement compared with PI control is not large.(iii)
σu with predictive control is approximately 67% of the value with PI control, while μ|u˙| 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 (reducingKc and increasing τi according to Skogestad’s formulas), both σu and μ|u˙| 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 μ|u˙| with predictive control and PI control was approximately equal. The IAEd 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 inSvsin. Consequently, the estimate of Svsin is increased, which in turn is used in the prediction by the predictive controller, causing a large overshoot or undershoot in Fmeth before it eventually reaches Fmethsp (plots of simulations not shown). From simulations it is found that Fmethsp is back at its setpoint during 1-2 d for 0.5≤ΔKu≤4.(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 Section6.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,τs=0.05d and τp=2 d. In the simulations in Section 6.4, τs=0.025 d and τp=1 d. τs=0.05 d has been tested in simulations, giving a slight change in performance, probably due to less accurate numerical integration (explicit Euler is used). Np=τp/τs=40 is the same both in the practical experiments and in the simulations.(ii)
Ffeed 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 factorCdu 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 Cdu in the simulations may be due to dynamic phenomena of the real reactor not encapsulated by the model. In any case, Cdu is typically a tuning parameter.


### 6.5.1. Results and Discussion

Figure8 shows the time-series of the practical experiments. Below are comments on this figure.(i)
Att=99.8 d, Fmethsp 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)
Att=102.3 d, a preset-ramped setpoint profile started. The predictive controller adjusts Ffeed before Fmethsp starts increasing. The tracking is accurate. The upper bound of Ffeed of 40 L/d is eventually reached, causing the rate of change of Fmeth to become less than the rate of change of Fmethsp.(iii)
Att=102.8 d, the rate of change of Fmethsp is instantly adjusted from +20  (L CH4/d)/d to -20  (L CH4/d)/d. The observed lag in Fmeth can be explained with the instant change of Fmethsp which prevents predictive control action.(iv)
Att=104.1 d, a preset step change of Fmethsp from 150 to 155 L CH4/d is applied. The predictive adjustment of Ffeed is obvious. Fmeth shows a clear overshoot, but it is expected that the response will stabilize.(v)
Att=104.4 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 Fmeth eventually would have been brought back to its setpoint.Time series from application of predictive control ofFmeth on the real reactor.
(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 Figure3. According to this plot, the setpoint values used in the experiments (cf. Figure 8) are actually feasible.


## 6.5.1. Results and Discussion

Figure8 shows the time-series of the practical experiments. Below are comments on this figure.(i)
Att=99.8 d, Fmethsp 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)
Att=102.3 d, a preset-ramped setpoint profile started. The predictive controller adjusts Ffeed before Fmethsp starts increasing. The tracking is accurate. The upper bound of Ffeed of 40 L/d is eventually reached, causing the rate of change of Fmeth to become less than the rate of change of Fmethsp.(iii)
Att=102.8 d, the rate of change of Fmethsp is instantly adjusted from +20  (L CH4/d)/d to -20  (L CH4/d)/d. The observed lag in Fmeth can be explained with the instant change of Fmethsp which prevents predictive control action.(iv)
Att=104.1 d, a preset step change of Fmethsp from 150 to 155 L CH4/d is applied. The predictive adjustment of Ffeed is obvious. Fmeth shows a clear overshoot, but it is expected that the response will stabilize.(v)
Att=104.4 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 Fmeth eventually would have been brought back to its setpoint.Time series from application of predictive control ofFmeth on the real reactor.
(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 Figure3. 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 Table2 (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)
Ffeed is controlled to a setpoint of Ffeedsp, which is35.3 L/d, assuming the operating point shown in Table 2. This control strategy is described in Section 7.2.(2)
Svfa is controlled to a setpoint of Svfasp, which is 0.8 g/L according to Table 2. In principle, this control requires feedback from the measurement of Svfa. Such sensors do exist [41, 42], but they are not in use on the present reactor. Instead, the estimate of Svfa calculated continuously with a state estimator (Kalman Filter) is used (cf. Section 5). This control strategy is described in Section 7.3.(3)
Fmeth is controlled to a setpoint of Fmethsp, which is174 L CH4/d according to Table 2. This control requires feedback from the measurement of Fmeth. 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,u=Ffeed (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 inSvsin is applied.


### 7.2. Control ofFfeed

F
feed 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.Figure9 shows the simulated response with constant Ffeed. Table 4 shows performance measures.Table 4
Performance measures for three control strategies for controllingSvfa.
Measure
Const.Ffeed
Predictive
PI
IAE
d
1.52
0.20
0.090
GM
N/A
≈

∞
5.5
PM
N/A
≈

∞
71.0°Control ofSvfa with different control functions: constant Ffeed, predictive control (MPC), and PI control. Simulated responses in Svfa, Fmeth, Ffeed, and Svsin are shown.
(a)
(b)
(c)
(d)


### 7.3. Control ofSvfa

S
vfa is controlled to its setpoint, Svfasp=Svfamax⁡, using feedback from Svfaest 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(16)min⁡ufobj,
where
(17)fobj=∫tt+τh‍[e2(θ)+Cduu˙2(θ)]dθ+Che2(τh)
with constraint umin⁡≤u(t)≤umax⁡. The control error is e=Svfasp-Svfaest. Comparing with the criterion of predictive control of Fmeth, (11), the term e2(τh), which is e2 at the end of the prediction horizon, is now included. The term brings e(t+τh) approximately to zero. Without this term, e(t+τh) is 0.1 g/L, and the control signal is actually constant. It is found that Cdu=0.2 and Ch=20 are proper settings.It is found that the predictive control is considerably smoother withτh=4d than with τh=1 d which is used in Section 6. Increasing τs from 0.025 d, which is used in Section 6, to 0.1 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 areKc=50.9 (L/d)/(g VFA/L) and τi=0.9 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 Table2. The setpoint is Svfasp=Svfamax⁡=0.8 g/L. At t=10 d, the disturbance Svsin is changed as a ramp of slope 2 (g/L)/d during 1 d, which is the same variation as in Fmeth control (cf. Section 6). This is a reasonable variation for the real reactor. Measurement noise is not included in simulations.Figure9 shows simulated responses in Svfa, Fmeth, Ffeed, and Svsin with predictive control and PI control and with constant Ffeed. Table 4 shows performance measures.


#### 7.3.4. Results and Discussion

Table4 shows performance and robustness measures with the three control strategies above. IAEd, defined by (13), is calculated over the simulated time interval. |e|max⁡ is the maximum control error. GM and PM are found as explained in Section 6.4.Comments
(i)
In Figure9 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 Svfa.
(ii)
The lower-right plot in Figure9 shows the manipulated Ffeed which is adjusted by the controllers. The control action is smoother with PI control than with predictive control.
(iii)
The performance measures shown in Table4 indicate that PI control of Svfa, based on feedback from UKF, is the best control strategy here.
(iv)
Also using a constantFfeed can be regarded as acceptable with the disturbance change simulated.
(v)
The upper-right plot in Figure9 illustrates that Fmeth is not under control. Although not shown here, Fmeth settles at steady state at approximately t=120 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 inSvsin. The relatively large estimate of Svsin is used in the prediction by the predictive controller, causing a large overshoot in Fmeth before it eventually reaches Fmethsp (plots are not shown here). This behavior is the same as with predictive control of Fmeth (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 ofFmeth

#### 7.4.1. Controllers

The third control strategy proposed in Section7.1 is controlling Fmeth to a setpoint, Fmethsp, set equal to the value of Fmeth at the ultimate operating point (cf. Table 2). Both predictive control based on feedback from UKF estimates and PI control based on measurement of Fmeth are simulated.


#### 7.4.2. Simulations

The simulation scenario differs from the scenario of the simulations in Section7.3 as Svsin is now decreased instead of increased. Decreasing Svsin is selected here because, in the corresponding response, Svfa increases (in steady state), and an increase of Svfa is more critical than a decrease.In the predictive controller,τs is set as0.05 d, and τh is 1 d.Figure10 shows simulated responses in Svfa, Fmeth, Ffeed, and Svsin with predictive control based on feedback from UKF estimates and PI control based on measurement of Fmeth.Predictive control and PI control ofFmeth. The process perturbation is a change in Svsin. Simulated responses in Svfa, Fmeth, Ffeed, and Svsin are shown.
(a)
(b)
(c)
(d)


#### 7.4.3. Results

(i)
As seen in Figure10, Fmeth is much closer to Fmethsp with predictive control than with PI control.
(ii)
With both predictive control and PI control,Svfa increases. Simulations over 400 d show that Svfa 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 Sections7.3, 7.2, and 7.4, it can be concluded that the best control strategy for safe reactor operation is controlling Svfa 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.


## 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 Table2 (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)
Ffeed is controlled to a setpoint of Ffeedsp, which is35.3 L/d, assuming the operating point shown in Table 2. This control strategy is described in Section 7.2.(2)
Svfa is controlled to a setpoint of Svfasp, which is 0.8 g/L according to Table 2. In principle, this control requires feedback from the measurement of Svfa. Such sensors do exist [41, 42], but they are not in use on the present reactor. Instead, the estimate of Svfa calculated continuously with a state estimator (Kalman Filter) is used (cf. Section 5). This control strategy is described in Section 7.3.(3)
Fmeth is controlled to a setpoint of Fmethsp, which is174 L CH4/d according to Table 2. This control requires feedback from the measurement of Fmeth. 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,u=Ffeed (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 inSvsin is applied.


## 7.2. Control ofFfeed

F
feed 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.Figure9 shows the simulated response with constant Ffeed. Table 4 shows performance measures.Table 4
Performance measures for three control strategies for controllingSvfa.
Measure
Const.Ffeed
Predictive
PI
IAE
d
1.52
0.20
0.090
GM
N/A
≈

∞
5.5
PM
N/A
≈

∞
71.0°Control ofSvfa with different control functions: constant Ffeed, predictive control (MPC), and PI control. Simulated responses in Svfa, Fmeth, Ffeed, and Svsin are shown.
(a)
(b)
(c)
(d)


## 7.3. Control ofSvfa

S
vfa is controlled to its setpoint, Svfasp=Svfamax⁡, using feedback from Svfaest 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(16)min⁡ufobj,
where
(17)fobj=∫tt+τh‍[e2(θ)+Cduu˙2(θ)]dθ+Che2(τh)
with constraint umin⁡≤u(t)≤umax⁡. The control error is e=Svfasp-Svfaest. Comparing with the criterion of predictive control of Fmeth, (11), the term e2(τh), which is e2 at the end of the prediction horizon, is now included. The term brings e(t+τh) approximately to zero. Without this term, e(t+τh) is 0.1 g/L, and the control signal is actually constant. It is found that Cdu=0.2 and Ch=20 are proper settings.It is found that the predictive control is considerably smoother withτh=4d than with τh=1 d which is used in Section 6. Increasing τs from 0.025 d, which is used in Section 6, to 0.1 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 areKc=50.9 (L/d)/(g VFA/L) and τi=0.9 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 Table2. The setpoint is Svfasp=Svfamax⁡=0.8 g/L. At t=10 d, the disturbance Svsin is changed as a ramp of slope 2 (g/L)/d during 1 d, which is the same variation as in Fmeth control (cf. Section 6). This is a reasonable variation for the real reactor. Measurement noise is not included in simulations.Figure9 shows simulated responses in Svfa, Fmeth, Ffeed, and Svsin with predictive control and PI control and with constant Ffeed. Table 4 shows performance measures.


### 7.3.4. Results and Discussion

Table4 shows performance and robustness measures with the three control strategies above. IAEd, defined by (13), is calculated over the simulated time interval. |e|max⁡ is the maximum control error. GM and PM are found as explained in Section 6.4.Comments
(i)
In Figure9 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 Svfa.
(ii)
The lower-right plot in Figure9 shows the manipulated Ffeed which is adjusted by the controllers. The control action is smoother with PI control than with predictive control.
(iii)
The performance measures shown in Table4 indicate that PI control of Svfa, based on feedback from UKF, is the best control strategy here.
(iv)
Also using a constantFfeed can be regarded as acceptable with the disturbance change simulated.
(v)
The upper-right plot in Figure9 illustrates that Fmeth is not under control. Although not shown here, Fmeth settles at steady state at approximately t=120 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 inSvsin. The relatively large estimate of Svsin is used in the prediction by the predictive controller, causing a large overshoot in Fmeth before it eventually reaches Fmethsp (plots are not shown here). This behavior is the same as with predictive control of Fmeth (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.3.1. Predictive Control

The optimization criterion of the predictive controller is selected as(16)min⁡ufobj,
where
(17)fobj=∫tt+τh‍[e2(θ)+Cduu˙2(θ)]dθ+Che2(τh)
with constraint umin⁡≤u(t)≤umax⁡. The control error is e=Svfasp-Svfaest. Comparing with the criterion of predictive control of Fmeth, (11), the term e2(τh), which is e2 at the end of the prediction horizon, is now included. The term brings e(t+τh) approximately to zero. Without this term, e(t+τh) is 0.1 g/L, and the control signal is actually constant. It is found that Cdu=0.2 and Ch=20 are proper settings.It is found that the predictive control is considerably smoother withτh=4d than with τh=1 d which is used in Section 6. Increasing τs from 0.025 d, which is used in Section 6, to 0.1 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 areKc=50.9 (L/d)/(g VFA/L) and τi=0.9 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 Table2. The setpoint is Svfasp=Svfamax⁡=0.8 g/L. At t=10 d, the disturbance Svsin is changed as a ramp of slope 2 (g/L)/d during 1 d, which is the same variation as in Fmeth control (cf. Section 6). This is a reasonable variation for the real reactor. Measurement noise is not included in simulations.Figure9 shows simulated responses in Svfa, Fmeth, Ffeed, and Svsin with predictive control and PI control and with constant Ffeed. Table 4 shows performance measures.


## 7.3.4. Results and Discussion

Table4 shows performance and robustness measures with the three control strategies above. IAEd, defined by (13), is calculated over the simulated time interval. |e|max⁡ is the maximum control error. GM and PM are found as explained in Section 6.4.Comments
(i)
In Figure9 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 Svfa.
(ii)
The lower-right plot in Figure9 shows the manipulated Ffeed which is adjusted by the controllers. The control action is smoother with PI control than with predictive control.
(iii)
The performance measures shown in Table4 indicate that PI control of Svfa, based on feedback from UKF, is the best control strategy here.
(iv)
Also using a constantFfeed can be regarded as acceptable with the disturbance change simulated.
(v)
The upper-right plot in Figure9 illustrates that Fmeth is not under control. Although not shown here, Fmeth settles at steady state at approximately t=120 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 inSvsin. The relatively large estimate of Svsin is used in the prediction by the predictive controller, causing a large overshoot in Fmeth before it eventually reaches Fmethsp (plots are not shown here). This behavior is the same as with predictive control of Fmeth (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 ofFmeth

### 7.4.1. Controllers

The third control strategy proposed in Section7.1 is controlling Fmeth to a setpoint, Fmethsp, set equal to the value of Fmeth at the ultimate operating point (cf. Table 2). Both predictive control based on feedback from UKF estimates and PI control based on measurement of Fmeth are simulated.


### 7.4.2. Simulations

The simulation scenario differs from the scenario of the simulations in Section7.3 as Svsin is now decreased instead of increased. Decreasing Svsin is selected here because, in the corresponding response, Svfa increases (in steady state), and an increase of Svfa is more critical than a decrease.In the predictive controller,τs is set as0.05 d, and τh is 1 d.Figure10 shows simulated responses in Svfa, Fmeth, Ffeed, and Svsin with predictive control based on feedback from UKF estimates and PI control based on measurement of Fmeth.Predictive control and PI control ofFmeth. The process perturbation is a change in Svsin. Simulated responses in Svfa, Fmeth, Ffeed, and Svsin are shown.
(a)
(b)
(c)
(d)


### 7.4.3. Results

(i)
As seen in Figure10, Fmeth is much closer to Fmethsp with predictive control than with PI control.
(ii)
With both predictive control and PI control,Svfa increases. Simulations over 400 d show that Svfa 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.4.1. Controllers

The third control strategy proposed in Section7.1 is controlling Fmeth to a setpoint, Fmethsp, set equal to the value of Fmeth at the ultimate operating point (cf. Table 2). Both predictive control based on feedback from UKF estimates and PI control based on measurement of Fmeth are simulated.


## 7.4.2. Simulations

The simulation scenario differs from the scenario of the simulations in Section7.3 as Svsin is now decreased instead of increased. Decreasing Svsin is selected here because, in the corresponding response, Svfa increases (in steady state), and an increase of Svfa is more critical than a decrease.In the predictive controller,τs is set as0.05 d, and τh is 1 d.Figure10 shows simulated responses in Svfa, Fmeth, Ffeed, and Svsin with predictive control based on feedback from UKF estimates and PI control based on measurement of Fmeth.Predictive control and PI control ofFmeth. The process perturbation is a change in Svsin. Simulated responses in Svfa, Fmeth, Ffeed, and Svsin are shown.
(a)
(b)
(c)
(d)


## 7.4.3. Results

(i)
As seen in Figure10, Fmeth is much closer to Fmethsp with predictive control than with PI control.
(ii)
With both predictive control and PI control,Svfa increases. Simulations over 400 d show that Svfa 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 Sections7.3, 7.2, and 7.4, it can be concluded that the best control strategy for safe reactor operation is controlling Svfa 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,Sbvs, Svfa, Xacid, Xmeth, and the assumed unknown organic content, Svsin, 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 Svfa 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 retainingFmeth 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, whereSvfa has a critical maximum value. This operating point is characterized by three corresponding values of Ffeed, Svfa, and Fmeth, 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 Svfa 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.



---

*Source: 572621-2014-04-22.xml*