On-off and PI Control of Methane Gas Production of a Pilot Anaerobic Digestion Reactor

A proposed feedback control system for methane flow control of a real pilot anaerobic digestion reactor fed with dairy waste is designed and analyzed using the modified Hill model, which has previously been adapted to the reactor. Conditions for safe operation of the reactor are found using steady-state responses of dynamic simulations, taking into account the upper limit of the volatile fatty acids (VFA) concentration recommended in the literature. The controllers used are standard process controllers, namely the on-off controller and the PI controller. Several PI controller tuning methods are evaluated using simulations. Two methods are favoured, namely the Skogestad method, which is an open loop method, and the Relaxed Ziegler-Nichols closed loop method. The two methods give approximately the same PI settings. Still, the Skogestad method is ranged first as it requires less tuning time, and because it is easier to change the PI settings at known changes in the process dynamics. Skogestad’s method is successfully applied to a PI control system for the real reactor. Using simulations, the critical operating point to be used for safe controller tuning is identified.


Introduction
This paper attempts to answer the following questions related to a real pilot upflow anaerboic sludge blanket (UASB) reactor fed with dairy waste: What are the benefits and drawbacks of feedback, or closed loop, control of the produced methane flow compared to using open loop control, i.e. a constant feed rate?Assuming the use of standard process controllers, namely on-off and proportional-integral (PI) control, how do the control systems perform?How should the PI controller be tuned?Some of these questions are addressed using both simulations and practical experiments, while some are addressed only using simulations.

The pilot plant
The reactor is a 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 waste diluted with 25% water and filtered with a sieve, and outputs are fertilizer and biogas consisting of 70-75% methane.The reactor temperature is kept fixed at its setpoint with an automatic temperature control system, Haugen et al. (2013b).A description of the plant, including its monitoring and control system, is provided in Haugen et al. (2013c).

Anaerobic digestion (AD) of animal wastes
AD of animal wastes can produce biogas with methane to be used as an energy source, and a liquid effluent containing valuable nutrients.Moreover, AD reduces methane emission, odours and contaminants.AD bioreactors are effective as they allow for relatively high load rates (feed rates) and small reactor volumes.Various theoretical and practical aspects of AD processes are described e.g. in Tchobanoglous et al. (2003) and Deublein and Steinhauser (2010).A presentation of AD of animal wastes, from dairy, beef, poultry, and swine, is provided e.g. in Husain (1998).
Literature review of reactor control Bernard et al. (2001b) have implemented a modelbased 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 avoid VFA accumulation despite organic load disturbances.The so-called AM2 model, Bernard et al. (2001a), is used for design and simulation.Puñal et al. (2003) 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. (2005) have designed a modelbased controller for maintaining the COD (chemical oxygen demand) of the reactor effluent at its setpoint, using the AM2 model, Bernard et al. (2001a).Méndez-Acosta et al. (2010) 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.
In neither of the above control systems, the biogas flow is controlled.Focus is on reactor stability rather than on energy production (1 Nm 3 methane at NTP = 9.95 kWh).In the papers referred below, the biogas flow is controlled.Strömberg (2010) has 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 are based on the ADM1 model, Batstone et al. (2002).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.A short description of these controllers follows.
The extremum-seeking variable gain controller by Liu et al. (2006) has the structure of a cascade control system where the primary loop implements biogas flow control, and the secondary loop implements pH control to stabilize pH.
In the disturbance monitoring controller by Steyer et al. (1999), disturbances in the form of pulses are added to the feed rate, and from the gas flow response the feed rate is adjusted to obtain maximum production.Pulsing is stopped if measured pH is below a critical value.
In the hydrogen-based variable gain controller by Rodriguez et al. (2006), online measurements of methane and hydrogen concentrations of the biogas are measured online and are used by the controller for approaching a preset maximum methane gas flow.The controller is based on a relation between hydrogen concentration and effluent COD (Chemical Oxygen Demand) concentration as found from the ADM1 model, Batstone et al. (2002).Strömberg et al. (2013) note that no uniform tuning method could be derived to tune the three controllers.Instead, the controllers were tuned by first running a large number of simulations to become familiar with the controller performances, and then the parameters were tuned manually.

Paper outline
This paper is organized as follows.Section 2 gives a short description of the proposed methane gas flow control system.Section 3 covers on-off control.Section 4 covers PI control, including an analysis of the control system robustness against process parameter changes.A discussion is given in Section 5, and conclusions are given in Section 6. Appendix A describes mathematical models used, and abbreviations and nomenclature are given in Appendix B.

Computing tools
MATLAB and SIMULINK (MathWorks, Inc.) are used for numerical computations and simulations.The algorithms of the real control system are implemented in LabVIEW (National Instruments, Inc.) running on a laptop PC.

Control system objective
As is clear from the control systems referred in Section 1, there are alternative objectives for reactor control, i.e. obtaining a specified VFA concentration, or alkalinity, or obtaining a specified biogas production.
For the present reactor, the control objective is proposed as follows: F meth is maintained at its setpoint, F methsp , assuming safe operation conditions, defined below.
Generally, the specific value of F methsp may be calculated as the solution of a model-based optimization problem with a proper optimization criterion, e.g.maximum gas production, or economic optimization where power loss, energy prices, and value of money are taken into account.However, formulation of the optimization problem is not discussed further in the present paper, but in a forthcoming paper.

Selection of control variable
F meth is the process output variable to be controlled.An obvious candidate control variable is F feed as it has a clear impact on F meth .Also T reac , which in the present reactor is controlled to its setpoint using a feedback control system, has a clear impact on F meth , cf.Haugen et al. (2013b).Therefore, T reac is a candidate as control variable.However, as pointed out in Tchobanoglous et al. (2003), T reac should be kept ideally constant, variations within ±0.5 o C being acceptable, to avoid stressing the methane-generating microorganisms (methanogens).Therefore, T reac is not considered a usable control variable.

Implementation of the control system
In the methane gas flow control systems studied in the present paper, F feed is the control variable while F meth is the process output variable.The AD reactor with the methane flow control system is illustrated in Figure 1.The control signal from the gas flow controller, FC, acts on the feed pump (actuator) which is assumed to provide the feed flow, F feed , demanded by the controller.The controller adjusts F feed based on the control error which is the difference between the F meth measurement and its setpoint.This measurement is provided by sensor FT.In practice, 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 in filters which are described in Appendix A.3.
The gas flow controller manipulates the peristaltic feed pump using PWM (Pulse-Width Modulation).It is found that proper PWM settings are as follows: Fixed cycle time of 700 sec, on-value of control signal corresponding to 714 L/d, and off-value corresponding to zero L/d.Two factual benefits of using PWM control compared with analog control are (1) the calibration of the pump is needed only at the on-state of the flow rate, and (2) blockings in the feed pipeline are reduced.

Control functions
The controllers reviewed in Section 1 can be regarded as non-standard process controllers (though it can be claimed that a fuzzy logic controller is a standard controller).For the present reactor, it is proposed to use standard process controllers as on-off controllers and PI controllers, which are relatively simple controllers.For the latter there are many tuning procedures available, but there is no guarantee that any method gives successful PI settings.In this paper, several tuning methods are tested to identify the most suitable methods.In general, feedforward control can be a very ef-ficient control method to compensate for severe disturbances, assuming that these disturbances are measured continuously or are estimated continuously using a soft-sensor, e.g. a Kalman Filter algorithm, Dochain (2008), Simon (2006).For the present pilot reactor, variations in the volatile solids (VS) concentration of the feed, S vsin [g/L], are regarded as the most important disturbance acting on F meth .Figure 2 shows a plot of S vsin over a period of one year.The largest change between the samples occurs around t = 150 d, and is approximately ∆S bvsin = 6 g/L.
A simulation is here used to indicate the response in F meth due to a change in S bvsin .It is assumed that the change is a step of amplitude 6 g/L.F meth is controlled with a PI controller tuned with the Skogestad method, cf.Section 4.4.1.Figure 3   The maximum transient offset of F meth from the setpoint is 1.9 L CH 4 /d, which volume-normalized is 1.9 L CH 4 /d/250 L = 0.0075 L CH 4 /d/L.Assuming this transient offset is acceptable in a practical application, feedforward control is not needed.
Assuming the offset due to the variations in S vsin is problematic, there is still a practical problem about implementing feedforward control because of the lack of appropriate online sensors.A soft-sensor for S vsin is an alternative, but is here regarded an advanced method, and is therefore not discussed in the present paper.Using a soft-sensor for S vsin is relevant for model-based control, to be addressed in a forthcoming paper.
2.5.Safe reactor operation and attainable operating points Hill et al. (1987) have found, from a comprehensive study of literature reporting operational data for reactors fed with swine and beef waste, and confirmed by their own laboratory experiments, that S vfa = 0.8 g/L is a good indicator of the reactor health.S vfa ≥ 0.8 g/L indicates an impending reactor failure, i.e. a reduction of methane production, while S vfa ≤ 0.8 g/L indicates that the reactor is operated successfully, i.e. that the reactor is healthy.Hill et al. found that also the proprionic to acetic acid (P/A) ratio is a good indicator.However, this ratio can not be calculated from the mathematical model used in this paper, and therefore, the analysis here is not based on this ratio.Hill et al. did not use dairy waste in their analysis since reliable data for such waste were not available.Nevertheless, it is here assumed that the aformentioned limit applies approximately also for reactors fed with dairy waste.A support for this assumption is that the validated AD reactor model by Hill (1983) has the same parameters describing the AD process for dairy, swine, poultry, and beef waste, except for parameters expressing the fraction of the organic feed that is degradable, but the AD process dynamics are independent of the latter parameters.
There are several alternative AD process indicators and control parameters.Angelidaki et al. (1993) identify e.g.ammonia as an important parameter for AD process control, particularly for animal waste rich on ammonia.Bernard et al. (2001b) reports that the internal alkalinity to total alkalinity ratio is an important indicator, and control parameter.According to Tchobanoglous et al. (2003), a pH level lower than 6.8 is inhibitory on the methanogenesis.The importance of monitoring and controlling these, and other, parameters, depends on the type of fed substrate, e.g.food waste, industrial waste, etc.However, from the literature, the ammonia content is relatively small and therefore hardly inhibitory for dairy waste (but it is for swine and poultry waste), the alkalinity is relatively large and not subject to large variations, and the buffer capacity is high implying that a proper pH level is maintained.Thus, VFA remains the main AD process parameter to be used here.
It is assumed that the reactor is represented by the  Figure 4 shows simulated static (steady-state) responses to a range of constant feed rates (F feed ).T reac is 35 o C which is typical for AD reactors.The global, unconstrained maximum of F meths is found as 214 L CH 4 /d which is obtained with F feed = 83.4L/d, represented with right green vertical lines in Figure 4.However, at this operating point, X acids is virtually zero, which is coherent with the middle left plot showing that S bvsin is not being degraded.Therefore, this unconstrained maximum is not regarded as a viable operating point.
In Figure 4, the cyan horizontal line in the S vfa plot represents S vfa = 0.8 g/L.At this value, F feed = 35.3L/d, which is represented by red vertical lines in the plots.At this feed flow, F meths = 174 L CH 4 /d which is then the maximum attainable F meths under safe con-ditions.The corresponding hydraulic retention time (HRT) of the reactor is 174/35.3= 4.9 d.
Assume that the controller output range, i.e. the range of F feed , is restricted by the user to ensure safe operating conditions.To continue the above example, assume that the upper limit of F feed is set to 35.3 L/d which, according to the model, corresponds to F meths = 174 L CH 4 /d.The setpoint is set to F methsp = 174 L CH 4 /d.Assume that for the practical reactor, the factual F meths that is obtained with F feed = 35.3L/d, is less than 174 L CH 4 /d.Then, obviously, the steadystate control error is non-zero.To obtain zero steadystate error, either the upper limit of F feed should be set higher than 35.3 L/d, or F methsp should be set smaller than 174 L CH 4 /d, of which the latter alternative is the safest.

Comparing feedback control with open loop control
To demonstrate the effect of F meth control, Figure 5 shows experimental time-series of F meth and F feed , and T reac , with (automatic) control and without control.
It is clearly demonstrated that F meth drifts less with control than without control.F meth remains close to F methsp even after the setpoint is changed.In the case of feedback control, F feed is of course varying, while it is constant in open loop control.T reac is actually different in the two cases, but it is assumed the difference between the two cases is independent of the temperature difference.If the drift of F meth is acceptable when a constant F feed is used, feedback control may be superfluous.

On-off control
The on-off controller can be regarded as the simplest feeedback controller available, Johnson (2000).It can be used without any tuning, except deciding the onvalue and the off-value of the controller output as well as the threshold level, d e .With on-off control, the control system oscillates, and typically, the offset of the mean value of the control error from the setpoint is non-zero.
The on-off controller function can be defined as where e is the control error, and d e is an adjustable dead-band to avoid switching of u due to (measurement) noise in e. u on and u off are constant control signal levels.In applications with the pilot reactor, d e = 0.5 L CH 4 /d.Table 1 summarizes a number of settings and absolute and volume-normalized characteristics found from simulated and experiments on the real reactor further described in the following.

Simulations
In simulations, though not displayed here, F methsp (t) has the typical form of a sinusoidal oscillation, while u(t) is a square wave.

Responses of the real reactor
Figure 6 shows F meth from an experiment with on-off control on the real reactor.
Comments and conclusions • For both the simulated and real reactor the mean control error are nonzero, but probably acceptable for both the simulated and the real reactor.• The on-off behaviour of the control signal, and hence the feed flow, may be acceptable in practical applications.Actually, it is observed that for the present reactor, discontinuous feeding reduces the frequency of the blockings.

Controller function
PID control is prevalent in industrial applications.The PID controller provides smooth control as opposed to on-off control, but its parameters must be tuned to fit the dynamics of the process to be controlled.The applied PID controller is based on Euler backward discretization of the following continuous time PID controller, with time-step τ s = 2 s: Typically, the derivative term may provide control stability and agility, but is nevertheless often deactivated, i.e. τ d = 0, in practical applications because of its problematic propagation of measurement noise causing a noisy control signal.For the pilot reactor, it is decided to not use the derivative term.

Selection of controller tuning methods
It is of interest to compare different PI(D) tuning methods to arrive at a conclusion about recommended methods.The following methods are applied to the simulated reactor, and some of them are applied to the real reactor: Among open loop methods, the SIMC method (Simple IMC), Skogestad (2004), here denoted the Skogestad method, is selected.2004), but it is not clear whether these methods have important benefits compared with the Skogestad method.Among closed loop tuning methods, the famous Ziegler Nichols (ZN) closed loop method is applied, although it is expected to give small stability margins.As an alternative, the Relaxed Ziegler Nichols (R-ZN) method proposed by Haugen and Lie ( 2013) is tested.The Tyreus and Luyben (TL) method, Tyreus and Luyben (1992), is probably the best known method to modify the ZN closed loop PI settings to obtain more relaxed control.However, the R-ZN method compares favourably with the TL method, cf.Haugen et al., (2013d).Therefore, the TL method is not tested here.

Summary of results
Table 2 gives a summary of results for a simulated reactor based on the model presented in Appendix A, and for the real reactor.The table shows controller settings, the gain margin (GM), the phase margin (PM), and the closed-loop response-time, τ r [d], which is estimated as  13.5 0.54 2.8 = 8.9 [dB] 37.4 0.25 • The values of K c for the simulated reactor are roughly half of the values for the real reactor, while the values of τ i differ little.The variation in K c is of course due to modeling errors.Actually, the modified Hill model was adapted to various timeseries from time intervals relevant to controller tuning, but with no significant changes in the pertinent model parameters, nor in K c .From this it can be concluded that using a phenomenological mathematical model as the basis for tuning a controller for the real reactor, is dubious.However, for the present reactor, the difference in values of K c is safe.
• The PI settings calculated from the estimated process transfer function have approximately double the value of K c and half the value of τ i comparing with the respective settings found by the Skogestads method and the R-ZN method applied directly to the real reactor, and thus, the former settings are expected to give more aggressive control.Due to practical obstacles, the PI settings found from the estimated process transfer function have not been applied to the real reactor, so it is not known if they are applicable on the real reactor, or not.

Applications of controller tunings 4.4.1. Skogestad tuning
Simulations and real experiments indicate that the reactor dynamics can be characterized approximately as "time-constant with time-delay" with a time-constant of a few days and a time-delay of a few hours.As pointed out in Haugen and Lie (2013), the Skogestad method, with a proposed justified reduction of the integral time setting, gives the following PI settings for processes where the time-delay is less than one eighth of the time-constant, which is the case for the pilot reactor: where K ip , the integrator gain, and τ delay , the (apparent) time-delay, can be found from a simple process step response.K ip can be calculated as  2. The stability margins are within the acceptable limits given by ineqs.( 5)-( 6).
Figure 3 shows simulated responses of the control system due to a setpoint step and a disturbance step in S vsin .Figure 9 shows responses on the real reactor with the above PI settings.There is no indication of poor stability.The relatively large response in F meth between t = 446.0 and 446.5 d is assumed being due to a disturbance.

Comments and conclusions
• The Skogestad PI tuning method, here used as an open loop step response method, gives good results.
• The step response test can be accomplished within approximately half a day which is considerably shorter than e.g.relay-based tuning methods which may require more than two days, cf.Section 4.4.2.

Ziegler-Nichols PI tuning based on relay oscillations
Åstrøm and Hägglund (1995) suggested a relay or onoff controller to replace the P controller in the tuning phase of the ZN closed loop method, thereby avoiding the trial-and-error procedure since the oscillations come automatically.During the relay tuning the control signal becomes a square wave.Assuming the oscillation in the process output is (approximately) sinusoidal, as is the case with the present reactor, the  ultimate gain, K cu , is calculated as where A u is the amplitude of the on-off control signal, and A e is the amplitude of the control error and the process output.For a PI controller, the ZN settings are K c = 0.45K cu and τ i = P u /1.2 where P u is the period of the oscillation.

Results
Using simulations with on-off controller (not shown) as the basis for relay tuning, gives A u = (45 − 5)/2 = 20 L/d, A e = 3.3 L M , P u = 1.16 d.The resulting PI settings are shown in Table 2 together with control system characteristics.The phase margin (PM) is 25.6 o which is less than the lower limit in eq. ( 6).Simulations, not shown here, confirm relatively small stability margins as responses are oscillatory.It was decided not to use the Ziegler-Nichols method with the real reactor since the theoretical results are not satisfactory.

Relaxed Ziegler-Nichols PI tuning
The Relaxed Ziegler-Nichols (R-ZN) PI tuning method is proposed by Haugen and Lie (2013) to give more relaxed control, i.e. improved stability, compared with the ZN closed loop method.The PI settings are: where K cu and P u can be found from relay oscillations.
The tuning parameter k r can be used for enhanced relaxation.Simulations with the modified Hill model, Haugen et al. (2013a), indicate that the default value k r = 1 works well with the present pilot reactor.

Simulations
PI settings are calculated from simulated relay oscillations, not shown here.Table 2 shows the resulting PI settings, and control system characteristics.The stability margins have acceptable values.Since these PI settings differ little from those found with the Skogestad method in Section 4.4.1, simulations with R-ZN settings are not shown here (Figure 3 shows responses with Skogestad PI settings).

Practical results
The relay oscillations shown in Figure 6 2. Since these settings differ little from those found with the Skogestad method, cf.Section 4.4.1, it was decided not to perform separate experiments on the real reactor with R-ZN settings.

Comments and conclusions
• Both simulations and real experiments indicate successful controller tuning using the R-ZN method.The PI settings become close to those obtained with the Skogestad method, which is not a surprise since this method is designed from a combination of the Skogestad method and the ZN method.
• The Skogestad method is here favoured compared with the R-ZN method, due to the following observations: Firstly, the Skogestad method has a shorter tuning phase, namely approximately 0.7 d, cf. Figure 8, while the tuning phase of the R-ZN method is 2-3 days.Secondly, in the Skogestad method, retuning the controller in the case of a known process parameter change, e.g. an increase of the apparent τ delay due to an increased filter time-constant, can be accomplished without performing any new experiment.With the R-ZN method, a new experiment is needed.

Optimal PI tuning based on the modified Hill model
In this method, PI controller parameter vector where e is the control error, u = Ḟfeed is the rate of change of the control signal, and R is a user-selected cost coefficient.With R = 0, eq. ( 14) is identical to the well-known IAE index, Seborg et al. (2004).The larger R, the more cost of control signal variations, and smoother, but also slower, control actions can be expected.Having only one parameter, R, to be tuned is a much easier tuning problem than having two parameters, K c and τ i .Furthermore, R has an intuitive interpretation.Many alternative objective functions are possible, e.g.quadratic functions instead of absolute values, frequency response-based functions, etc. Equation ( 14) is here selected since it is an enhancement of the IAE index.A quadratic objective function was tested, but no benefits were identified compared to the selected function.
In eq. ( 13), C is a constraint on the stability margins in terms of |S(jω)| max , where S(s) = 1/[1 + L(s)], where L(s) = H c (s)H ADm (s) is the loop transfer function.H ADm (s) is given by eq. ( 25).The acceptable range of |S(jω)| max is set as [1.2, 2.0], according to Seborg et al. (2004).However, this constraint was not active at the optimal solution (its value was 1.87).
f obj (p c ) is calculated from simulations with the (nonlinear) modified Hill model, Haugen et al. (2013a).The simulator is based on numerical integration of the differential equations using the Euler explicit numerical method implemented in native for-loops in a MATLAB script. 1  The optimization problem is here solved using "brute force" (BF), i.e. f obj (p c ) is calculated over a grid of equidistant values of K c and τ i defined in respective arrays (MATLAB), and the optimal p c is found by searching the matrix of stored values of f obj for the minimum.This gives a global, approximate solution.
If a more precise value is desired, either BF optimization can be repeated but with the new grid cells covering the original grid cells containing the global optimum candidate, or a local optimizer, Edgar et al. (2001), can be applied with the global optimum candidate as the initial guess. 2Both these alternatives were tested, with approximately the same optimum, but the repeated BF method being considerably easier to implement.

Application to simulated reactor
In eq. ( 14), t 1 = 0 d, and t 2 = 5 d.The reactor is initially in steady-state.F methsp is constant (88 L/d).At t = 1 d, S vsin is increased as a step of amplitude 2 g/L.The arrays of K c and τ i are equidistant with 100 elements each.By trial-and-error, a proper value of R is found as 0.3.Table 2 shows the resulting optimal PI settings.The stability margins and responsetime as calculated from the linear model presented in Appendix A.3.The stability margins have acceptable values.Figure 10 shows a simulation with the optimal PI settings, indicating acceptable stability.
1 Comparing with implementation of the simulator in SIMULINK, the computational time is reduced by a factor of about 100 with for-loops. 2One example of a local optimizer is MATLAB's fmincon function.

Comments
• The optimal PI settings do not differ much from those found with the Skogestad method and the R-ZN method.
• Optimal tuning has not been applied to the real reactor.However, the optimal PI settings found from simulations will probably work well on the real reactor since K c is smaller than K c found from the Skogestad method applied to the real reactor, and the values of τ i do not differ much.
• Optimal tuning is a flexible tuning method since it allows for alternative types of models and alternative objective functions.

PI tuning using estimated transfer function
Figure 11 shows real F meth and F feed , and simulated F meth using the real F feed , over time interval of 4 d.The simulation is based on transfer function H p (s) estimated from the shown real F meth and F feed .H p (s) is estimated using the n4sid function in MATLAB with automatic detection of the best model order3 , and using the delayest function to estimate the time-delay Models estimated from other time-series do not differ much from eq. ( 15).Among the large number of controller tuning methods which can be used for eq.( 15), the Skogestad method is selected.Since τ de is considerably smaller than τ e , the Skogestad PI setting for "integrator with time-delay" processes given in Section 4.4.1 with K pi e = K e /τ e = 0.27 [(L CH 4 /d)/(L/d)]/d can be applied.The resulting PI settings are given in Table 2.The PI control system with eq. ( 15) as controlled process has stability margins of GM = 2.8 = 8.9 dB and PM = 37.4 o , which are within the acceptable ranges ineqs.( 5)-( 6).K c and τ i found here are, respectively, larger and smaller compared with the values found using the Skogestad method in Section 4.4.1 applied directly to the real reactor, cf. 2. This is due to τ de = 0.135 d being smaller than the time-delay of 0.

Conclusions about PI tuning method
From the results in the above sections, the Skogestad method is favoured among the various tuning methods due to the following benefits.The step response experiment is simple, and the experimental period is short which is an important benefit with slow processes such as bioreactors.Known changes in the (apparent) time-delay can be accounted for in the PI(D) settings without new experiments.The control agility can easily be adjusted via the closed loop time-constant.The method can be applied without any prior mathematical model.Finally, the method has proven, in the present and in other application, to give good tuning results.
The R-ZN PI tuning method also works well, and can be expected to give tuning results similar to the Skogestad method.However, accomplishing the former method may take three or more times longer time compared with the Skogestad method.Furthermore, known changes in the (apparent) time-delay are not easily accounted for, without a new tuning.

Control system robustness against
process parameter changes 4.5.1.Introduction Figure 12 shows the static (steady-state) F meth , here denoted F meths , as a function of constant F feed for three different T reac found by simulations with the modified Hill model, Haugen et al. (2013a).The static process gain is defined as K p is the slope of the curve in Figure 12.Depending on F feed , K p is positive, zero or negative.In Figure 12 it can be seen that for a given T reac , there is a maximum achievable F meths , which defines the feasible setpoint of F meth .

Dependency on T reac
To illustrate temperature dependency, the operating point value is set to   assumed.Assume that the controller will have fixed settings.If a controller is tuned at the critical operating point, the control system will remain stable in any other operating point.On the other hand, if the controller is tuned at a non-critical operating point, the control system may become unstable at other operating points.Some alternatives to using fixed controller settings to handle varying process dynamics are: • Continuous adaptive tuning based on continuous estimation of a transfer function model, Åstrøm and Wittenmark (1994).
Implementation of experimental Gain scheduling is straightforward.T reac and F feed may be used as input variables to the table, and the PI settings are the output variables.Each of the PI settings can be found experimentally using e.g. the Skogestad method.However, implemention of the above mentioned alternatives are not designed nor analyzed here.
For the analysis of the control system stability the transfer function model of the methane flow control system presented in Appendix A.3 is used.The analysis is accomplished as follows: A PI controller is tuned using the R-ZN method (with tuning parameter k r = 1) at one specific operating point which is here denoted the basic operating point: (F feed = 10 L/d, T reac = 25 o C).In this tuning method the ultimate gain, K cu , and the ultimate period, P u , are needed to calculate PI settings K c and τ i .In the tuning, K cu = 1 and τ i = ∞ initially.Then, K cu = GM and P u = 2π/ω 180 where ω 180 [rad/d] is the gain margin crossover frequency.The PI settings are used at different operating points, and the stability margins, namely the gain margin, GM, and the phase margin, PM, are calculated from the transfer function model.
From the results, a conclusion is made in the following about what is the critical operating point for controller tuning, and an attempt is made to explain the results using the (nonlinear) mathematical reactor model.
Table 3 summarizes the results.The fixed PI settings found at the basic operating point are K c = 1.29 (L CH 4 /d)/(L/d), and τ i = 0.91 d.The upper left cell of Table 3 represents the basic operating point.In Table 3 the following observations are made: 1.The stability margins decrease with increasing T reac .
2. Except at T reac = 35 o C, where the stability margins are almost independent of feed , the stability margins decrease with decreasing F feed .
From these observations, the following general guideline is proposed, at least for T reac less than 35 o C: The critical operating point regarding controller tuning is maximum T reac and minimum F feed .A PI controller with fixed tuning should be tuned in this operating point.(At T reac = 35 o C, the controller tuning seems to become independent of F feed .) Below is an attempt to explain the above two observations using the mathematical reactor model.
1. Regarding observation 1: Figure 12 shows the steady-state F meth as a function of constant F feed for different T reac .The static process gain is defined with eq. ( 16).At least for F feed = 10 L/d and 25 L/d which are assumed above, K p increases with T reac , and hence, a reduction of control system stability margins can be expected.
An alternative explanation may be found directly from the modified Hill model: From model eqs.( 25), ( 27), (28) in Haugen et al. (2013a), it is seen that F meth becomes more sensitive to S vfa as the temperature increases, hence the process gain increases.This increased sensitivity may explain the reduced stability (margins) as T reac is increased.
2. Regarding observation 2: From Figure 12, K p becomes larger, and hence the control system stability margins are reduced, if F feed is reduced.

Discussion
Ideally, all the questions stated in the Introduction should be addressed with practical experiments.However, this was not possible for practical reasons, so some questions are addressed using simulations only.Since the modified Hill model has shown to represent the pilot AD reactor well, Haugen et al. (2013a), it is assumed that the results obtained from simulations hold qualitatively, and to some extent, quantitatively.Both on-off control and PI control are found being successful for controlling the methane gas flow, on a simulated reactor as well as on the practical reactor.For PI controller tuning, the Skogestad method, which is an open loop tuning method, is identified as the favoured tuning method.Also the R-ZN method, which is a closed loop method based on relay oscillations, works well.It is believed that the identification of these tuning methods can reduce time and efforts in controller tuning for AD processes.
Dairy waste as AD feedstock has large buffering capacity, and its composition is relatively constant.If the feedstock is more complex, as with poultry and swine waste and food waste, a richer mathematical model able to predict other AD variables than those of the Hill model, e.g.pH, alkalinity, partical alkalinity (PA), pH, ammonia, and carbon dioxide may be useful.Two model candidates are the AM2 model by Bernard et al. (2001a) and the ADM1 model, Batstone et al. (2002).Overviews over AD models are given in e.g.Gavala et al. (2003) and Strömberg (2010).

Conclusions
Using a mathematical model of the AD reactor and the specific upper limit of the concentration of volatile fatty acids, known from the literature, safe operating conditions for the reactor can be found.These conditions imply an upper limit of the feed rate, and an upper limit of the gas flow setpoint.These limits are theoretical, and should be adjusted on a practical reactor to avoid non-zero steady state control error.
For the present pilot reactor, both simulatons and practical experiments indicate that on-off control is a viable feedback controller if the oscillation in the feed rate and biogas flow can be tolerated.If smooth control is important, PI control is appropriate.The Skogestad method is favoured as a PI controller tuning method since it is easy to apply and gives good tuning results with the present reactor.Also, the R-ZN closed loop tuning method works well, but the time needed to accomplish the tuning is longer than with the Skogestad method.
Simulations indicate that it is safe for control loop stability to tune a PI controller with fixed parameters at low feed flow.

A.1. The modified Hill model
In this paper, simulations are based on the modified Hill model adapted to the pilot AD reactor.The model is derived in Haugen et al. (2013a).where x is the state vector:

A.3. Transfer function model
The transfer function from from ∆F feed to ∆F meth can be calculated from the linear state space model, eqs.( 17) and ( 18 In eq. ( 20), ∆F meth represents the (deviation of) the "raw" CH 4 gas flow.In practice, the CH 4 gas flow is known from the multiplication of the biogas flow online measurement and the CH 4 concentration online measurement.These measurements are smoothed with lowpass filters with the following respective timeconstants: • Time-constant of main biogas flow measurement filter: τ f1 = 0.2 d (21) • Time-constant of additional biogas flow measurement filter: Taking the above mentioned dynamic elements into account, the following transfer function from the feed flow to the methane gas flow measurement is obtained: where, for i = 1, 2, and 3,

B.2. Nomenclature
In the paper, but not shown in the list below, subindex "s" is used to represent "steady-state" or "static".The list below contains only symbols which are used in this paper.A complete list of symbols for the modified Hill model is in Haugen et al. (2013a).
A e : Amplitude of the control error and the process output.
A u : Amplitude of the on-off control signal.
c: Factor used to define lower limits of biodegraders X acid and X meth .
f obj : Objective function.ω c [rad/d]: Amplitude crossover frequency, which may be defined as the control system bandwidth.X acid [g acidogens/L]: Concentration of acidogens.

Figure 2 :
Figure 2: Plot of S vsin from laboratory analysis over one year.Mean value: µ Svs in = 29.7 g/L.Standard deviation: σ Svs in = 2.0 g/L.
shows simulated responses of the control system due to a step change in S vsin .(The response due to a step change in the setpoint of F meth at t = 5 d, is relevant in Section 4.4.1.)

Figure 3 :
Figure 3: Simulated responses of the control system due to a setpoint step and a disturbance step in S vsin .The PI controller is tuned with the Skogestad method.

Figure 4 :
Figure 4: Simulated static (steady-state) values of a number of variables versus F feed (constant) at T reac = 35 o C. Vertical and horizontal lines are explained in the text.

Figure 6 :
Figure 6: Experimental F meth with on-off methane gas flow control.
) τ r is approximately the time-constant of the control system.The above frequency response characteristics are based on the transfer functions model described in Appendix A.3, except for the tuning method denoted "Skogestad with estimated transfer function" where the frequency response characteristics are based on the estimated transfer function.Seborg et al. (2004) present the following ranges of the stability margins: 1.7 = 4.6 dB ≤ GM ≤ 4.0 = 12.0 dB (5) and 30 o ≤ PM ≤ 45 o (6)Comments to Table2: where ∆F feed [L/d] is the applied step amplitude, and S [(L CH 4 /d)/d] is slope of the time-delayed rampformed response in F meth .Simulations K ip and τ delay are found from an open-loop step response where F feed is changed as a step of amplitude ∆F feed = 1 L/d.From the step response shown in Figure 7, where the red line is the "time-delayed ramp" step response of the assumed integrator with timedelay, S = 0.883 (L CH 4 /d)/d and τ delay = 0.23 d.This gives K pi = S/∆F meth = 6.6/20 = 0.33 [(L CH 4 /d)/(L/d)]/d.The resulting PI settings and frequency response characteristics calculated from the linear model presented in Appendix A.3 are presented in Table

Figure 8
Figure 8 shows the response in F meth due to a step in F feed from 39 L/d to 19 L/d, i.e. the step amplitude is ∆F feed = −20 L/d.A negative step is used because the initial value of F feed is relatively large.From the step response, τ delay = 0.31 d and S = −6.6 (L CH 4 /d)/d are estimated.The resulting PI settings are shown in Table2.

Figure 7 :
Figure 7: Simulated open-loop step response where F feed is changed as a step of amplitude 1 L/d.

Figure 8 :
Figure 8: Open-loop step response for the real reactor.

Figure 9 :
Figure 9: Responses on the real reactor with PI methane flow controller tuned with Skogestad's method.
are used as the basis for the PI controller tuning.The and off values of the controller are F feedon = 45 L/d and F feed off = 5 L/d, respectively, giving A u = (45 − 5)/2 = 20 L/d.From the oscillations, A e = 2.0 L/d and P u = 1.1 d.The resulting PI settings are shown in Table

Figure 10 :
Figure10: Simulation of control system with optimal PI settings.

Figure 11 :
Figure 11: Real and simulated F meth .

Figure 12 :
Figure 12: Steady-state F meth as a function of constant F feed for different T reac .

F
meth due to a unit step in F feed [L/d].Blue: T reac =25 C. Red: T reac =35 C.

Figure 13 :∆Figure 14 :
Figure 13: ∆F meth (t) due to a step change of F feed of amplitude 1 L/d at T reac = 25 o C and 35 o C.

For
linear analysis and design, e.g.frequency response based method, a local linear version of the modified Hill model is used.The model has the standard form: ∆ ẋ = A∆x + B∆F feed (17) ∆F meth = C∆x + D∆F feed (18) ), with∆F meth (s) ∆F feed (s) ≡ H AD (s) = C (sI − A) −1 B + D (20)where matrices A, B, C, and D are given in Section A.2.
is observed that there is a time delay in the observed measured responses in the CH 4 gas flow which is approximately τ d = 0.05 d (24)

F
feed [L/d]: Influent or feed flow or load rate, assumed equal to effluent flow (constant volume).F meth [L CH 4 /d]: Methane gas flow.GM: Gain margin.k r : Parameter of the Relaxed Ziegler-Nichols tuning method.K p [(L CH 4 /d)/(L/d)] = ∂F meths /∂F feed : Static process gain.

P
u [d] Period of oscillation.PM [degrees]: Phase margin.S vfa [g VFA/L]: Concentration of VFA acids in reactor.S vfain [g VFA/L]: Concentration of VFA in biodegradable part of influent.S bvs [g BVS/L]: Concentration of BVS in reactor.S bvsin [g BVS/L]: Concentration of BVS in influent.S vsin [g VS/L]: Concentration of volatile solids in influent.T reac [ • C]: Reactor temperature.τ i [d]: Controller integral time.τ d [d]: Controller derivative time.V [L]: Effective reactor volume.

Table 1 :
Settings and steady-state characteristics for simulated and practical experiments with onoff control.

Table 2 :
Results with various PI tuning methods for simulated and for real reactor.

Table 2 .
and the response in ∆F meth (t), which is the deviation from the operating point due to a step change of F feed of amplitude ∆F feed = 1 L/d, is simulated for T reac = 25 o C and 35 o C. Figure 13 shows ∆F meth (t) simulated with the linearized model presented in Appendix A.2.The simulations indicate that the dynamics of the reactor is faster the larger T reac .To quantify the agility of the dynamics, the 63% response time, resembling the timeconstant, is approximately 7.5 d with T reac = 25 o C and approximately 4.1 d with T reac = 35 o C.

Table 3 :
Stability margins of the methane control system at various operating points.The upper left cell is the basic operating point.