Temperature Control of a Pilot Anaerobic Digestion Reactor

Results of analysis and design and implementation of a temperature control system for a practical pilot anaerobic digestion (AD) bioreactor fed with dairy waste are presented. A dynamic model of the reactor temperature is used as the basis for theoretical results, including simulations. Controller functions include on-off control, proportional plus integral (PI) control, and feedforward control. Various PI controller tuning methods are compared. The need for adaptivity of PI settings is investigated. Results for a simulated full-scale reactor are given.


Introduction
The main aim of this paper is to present the results of analysis and design and implementation of a temperature control system for a practical pilot anaerobic digestion (AD) bioreactor fed with dairy waste.The bioreactor is a part of the biological plant for nutrient and energy recovery at Foss Biolab, Skien, Norway.A description of the plant and the monitoring and control system can be found in Haugen et al. (2013a) and Haugen et al. (2013b).Results of analysis and design of a temperature control system for a hypothetical fullscale reactor is also included.This paper focuses on feedback control with on-off control and PI control (proportional + integral), and feedforward control.An advanced alternative to these traditional control methods is MPC (model-based predictive control) which is possible to use for the present reactor control since a fairly accurate dynamic reactor temperature model exists and the crucial variables are measured.However, MPC is not included in this paper since the temperature setpoint is constant and disturbances as ambient temperature changes can be compensated for effectively using the aforementioned control methods.However, MPC may be the preferred control method in applications where the setpoint (profile) is changing, and in applications where the trade-off between small control error and smooth control actions should be directly adjustable by the user.
The outline of this paper is as follows.A process description is given in Section 2. In Section 3, rationales of bioreactor reactor temperature control are given.On-off temperature control is described in Section 4. Smooth control with PI control is covered in Section 5. Model-based and model-free feedforward control of a simulated reactor are described in Section 6. Results for temperature control of a full-scale reactor are provided in Section 7. A discussion is given in Section 8, and conclusions are given in Section 9. A dynamical mathematical model used for analysis and simulation is presented in Appendix A. Abbreviations and nomenclature are given in Appendix B.
Energy recovery design, e.g. using relatively warm reactor liquid effluent to preheat cold influent, and using biogas combustion to heat the reactor, is not covered in this paper.Energy recovery and other issues in optimal reactor design and operation will be addressed in future publications.
MATLAB and SIMULINK (MathWorks, Inc.) are used for numerical computations and simulations based on the models described in Section A. The software of the real temperature control system is implemented in LabVIEW (National Instruments, Inc.).The feed temperature, T feed , is assumed to be the same as the ambient -here: room -temperature, T amb , since the feed resides for several days in an intermediate storage in the room where the reactor is placed.

Process description
Detailed information about the system components are given in the following.
Reactor is cylindrical with 250 L effective liquid volume.Height (0.40 m) is 5 times reactor diameter (2.00 m).Gas volume is assumed negligible compared to liquid volume.
TT-1 is a Pt100 reactor temperature sensor.It has an accuracy of approximately ± 0.3 • C. The repeatability is not known.
The temperature measurement signal is noisy.It is assumed that the noise, n, is a random stochastic signal with zero mean value.From a representative unfiltered measurement time-series, TC-1 is the temperature controller implemented in LabVIEW (National Instruments) running on a PC.The available controller functions are manual control, PID control, and on-off control.The timestep of the control loop is 2 s.
H-1 is an electrical heater for the AD reactor which is controlled using pulse-width modulation (PWM) option in TC1.The heater comprises an electrical resistor wound around the reactor inside the thermal insulation jacket.The maximum power delivered by the heater is 200 W.
PWM is a pulse-width modulation element implemented with the Square Wave Point-by-Point function in LabVIEW.The PWM element operates with a fixed cycle time of 30 sec which is negligible compared to the dynamics of the temperature control loop.The control signal (in percent) calculated by the temperature controller is the duty cycle of the PWM element.PWM control emulates analog control.
SSR is a solid state relay (semiconductor) which is turned on-off with a voltage (5/0 V) corresponding to the state of the PWM element (on-off) which controls the SSR.The SSR switches the 220 VAC mains voltage onto/off the heater.
A secure communication between the PC and the Internet with the LogMeIn software facilitates remote access to the computer screen and to the file system on the lab PC.

Rationales of bioreactor temperature control
In the following, rationales of bioreactor temperature control are given in terms of temperature dependency of methane production and temperature disturbance compensation.

Temperature dependency of methane production
For a bioreactor, the produced methane flow depends on the reactor temperature, T reac .The temperature dependency is expressed in e.g.Hill's model of AD of animal wastes, Hill (1983).In Hill's model the maximum reaction rate, µ max , of the acidgenerating microorganisms -acidogens -and the methane-generating microorganisms -methanogensis temperature-dependent.This dependency is represented by the following linear function, Hashimoto et al. (1981): In the ADM1 model (Anaerobic Digestion Model No. 1), Batstone et al. (2002), the temperature dependency of a number of model parameters is expressed in terms of Arrhenius-like functions.As an illustration of the dependence of F meth on T reac , Figure 2 shows measured (and simulated) timeseries in F meth due to changes in T reac and in F feed for the reactor ADR1 which was in use at Foss Biolab from August 17, 2011until April 19, 2012.The rest of the present paper focuses on reactor ADR2 which has been in use from April 19, 2012.However, it is assumed that the temperature dependency holds equally well for ADR2 as for ADR1 since the physical appearances of the two reactors are similar and the operation and feedstock (waste from the same dairy livestock) are similar.
During the time interval in Figure 2 both T reac and the feed flow F feed were changed, but only the variations caused by the temperature change is of interest here.T reac is increased twice (implemented as step-wise changes of the temperature control setpoint):  The simulations of F meth shown in Figure 2 are based on the modified Hill's model adapted to ADR1.Adaptation of modified Hill's model to ADR1 is not published.(Model adaptation to ADR2 is presented in Haugen et al. (2013a).) The temperature effect on gas production is a result of gas solubility changes, reduced microbial growth rates, and stress caused by the temperature transition, Tchobanoglous et al. (2003).
Tchobanoglous et al. state that methanogens are sensitive to temperature changes, and that these microbes should not be excited to temperature variations larger than ± 0.5 • C. Consequently, a temperature control system should be designed to be able to keep the temperature offset from the specified temperature setpoint less than ± 0.25 • C.
Tchobanoglous et al. also point out that most AD processes are designed for operation at mesophile conditions, i.e. at temperatures in the range 30-38 o C. While it is important to determine the optimal reactor temperature, this is not address in the present paper, but instead in a forthcoming paper based on theoretical optimization methods applied to mathematical models of the AD reactor.
The most important disturbances acting on the reactor temperature are • the ambient temperature, T amb , • the temperature of the feed flow, T feed (assumed to be the same as T amb for the practical reactor), • the feed flow, F feed .

Disturbance compensation
A well operating temperature control system will compensate for changes in these disturbances automatically.To demonstrate the importance of temperature control, Figure 3 shows the responses in T reac , u, and T amb with (automatic) control with and without control for reactor ADR2.The setpoint, T sp , is 30

On-off control
The on-off controller can be regarded as the simplest feedback controller available.The controller function is where e is control error: and d e an adjustable dead-band to avoid switching of u due to (measurement) noise in e. Noise-triggered switching is also counteracted using a measurement lowpass filter with a properly adjusted time-constant, τ f , to attenuate the noise properly.In the present application, τ f = 10 min, found by trial-and-error, and d e = 0.

Simulations
Figure 4 shows simulated responses of T reac and u with T sp = 30 o C. T amb is set to 17 o C which is representative for the room temperature in the real experiment reported below.u on = 100%.u off = 0%.The time interval of the plot is 0.32 d which is the same as for the real time-series plotted in Figure 4.
• Amplitude of oscillation of T reac is 0.04 o C.

Practical results
Figure 5 shows experimental T reac , T sp , u, and T amb .
• Amplitude of oscillation of T reac is 0.05 o C.

Comments and conclusions
• The simulated responses are in good accordance with the real responses, which indicates that the dynamic model used for simulation is quite accurate.
• For the real responses: |e| max ≈ 0.07 o C is acceptable.Also, µ e = 0.019 o C is acceptable.The oscillatory behaviour of T reac is acceptable since the variation is within ±0.25 o C.
• The on-off behaviour of the control signal, u, is also acceptable in our application since the actuator is an electrical heater with no moving parts.However, in applications with a mechanical actuator like a pump or valve used to manipulate the flow of e.g.hot water or steam, smooth or continuous control with PI(D) control may be preferred, cf.Section 5.

Controller function
PID control is prevalent in industrial applications, Seborg et al. (2004).The PID controller provides smooth control as opposite to on-off control.The PID controller used in this paper is based on discretization of the following continuous-time PID controller: The discretization is based on the implicit Euler method with time-step τ s = 2 s.Typically, the derivative term provides control stability and agility, but it also propagates measurement noise which may cause too noisy control signal.In the present application, representative time-series of the raw (unfiltered) temperature measurement show a control signal standard deviation of 1.3 K with PI control and 4.5 K with PID control.Although the actuator in the present application is an electrical heater with no moving parts so the control noise can not make any mechanical problems, it was decided to not use the derivative term.This decision is made to increase the relevancy of controller tuning results in the present paper to systems with mechanical actuators.
A great number of controller methods for tuning controller parameters exist, Seborg et al. (2004), O'Dwyer (2003).Selected open loop controller tuning methods applied to the reactor are presented in Section 5.3, while applications of selected closed loop methods are presented in Section 5.4.Summaries of tuning results are given in Section 5.2.

Summary of results with various tuning methods
Sections 5.2.1 and 5.2.2 below summarize the results of controller tuning for the simulated temperature control system and for the real system, respectively.Tuning details are in Sections 5.3 and 5.4.

Simulated temperature control system
Table 1 summarizes the results for the simulated temperature control system based on the model presented in Appendix A.1.The table shows controller settings, GM (gain margin), phase margin (PM), and the closedloop response-time τ r [d] which is calculated as the inverse of the bandwidth which is here defined as the amplitude crossover frequency, ω c [rad/d]: τ r indicates the speed of the response of the control system due to a setpoint step change.τ r is approximately the time-constant of the control system.The above frequency response characteristics are based on the transfer functions model in Appendix A.2. Seborg et al. (2004) recommend the following ranges for the stability margins, where the lower limits can be  2004), the Hägglund-Åstrøm Robust Tuning method, Hägglund and Åstrøm (2002), and the SIMC method (Simple IMC) by Skogestad (2004), here denoted the Skogestad method.The Ziegler and Nichols open loop method has no adjustable settings, and typically give very fast control but with relatively small stability margins.The Hägglund-Åstrøm method has no adjustable settings.The IMC, Lambda, and the Skogestad method each has one tuning parameter which determines the closed loop time constant, and typically for these methods the setpoint step responses are without oscillations indicating relatively large stability margins.
It is convenient to use a tuning method.This leaves out the Ziegler and Nichols open loop method and the Hägglund-Åstrøm method.Among the remaining candidates, the Skogestad method is selected as we are not aware of important benefits with the other methods over the Skogestad method.It is evaluated favourably in Haugen (2010) comparing with a number of open loop methods, and closed loop methods.

The Skogestad method (SIMC method)
Skogestad ( 2004) has developed PID controller tuning formulas for a number of processes given by their transfer functions.As shown in Appendix A.2, the reactor can be represented by a transfer function comprising a dominant time-constant term representing the energy balance of the reactor liquid with some additional lag.In controller tuning, it is safe regarding control system stability to assume that this lag is a time-delay of the same amount as the lag, Skogestad and Postlethwaite (2007).Thus, it is safe to use the Skogestad PI tuning formulas for the following "time-constant with time-delay" model of the reactor: For this type of process, Skogestad designates a PI controller.As pointed out in Haugen and Lie (2013), cf. also DiRuscio (2010), the Skogestad PI tuning formulas for eq.( 9) become identical with the tuning formulas for the following "integrator with time-delay" process which approximates eq. ( 9) in the transient phase: where The Skogestad PI settings for the process model eq.
(10) become Comments: In eqs.( 13) and ( 15), Skogestad's rule-ofthumb τ c = τ delay are used.The factor 2 in eq.14 and 4 in eq.15 are due to the modification of the τ i setting introduced in Haugen and Lie (2013) to give faster disturbance compensation, while retaining acceptable stability margins.In Skogestad's original settings, the factors are 4 and 8, respectively.The process parameters K ip and τ delay can be found experimentally, or from a model, as explained in the following.
Estimating τ delay from an experimental response Figure 6 shows the response in the T reac due to a step change in the control signal, u, from 62% to 82%, hence a step amplitude of ∆u = 20%.The measurement filter which is normally in use, is bypassed in this experiment to obtain parameter values that are independent of the filter dynamics.Figure 6 also shows (in red colour) the ramp-like response adapted to the real T reac .The slope of this ramp-like response determines parameter K ip , as explained below.
From the response shown in Figure 6 which is without measurement filtering, a lag of approximately 0.01 d can be observed.Under normal operation of the reactor a measurement filter of time-constant 15 min = 0.0069 d is used.Thus, the total lag is approximately 0.01 + 0.0069 = 0.0169.Furthermore, the sampling time of the data shown in Figure 6 is 15 min = 0.0104 d which adds uncertainty to the estimation of the aforementioned lag of 0.01 d.Consequently, a total lag, τ lag , is estimated visually as The estimate eq. ( 16) is in good accordance with the lag estimated with a nonlinear least square method as 0.023 d in Haugen et al. (2013a).There, the lag was estimated at a lower feed rate, namely 45 L/d, while the feed rate in the present study is 65 L/d.As argued in the beginning of the present section, Estimating K ip from an experimental response K pi can be found as the normalized initial slope of the step response in the reactor temperature: where S is slope and ∆u is amplitude of step change of u.A step response test can be accomplished during a few hours, while it may take several days to obtain K pi from eq. ( 11) if K and τ reac are estimated from a step response since the (theoretical) τ reac for the reactor is typically several days (in the operating point defined in Table 6 it is 2.24 d).
For the present reactor, K pi is found as follows.Figure 6 shows (in red colour) the ramp-like step response adapted to the real T reac over the time-interval t 0 = 340.94-341.20 d using the following assumed model for this ramp: where t [d] is time.The coefficients a and b are estimated with the least squares method.However, only a is of interest here.It is estimated as Now, K ip can be calculated from eq. ( 11): Calculating K ip from the reactor model K ip can be calculated from eq. ( 18) where K and τ reac can be calculated from the transfer function derived from the energy balance of the reactor, cf.Appendix A.2. From eq. ( 11), using eqs.( 53) and ( 54), which agrees very well with the experimental value in eq. ( 23).

Simulations
The PI settings are calculated with K pi given by eq. ( 24) and τ delay given by eq. ( 17).The PI settings are shown in Table 1. Figure 7 shows simulations of the control system.The responses indicate acceptable stability.The stability margins shown in Table 1 have acceptable values, though GM is large.

Practical results
PI settings are calculated with eqs.( 13)-( 15) with K pi given by eq. ( 23) and τ delay given by eq. ( 17).The resulting settings are shown in Table 2. Figure 8 shows responses on the real reactor with these settings.
With the above Skogestad PI settings the standard deviation of T reac (10 min time-constant filter) over 20 days is 0.015 K.The mean of T reac is very close to its setpoint.The variations of T reac are within approximately ±0.05 o C. 1

Comments and conclusions
With the Skogestad tuning method: • The tuned control loop shows good stability.
• The tuning experiment does not involve any trialand-error, i.e. iterations are not needed, which is beneficial from a practical point of view.
1 Due to external demands for the operation of the reactor,  Closed loop tuning methods are applied with the controller in place (in the loop).The following closed loop methods are considered: • The well-known Ziegler-Nichols (ZN) closed loop method, Ziegler and Nichols (1942), with Åström-Hägglund's relay-method, Åstrøm and Hägglund (1995) to find the ultimate gain and period.(Section 5.4.2.) • The Relaxed Ziegler-Nichols (R-ZN) closed loop PI tuning method, proposed by Haugen and Lie (2013).This method is based on the same experiments as in the ZN (closed loop) method, but relaxes the PI settings to obtain a smoother control signal and to improve the stability compared with the original ZN method.The method is based on a combination of the Skogestad method and the ZN closed loop method.(Section 5.4.3.) • The Tyreus-Luyben method, Tyreus and Luyben (1992), which is, probably, the best known method to modify the ZN closed loop PI settings to obtain more relaxed control.However, it is shown by Haugen and Lie (2013) that the R-ZN method is beneficial compared with the Tyreus-Luyben  method.These benefits are confirmed in simulations of the reactor (detailed results are not shown here).
• The Good Gain method, Haugen (2012), which has similarities with the ZN method.Sustained oscillation in the tuning phase is avoided, and in addition the final stability of the control system is typically improved comparing with the ZN method.However, the method can be used reliably only if the noise and disturbances affecting the process measurement is small to make it possible to read off the tuning parameter T ou (time from overshoot to undershoot after a setpoint step with a P controller).On the real reactor the noise and disturbances are so prevalent that the Good Gain method is not applicable.This problem is confirmed in simulations containing realistic noise (responses are not shown).

The ZN closed loop method based on relay tuning
Åstrøm and Hägglund (1995) suggest a relay or onoff controller to replace the P controller in the tuning phase of the ZN closed loop (or Ultimate Gain) method, Ziegler and Nichols (1942), thereby avoiding the trial-and-error procedure since the oscillations come automatically.The ultimate controller gain is calculated as where A is the amplitude of the on-off control signal.
If u on = 100% and u of f = 0%, as in our application, A is 50%.E is the amplitude of the oscillations in the process measurement.
The PI controller settings are where P u is the period of the oscillation.
Due to external demands it was necessary to operate the reactor at approximately 30 o C in the experiments with the ZN method, while 35 o C was used in experiments with the Skogestad method.

Simulations
The simulations with on-off controller shown in Figure 4 are the basis for relay tuning.From the simulations, E = 0.04 K and P u = 0.055 d.Furthermore, A = 50%.This gives PI settings as shown in Table 1, where also stability margins, and response-time are shown.The resulting stability margins, cf.Table 1, are very small.Although not shown here, simulations show oscillatory responses, with little damping.

Practical results
The real responses with on-off controller shown in Figure 5 are used for relay tuning.From the responses, E = 0.05 K and P u = 0.045 d.Furthermore, A = 50 %.This gives The resulting PI settings are calculated with eqs.( 27)-( 28) to give PI as shown Table 2. Figure 9 shows responses on the real reactor.

Comments and conclusions
• Both the model and the real system shows poor stability.This poor stability is actually typical when the ZN tuning is applied to a process where there is a small or no pure time-delay, as is the case here.
• It is concluded that the ZN closed loop method is inappropriate for tuning the temperature controller.

Relaxed ZN PI tuning
Relaxed ZN PI settings, as proposed by Haugen et al. (2012b), are calculated from the ultimate gain, K cu , and the ultimate period, P u , found from e.g.relay oscillations: and where k r is a parameter determined by the user to obtain a proper closed loop system time-constant, where τ delay is the process time-delay.k r = 1 kan be regarded as the default value.With k r = 1 eq.( 32) is the same as Skogestad's rule-of-thumb: τ c = τ .Enhanced relaxed control can be obtained with k r > 1. Haugen et al. (2013b) recommend k r = 1 in eqs.( 30)-(31) if the process has a dominating lag or integrator, due to energy or material balance, plus a noteable time-delay, and k r = 4 if the process has zero or neglible time-delay, but some lag, in addition to the dominating lag or integrator.The bioreactor has a dominating lag -approximately an integrator -due to the energy balance of the liquid, and an additional relatively small lag due to dynamics in the heater and the reactor wall.There is also a relatively small lag due to the measurement filter.A physical reason for a clear pure time-delay is not obvious.Thus, k r is set to 4 in the PI settings given by eqs.( 30

Practical results
k r = 4 is used in eqs.( 33)-( 34) with K cu = 1273 %/K and P u = 0.045 d from Section 5.4.2.The resulting PI settings are shown in Table 2. Figure 10 shows responses on the real reactor with these PI settings.
(The PI settings with k r = 4 were applied just before the setpoint step.)The responses indicate acceptable stability.Also, k r = 1 is applied on the real reactor, but responses (not shown here) indicate poor stability.

Comments and conclusions
• Both theoretical analysis, i.e. simulations (though not shown here) and stability margins, and practical results indicate successful controller tuning using enhanced R-ZN settings with k r = 4.The stability margins are large, cf.Table 1, and the simulated and real responses are smooth.
• R-ZN settings with k r = 1 are not recommended here.
5.5.Control system robustness against process parameter changes

Introduction
The transfer function model of the temperature control system presented in Appendix A.2 forms a good basis for a stability robustness analysis of the control system.It is assumed that the controller is a PI controller tuned with the Skogestad method at one specific operating point.The Skogestad model-based PI settings formulas The PI settings are given by eqs.( 13) and ( 15).Assuming K ip is given by eq. ( 24), the PI settings become In the following subsections the control system robustness against changes in parameters assumed most apt to changes, is discussed.The changes are: • Changes in the feed flow, F feed .
• Changes in the reactor lag, τ lag .
The impact that changes in these two parameters have on the dynamic properties of the control system is analysed.To this end, it is assumed that the controller, which is assumed a PI controller, is tuned with the Skogestad method as in Section 5.3.2.

Changes in feed flow
The tuning is based on the process having "integrator with time-delay" dynamics, cf.Section 5.3.2.These settings are valid as long as the time-constant is larger than four times the time-delay, and this assumption is always valid for a practical reactor -even with a varying F feed .As an example, assume the relatively high value F feed = 87 L/d which is the feed flow which gives the maximum methane gas flow in steady state as calculated from Hill's model adapted to the present bioreactor by Haugen et al. (2013a).The reactor timeconstant is then which is much more than four times the effective lag of 0.02 d used as a time-delay in the Skogestad tuning method.So, the above PI settings, which are independent of F feed , apply even if F feed has its largest value.Obviously, they also apply for the smallest resonable value of F feed since a small value makes a relatively large value of τ reac .In other words, the stability of the control loop is essentially independent of F feed since the assumptions of Skogestad's PI tuning rules remain valid.
Assuming that variations of F feed are the only parameter variations which affect the reactor dynamics, it is concluded that there is no need for adjusting the PI settings as functions of the varying F feed .This is also confirmed in simulations.

Changes in lag or time-delay
Although the transfer function model presented in Appendix A.2 does not contain any pure time-delay transfer function, it is useful to assume that such a timedelay is present since the time-delay margin is a safe (conservative) estimate of the lag margin.It can be shown, see e.g.Haugen and Lie (2013), that the timedelay margin (increase), ∆τ , can be calculated from the phase margin, PM, with eq. ( 38) below.Inserting numbers related to PI controller tuning with the Skogestad method given in Table 1, yields the results given in eq. ( 40) below.
One implication of this value of ∆τ is that stability problems may occur if the measurement filter timeconstant, τ f , is increased by an amount approximately 19.5 min.If it is necessary to increase τ f , it should be accompanied by an equal increase in τ delay used in the Skogestad tuning formulas, eqs.( 13) and (15).

Feedforward control 6.1. Introduction
Feedforward control can compensate very effectively for variations in process disturbances, Seborg et al. (2004).The following variables are regarded as disturbances acting on the bioreactor here regarded as a thermal system: T amb , T feed , and F feed .
Figure 11 shows the structure of a temperature control system for the bioreactor with both feedforward and feedback control.The total control signal is calculated as the sum of the feedback and feedback control terms: Results presented in previous sections indicate that for the reactor studied, feedback control is sufficient to keep the reactor temperature close to the setpoint.Hence, there is hardly any need for feedforward control.However, in other cases with severe, varying disturbances due to e.g.large ambient temperature variations, feedforward control may give a substantial improvement of the control.It will be shown how to design feedforward control for the present bioreactor, and the results should be transferable to other reactors or similar thermal systems.
In the following respective sections, two alternative feedforward controllers are developed: • Model-based feedforward controller using a phenomenological model, i.e. an energy balance of the reactor.
• Model-free feedforward controller using steadystate operational data only.
Simulation results are shown in the following.However, no practical results are shown since feedforward control is not implemented on the real system.

Model-based feedforward control
The feedforward controller can be designed from the process model, eq. ( 46), as follows: First, the reactor temperature T reac is substituted by its setpoint T sp .Then the resulting model is solved for the control variable u, now denoted u ff , to get the feedforward controller: which can be implemented assuming all parameters and variables on the right-side of eq. ( 42) are known apriori or from measurements, which is a realizable assumption here.

Simulations
The feedforward controller, eq. ( 42), is applied to a simulated reactor having model parameter values as shown in Table 6.The setpoint is constant.T amb = T feed (as in all subsequent simulations) is varied as a sinusoid of amplitude 10 o C with a period of 1 d, which assumed a representative variation if the reactor is outdoors, and with mean value 15 o C. In the simulation, T amb pass through a lag of time-constant 0.01 d before it enters the contents of the reactor.This lag is meant to represent additional thermal dynamics of a real reactor.However, this lag is not included in the feedforward controller.Hence, a (relatively small) model error is included.
Figure 12 shows responses with and without feedforward control.Both cases include feedback PI control with the controller tuned with the Skogestad method with K c = 152 %/ o C and T i = 6912 s.
Table 3 shows analysis results.|e| max is the maximum control error.The IAE index is calculated from t = 0.5 to 5 d.Note: When the feedforward controller is active, it is necessary to set the output range of the PI controller to cover both positive and negative values to make the PI controller be able to compensate for the imperfect feedforward control signal, which is due to model error, with both positive and negative values.If the output range is only positive, the compensation may be insufficient, and the result may be a nonzero steady state control.If the output range of the PI controller can cover positive values only, an alternative solution is to subtract a proper negative constant, e.g.20%, from the PI control signal, thereby forcing the PI control signal to become positive.

Comments and conclusions
• Feedforward control improves the control performance considerably.
• The simulations show that the control signal timeseries appears very similar in feedback (only) control and in feedforward control, indicating that the "timing" of the control action is crucial for good control performance, and good timing is provided by feedforward control.
• The reason why there is a nonzero control error with feedforward control is the inclusion of the assumed realistic thermal dynamics in terms of a lag of 0.01 d.Without this model error, the control error would have been zero with feedforward control.

Model-free feedforward control
A feedforward controller may be designed from steadystate operational data.It is assumed here that T amb = T feed is the most important varying disturbance for our reactor.This is an approximate design method since it is based on only steady-state data, but it can improve the disturbance compensation substantially.If the disturbance is a so-called input disturbance, i.e. the disturbance enters the process dynamically at the same "position" as the control variable does, the model-free feedforward may perform as well as the model-based feedforward.This is the case in the simulation example described below.

Simulations
Table 4 shows the corresponding steady-state values of T amb = T feed and u found under steady-state conditions with PI control (subindex "s" means steady-state).The conditions are the same as in the simulation in Section 6.2.The simulated responses with model-free feedforward is virtually indistinguishable from the results with model-based feedforward shown in Figure 12 and Table 3.In this section results for the pilot reactor presented in previous sections are applied to a simulated full-scale reactor having the same form of mathematical model as for the pilot reactor, cf.Appendix A, with the following parameter and operational values: • The reactor volume is V = 10 m 3 .This size is assumed representative for reactors at farms using animal waste as feed.
• The reactor is assumed rectangular with height, H, and depth, D, being equal, and width, W , being twice the depth, as is approximately the case for anaerobic baffle reactors (ABR). 2 From the known volume, H = D = W/2 = 1.71 m.
• The area-specific heat conductivity is same as for pilot reactor.Hence, the heat conductivity of fullscale reactor is G fs = GA fs /A pilot = 2.08 • 10 6 (J/d)K.G is conductivity of pilot reactor.A fs and A pilot are conductive areas of the respective reactors.A pilot is calculated from the given volume and design of the pilot reactor, cf.Section 2 (detailes not given here).Assuming for simplicity that all areas are conductive, A fs is calculated as 10H 2 .
• The reactor lag is guessed as τ lag = 0.05 d, while it is 0.01 d for the pilot.
• The temperature measurement filter timeconstant is as for pilot reactor, τ f = 10 min = 0.0069 d.
• • Using the extreme operating point (above) in the static version of the dynamic energy balance eq.( 46), and allowing for 50% design margin, the maximum power to be delivered by the electrical heater is 29.8 kW.
• The controller output, u, is in unit of kW, not percent as for the pilot.
2 An ABR reactor of this size is being constructed at Skoglund farm, Porsgrunn, Norway.
The main specification of the temperature control system is: • The control error is limited, cf.Section 3: The following conditions are assumed for simulation and analysis: T sp = 38 o C, T feed = 0 o C, F feed = 10000 L/d.T amb is assumed sinusoidal with period 1 d, assumed to represent a relatively large, still realistic, outdoor temperature variation: where t is time [d].Table 5 summarizes the results of simulations and analysis for different controllers and controller settings.The results are commented in the following.

Simulations
Simulations have been run, but plots are not shown here, with the following three controllers: • PI controller tuned with the R-ZN method with tuning parameter k r = 4, cf.Section 5.4.3.
• PI controller tuned with the Skogestad method, cf.Section 5.3.2.

Frequency response analysis
Figure 13 shows for each of the two PI controller settings, a Bode plot of the magnitude of the sensitivity function of the control loop, S(jω) = S(j2πf ) = 1/ [1 + L(j2πf )] where L is the loop transfer function defined in Appendix A.2.The frequency f 1 = 1 d −1 is the frequency of the sinusoidal T amb in the simulations.The |S|-plot shows how much the response in the process output variable due to a sinusoidal process disturbance is reduced by using feedback control, compared with no feedback control, Seborg et al. (2004).The smaller the value of |S(j2πf j )|, the more effective the feedback disturbance compensation for a sinusoidal disturbance of frequency  43) if no other disturbances exist.However, in practice, feedback control is needed to compensate for static or low-varying disturbances.

Comments and conclusions
• The values of GM and PM indicate that the stability is acceptable with both PI settings.
• The Bode plot in Figure 13, and Table 5, shows that with Skogestad PI settings, |S(j2πf 1 )| = 0.88.This indicates that the feedback control loop reduces the impact of the assumed sinusoidal T amb on T reac by only 12%.With the R-ZN settings the reduction is far better, namely 65%.
• Inequality ( 43) is not satisfied with the on-off controller due to a permament mean offset from setpoint.However, ineq.( 43) may be satisfied if the value of u on is reduced.
• Assume that the ineq.( 43) is not satisfied with PI control.An attempt to optimize the controller tuning can be made using loop-shaping, Skogestad and Postlethwaite (2007), or optimization methods, Edgar et al. (2001).The controller settings must satisfy the following requirement: where A T amb is the maximum amplitude of T amb , e.g. 10 o C, E is given by ineq.( 43), and f 1 is the frequency [d −1 ] of the sinusoidal T amb , here assumed 1 d −1 .
Also, activating the derivative term should be considered to reduced the control error.
• The control error may be reduced considerably with feedforward control, cf.Section 6. Modelbased feedforward is implemented on the simulated full-scale reactor with the aformentioned result, but responses are not shown here.However, the results above indicate that feedback PI control is sufficient.

Model accuracy
For the pilot reactor, the practical performance of the reactor temperature control systems with on-off control and with PI control is in good accordance with the theoretical performances as seen in simulations.This indicates that the mathematical models used are sufficiently accurate to be used for analysis, design, and simulations.The accuracy of the models in the present study motivates for use of models for design of planned reactors having different physical dimensions.The present reactor is heated by an electrical resistor wound around the reactor inside the thermal insulation jacket.If the reactor is heated differently, e.g. by heating the influent, we think that just simple modifications of the model are necessary.

Sensor accuracy
According to technical specifications the reactor temperature sensor, a Pt100 sensor, has an accuracy of approximately ± 0.3 o C in the pertinent temperature range.In various experiments that are conducted in this study, the observed temperature responses vary less than this accuracy.Although we have no data for the repeatability of the sensor, the good accordance between measured and simulated responses indicate that the repeatability is sufficient for us to rely on the temperature measurements.

Conclusions
It is demonstrated that the produced methane gas flow depends clearly on the bioreactor temperature.Moreover, according to literature references, methanegenerating microbes should not be exposed to temperature changes larger than 0.5 o C in amplitude.Thus, the temperature should be controlled to a setpoint with a maximum control error of ± 0.25 o C.
On-off feedback control may be used for temperature control, given that the on-off operation is acceptable, which may not be the case with a mechanical actuator.The mean control error (offset) from the setpoint is typically non-zero.The maximum control error may be unacceptable.
A PI controller can be tuned successfully with the Skogestad method and with the R-ZN tuning method.The original ZN closed loop method is not appropriate because of poor resulting stability.
The robustness of the PI control system is investigated assuming model-based Skogestad PI settings.The PI settings are independent of the feed flow, so the tuning is robust against feed flow variations.Frequency response analysis shows a time-delay margin of approximately 20 minutes which is here assumed a safe value.
Both model-based feedforward controller designed from the energy balance of the reactor, and a modelfree controller using table lookup on operational data, are applied to a simulated reactor, with almost identical performances.Comparing with only feedback control, feedforward control improves the temperature control considerably.
A temperature control system for a simulated fullscale reactor is simulated.The self-regulation of the reactor is sufficient to limit the impact on the reactor temperature by an assumed large sinusoidal daily variation of the ambient temperature.In practice, feedback control is needed to compensate for static or lowvarying disturbances.
Figure 6 shows the response in T reac lag due to a step in the control signal, u.From this response the lag is estimated visually with the value given in Table 6.

A.2. Transfer functions model of the temperature control system
For simplicity, the same symbol is here used for a variable in both the Laplace domain and in the timedomain.
Figure 14 shows a block diagram of a transfer functions model of the control system.This model is used as the basis for the frequency response analysis in Section 5.
The individual transfer functions of eq. ( 50) are defined in the following.

Controller transfer function
Assuming

Filter transfer function
The raw temperature measurement, T mr , which here is the same as T reac lag , is quite noisy and is therefore filtered with a time-constant filter having the following transfer function model: where τ f is the filter time-constant.For the present reactor, τ f is set to 10 min.Note: Above, T mf represents the filtered measured reactor temperature.However, in most sections in this paper, it is practical to use symbol T reac for filtered measured reactor temperature.

Figure 2 :
Figure 2: Reactor ADR1: Responses in F meth (middle) due to changes in T reac (lower) and F feed (upper).For F meth : Measured is blue.Simulation is red.[This plot also appears in Haugen et al. (2013a)].

•
At time t = 60.5 d (days): From 24 o C to 30 o C. • At t = 67.5 d: From 30 o C to 35 o C.

Figure 2
Figure 2 illustrates clearly the dependence of F meth on T reac .The simulations of F meth shown in Figure2are based on the modified Hill's model adapted to ADR1.Adaptation of modified Hill's model to ADR1 is not published.(Model adaptation to ADR2 is presented inHaugen et al. (2013a).)Thetemperature effect on gas production is a result of gas solubility changes, reduced microbial growth rates, and stress caused by the temperature transition,Tchobanoglous et al. (2003).Tchobanoglous et al. state that methanogens are sensitive to temperature changes, and that these microbes should not be excited to temperature variations larger than ± 0.5 • C. Consequently, a temperature control system should be designed to be able to keep the temperature offset from the specified temperature setpoint less than ± 0.25 • C.Tchobanoglous et al. also point out that most AD processes are designed for operation at mesophile conditions, i.e. at temperatures in the range 30-38 o C. While it is important to determine the optimal reactor temperature, this is not address in the present paper, but instead in a forthcoming paper based on theoretical optimization methods applied to mathematical models of the AD reactor.The most important disturbances acting on the reactor temperature are

Figure 4 :
Figure 4: Simulated time-series for the reactor with onoff temperature control.Random measurement noise is included.

Figure 5 :
Figure 5: Real time-series for the reactor with on-off temperature control.T reac is filtered with time-constant 10 min.The sampling timestep of time-series is 15 min.

Figure 6 :
Figure 6: The response in T reac (upper plot, blue) due to a step change in the control signal, u, (lower diagram) from 62% to 82% (∆u = 20%), and the ramp-wise response (upper plot, red) adapted to T reac .(The sampling time of the data is 15 min = 0.0104 d.)

Figure 7 :
Figure 7: Simulations of temperature control system with the Skogestad PI settings.

Figure 8 :
Figure 8: Real time-series for reactor with PI temperature controller tuned with the Skogestad method.

Figure 9 :
Figure 9: Real time-series for reactor with PI temperature controller tuned with relay-based ZN method.

Figure 10 :
Figure 10: Real responses in the temperature control system with enhanced R-ZN PI settings with k r = 4.

Figure 11 :
Figure 11: Temperature control system for the bioreactor with both feedforward and feedback control.

Figure 12 :
Figure 12: Simulation of temperature control with and without feedforward control, with feedback PI control in both cases.
An extreme operating point, with maximum power demand, is assumed: Temperature setpoint is T sp = 38 o C. Ambient temperature is T amb = −20 o C. Temperature of liquid feed is T feed = 0 o C. Feed flow F feed = 1000 L/d giving hydraulic retention time HRT = 10000 L/10000 L/d = 1 d.

Figure 13 :
Figure 13: Bode plot of the magnitude of the sensitivity function of the control loop for two different PI controller settings.
f j [d −1 ].If |S(j2πf j )| ≈ 1 = 0 dB, the feedback makes no difference compared with open loop control (constant control signal).Let H d (s) be the open loop transfer function from T amb to T reac .H d (s) can be calculated from eq. (46).|H d (j2πf j )| expresses the self-regulation of the reactor for sinusoidal T amb of frequency f j [d −1 ].Assume f j = 1 d −1 = f 1 , corresponding to eq. (44).Calculations show that |H d (j2πf 1 )| ≈ 0.0079 for F feed between the assumed large value of 10 m 3 /d and small value of 1 m 3 /d.Hence, without feedback control, the amplitude of T reac is 0.0079 • 10 o C = 0.079 o C which satisfies ineq.(

Figure 14 :
Figure 14: Block diagram of a transfer functions model of the control system.
Seborg et al. (2004).2.Real temperature control systemTable2summarizes the controller settings for the real temperature control system.Table2: Results with various PI tuning methods for the real temperature control system: S = Skogestad.ZN = Ziegler-Nichols.R-ZN = Relaxed Ziegler-Nichols.Seborg et al. (2004), Internal Model Control (IMC) methods,Seborg et al. ( in the pertinent time-series, while Treac sp is varied around 35 o C in Figure8.Filtered meas T reac = blue.Setpoint T reac,sp = red.

T
• For each of N distinct values of the disturbance, T amb , observe the value of the control signal u which gives approximately zero steady state control error.This may be done during PI(D) feedback control.Typically, feedback control is used together with feedforward control, so no extra effort is needed to run the feedback control here.The result is a table of corresponding values of T amb and u s (steady-state value).• Use table lookup, i.e. some interpolation method, to calculate the instantaneous feedforward control signal u ff from the instantaneous measured T amb .

Table 4 :
N = 5 corresponding values of T amb = T feed

Table 5 :
Results with three different controllers.R-ZN PI = PI controller with R-ZN settings.S PI = PI tuned with Skogestad settings.On-off = On-off controller.