Bioreactor profile control by a nonlinear auto regressive moving average neuro and two degree of freedom PID controllers

Journal of Process Control 24 (2014) 1761–1777

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Bioreactor profile control by a nonlinear auto regressive moving average neuro and two degree of freedom PID controllers Ubaid Imtiaz a , Sudhanshu S. Jamuar b,∗ , J.N. Sahu c , P.B. Ganesan a a b c

Department of Electrical Engineering, University of Malaya, Faculty of Engineering, 50603 Kuala Lumpur, Malaysia School of Microelectric Engineering, University Malaysia Perlis, Malaysia Department of Petroleum and Chemical Engineering, Faculty of Engineering, Institut Teknologi Brunei, Brunei Darussalam

a r t i c l e

i n f o

Article history: Received 24 March 2014 Received in revised form 3 August 2014 Accepted 23 September 2014 Available online 16 October 2014 Keywords: Bioreactor profile Inverse neural network Process control NARMA neuro controller

a b s t r a c t This paper presents the use of nonlinear auto regressive moving average (NARMA) neuro controller for temperature control and two degree of freedom PID (2DOF-PID) for pH and dissolved oxygen (DO) of a biochemical reactor in comparison with the industry standard anti-windup PID (AWU-PID) controllers. The process model of yeast fermentation described in terms of temperature, pH and dissolved oxygen has been used in this study. Nonlinear auto regressive moving average (NARMA) neuro controller used for temperature control has been trained by Levenberg–Marquardt training algorithm. The 2DOF-PID controllers used for pH and dissolved oxygen have been tuned by MATLAB’s auto tune feature along with manual tuning. Random training data with input varying from 0 to 100 l/h have been obtained by using NARMA graphical interface. The data samples used for training, validation and testing are 20,000, 10,000 and 10,000 respectively. Random profiles have been used for simulation. The NARMA neuro controller and the 2DOF-PID controllers have shown improvement in rise time, residual error and overshoot. The proposed controllers have been implemented on TMS320 Digital Signal Processing board using code composure studio. Arduino Mega board has been used for input/output interface. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction The multifold growth in the biotechnology industry have taken place in last few decades due to its application in research and development and manufacturing of medical and health care products, agro based products, bio fuels and biodegradable plastics, etc. Biotechnologies have been instrumental in developing and manufacturing products by utilization of organisms and living systems, etc. The products developed or manufactured by biotechnology industry are used for life saving, improving general health standards, and providing alternatives to fossil fuels. It has contributed to the development of mankind in various other fields [1]. Bioprocess are biochemical reactions, which are mix of chemical and biological process. The reaction is carried out in vessels known as bioreactors, which are large stainless steel vessel with capacity ranging from 1000 to 100,000 l. The bioreactor supports a suitable biologically active environment by using appropriate control strategies. Three critical parameters that need to be controlled to follow a desired predefined profile for maximum product yield

∗ Corresponding author. Tel.: +60 49885525; fax: +60 49885510. E-mail address: [email protected] (S.S. Jamuar). http://dx.doi.org/10.1016/j.jprocont.2014.09.012 0959-1524/© 2014 Elsevier Ltd. All rights reserved.

and quality are temperature, pH and dissolved oxygen by affecting the performance of the cells with active enzymes [2,3]. According to Grubecki and Wojcik [4] optimal temperature also results in saving of enzymes which increases productivity. Bioreactor enables control engineers to employ process control techniques to control these parameters in a predefined manner for higher yields. A schematic diagram of a continuous bioreactor is shown in Fig. 1 [5,1]. Fresh stirred medium, which is fed into the bioreactor, is mixed thoroughly and aerated by using a stirrer. The liquid volume in the bioreactor is kept constant at all-times by removing the cell suspensions and desired product simultaneously while the feed is being added in the bioreactor [6]. There are several process parameters that affect the yield. But the variations in process parameters namely temperature, DO and pH greatly affect the yield and quality of product. Hence to increase the yield and to maintain the quality of ethanol these three parameters must follow a precise profile [7,8]. Appropriate sensors are used to measure temperature, dissolved oxygen and pH level. Acid or base are injected through the acid/base valve by means of a peristaltic pump for controlling pH. The dissolved oxygen monitoring controls the speed of stirrer by measuring the concentration of the liquid in the bioreactor. The temperature is controlled by adding a cooling agent in the outside jacket of the bioreactor.

1762

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

Fig. 1. Bioreactor schematic diagram.

Several control strategies like nonlinear model predictive control [9–12], adaptive control [13–15], fuzzy control [16,17] and neural network [18–20] have been employed in the past to replace the conventional PID controllers and to optimize the control of bioprocesses. But due to unavailability of not so accurate model of bioprocess with slow response time, the optimization of desired product yield is still a topic of research [18]. Another reason was unavailability of appropriate sensors and monitoring tools, which have hampered the implementation of such an idea. The main problems encountered with sensors include lack of equal operation dependability in different environments, maintenance of calibration over time, proper interpretation of sensor output, etc. [21]. A description of the sensors used for parameter monitoring is given in the subsequent section.

1.1. Sensors for parameter monitoring Platinum resistance thermometers are the standard thermometers used in bioreactor. These thermometers are usually fitted into thin walled stainless steel pockets which project into the bioreactor. These pockets are also filled with a heat conducting liquid to provide good contact and speed up the instrument’s response. The pH probe has three components namely pH sensor which includes a measuring electrode, a reference electrode and a temperature sensor, a preamplifier and an analyzer or transmitter. The whole pH measurement loop is essentially equivalent to a battery where the positive terminal is the measuring electrode and the negative terminal is the reference electrode. The measuring electrode is sensitive to the hydrogen ion and develops a potential directly proportional to the concentration of hydrogen ions. The reference electrode provides a stable potential against which the measuring electrode is compared. An electrode based oxygen probe is the most common sensor for measuring oxygen dissolved in a liquid. The electrode has a cathode at the center and the anode surrounds the cathode. The anode and the outer membrane are separated by an electrolyte of low density. Oxygen diffusing through the membrane is reduced due to the polarization voltage on the cathode. The result of this electrochemical reaction is an electrical current which is proportional to the partial pressure of the oxygen [22].

Recently nonlinear auto regressive moving average (NARMA) neuro controller have been used to control the nonlinear process by linearizing the nonlinear system dynamics [23], thus optimizing the process. In this paper, nonlinear auto regressive moving average (NARMA) for optimizing the temperature profile has been implemented. This has resulted in faster response time and lower steady state error. Table 1 highlights some of the properties of NARMA controller in contrast with the classical PID controller. Two degree of freedom PID (2DOF-PID) controllers has been designed for the control of pH via regulation of flow of base/acid and for the control of dissolved oxygen via change in stirrer speed. The combination of NARMA controller for temperature profile control and 2DOF-PID for pH and DO has resulted in better performance of bioreactor resulting in higher yield of ethanol as end product. Various profiles have been used for testing. A new profile has been developed for improved production by studying the response of the bioprocess model for various parameter changes by using the test results. The training set for the neural network has been extracted from the developed process model for temperature. Levenberg–Marquardt algorithm [24,25] have been used for training (Appendix I). Two degree of freedom PID controllers have been coarsely optimized with MATLAB’s interactive tuning feature and fine-tuned manually for optimizing the performance of pH and dissolved oxygen controllers. The results show improvement in accuracy and performance of the developed controllers. The controllers developed in Simulink have been converted into machine language code and implemented onto the digital signal processing Table 1 Difference between NARMA and PID. NARMA

PID

There is massive parallel processing. No knowledge of underlying process is required in design of an accurate controller It possesses the ability to perform nonlinear mapping. It has ability to process and determine relations in substantial amount of data with multiple constraints

Serial processing. Complete knowledge of underlying process is required to design a accurate controller. It does not have ability to perform nonlinear mapping. It does not possess the ability to determine relations between various constraints in a set of data.

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

board (DSP) by means of Texas Code Composer Studio. The controllers response have been compared with simulation results for checking the accuracy of the hardware implementation.

dissolved oxygen. In the liquid phase the equilibrium DO concen∗ ) is related to the DO concentration (c ) by tration (cO O2 2

dcO2 dt

2. Process model The block diagram of continuous stirring bioreactor process model is shown in Fig. 2. The process model has been derived from the model described by Nagy [19] and Imtiaz et al. [18]. The model has been implemented as a continuous process and outflow and feed rate has been kept constant. The controller controls three inputs namely flow of cooling agent, flow of acid/base and stirring speed. Temperature, pH level and dissolved oxygen are three outputs of controller. Desired profiles for temperature, pH and dissolved oxygen is achieved by regulating flow of cooling agent, flow of acid/base and stirring speed respectively. The process model is described briefly in the following sections. 2.1. Temperature model Based on energy balance the model for temperature (Tr ) [18,19] can be described as KT AT (Tr − Tag ) F dTr Fe rO2 Hr = i (Tin + 273) − (Tr + 273) + + V V 32r Cheat,r Vr Cheat,r dt (1) where Fi is the substrate flow rate entering the reactor, Tin is the substrate temperature, Fe is the outflow rate, V is the reactor volume, rO2 is the oxygen consumption rate, Hr is the fermentation reaction heat, r is the reaction mass density, Cheat,r is the heat capacity of the mass of the reaction, KT is the coefficient of heat transfer, AT is the heat transfer area, and Tag temperature of the cooling agent in the jacket. The following equation describes temperature of the cooling agent (Tag ) in the jacket dTag KT AT (Tr − Tag ) Fag (T − Tag ) + = Vj in,ag Vj ag Cheat,ag dt

(2)

where Fag is the flow of cooling agent, Tin,ag is the cooling agent temperature entering the jacket, ag is the cooling agent density, Cheat,ag is the cooling agent’s heat capacity and Vj is the cooling jacket volume. 2.2. Dissolved oxygen model [18,19] The dissolved oxygen concentration is varied by changing the stirrer speed and is linked with the equilibrium concentration of

1763





∗ = kl a cO − cO2 − rO2 2

(3)

where kl a is the product of gas specific area and mass-transfer coefficient for oxygen. Rate of oxygen consumption (rO2 ) is given by the following equation rO2 = O2

cO2 1 cX YO2 KO2 + cO2

(4)

where O2 is the maximum specific oxygen consumption rate, YO2 is the amount of oxygen consumed per unit biomass produced, cX is the biomass concentration, and KO2 is the oxygen consumption constant and is represented by the following equation dcX cX Fe = X cX e−Kp cp − cX Ks + cs V dt

(5)

where Ks and Kp represent the constant in the substrate term for growth and the constant of growth inhibition by ethanol respectively. The substrate glucose concentration is represented by cs while cp is the product concentration and X is the maximum specific growth rate. The mass balance for ethanol (product) in which KSI represents the constant in substrate term for ethanol production and p the specific rate of fermentation is given by dcp cs Fe = p cX e−Kp cp − cP KSI + cs V dt

(6)

2.3. pH model The H+ ion concentration defines the pH of the solution [18,19] and is given by pH = −log10 (Q +



Q 2 + 4KW )

(7)

where Q is the difference of hydrogen and hydroxide ions in water which defines the deviation from neutrality and KW is the water constant. The neutralization process of pH is given by the following equation V

dQ = −(FA − FB )Q + CA FA − CB FB dt

(8)

where FA is the flow rate of acid, FB is the flow rate of base, CA is the concentration of acid and CB is the concentration of base. The detail derivations of equations (1)–(8) are given in Nagy [19] and Imtiaz et al. [18].

Fig. 2. Bioreactor process model.

1764

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

3. Controller design

To overcome the problem of slow response, an approximate model of the system has been used [27] in which the dependence on the control input u(k) is linear as given by Eq. (11)

3.1. Nonlinear auto regressive moving average neuro controller The nonlinear auto regressive moving average (NARMA) neuro controller aims at linearizing nonlinear system dynamics by cancelation of the nonlinearities. Eq. (9) gives the standard model to represent discrete nonlinear systems.

yˆ (k + d) = f [y(k), y(k − 1), . . .y(k − n + 1), u(k), u(k − 1), . . .u(k − m + 1)] + g[y(k), y(k − 1), . . .y(k − n + 1), u(k), u(k − 1), . . .u(k − m + 1)] · u(k)

y(k + d) = N[y(k), y(k − 1), y(k − 2). . .y(k − n + 1), u(k + 1)u(k), u(k − 1), u(k − 2). . .u(k − n + 1)]

(11) (9)

where u(k) is the input to the nonlinear system (i.e. the flow of cooling agent for current process), y(k) is the output of the system

u(k) =

yr (k + d) − f [y(k), y(k − 1), . . .y(k − n + 1), u(k), u(k − 1), . . .u(k − m + 1)] g[y(k), y(k − 1), . . .y(k − n + 1), u(k), u(k − 1), . . .u(k − m + 1)]

(i.e. temperature of the current process) and system output y(k + d) is nonlinear function of previous outputs y(i) and previous inputs u(j) where i varies from k − n + 1 to k and j varying between k − n + 1 to k + 1 in which integer d represents the future step for which the output is to be calculated. A nonlinear controller represented by Eq. (10) must be used to make the future output y(k + d) equal to desired output yr (k + d), where n and m are integers. u(k) = G[y(k), y(k − 1), y(k − 2). . .y(k − n + 1), yr (k + d), u(k − 1), u(k − 2). . .u(k − m + 1)]

The above equation is said to be in companion form, where the control input u(k) is separated from the nonlinearity. The control input can be derived from Eq. (11) so that the output y(k + d) follows the reference trajectory yr (k + d) by using Eq. (12). (12)

However above equation cannot be directly implemented because control input is dependent on output at the same time. Hence the model in Eq. (11) is modified to following equation to avoid above conflict: yˆ (k + d) = f [y(k), y(k − 1), . . .y(k − n + 1), u(k), u(k − 1), . . .u(k − m + 1)] + g[y(k), y(k − 1), . . .y(k − n + 1), u(k), u(k − 1), . . .u(k − m + 1)] · u(k + 1)

(13)

(10)

Controller based on dynamic back-propagation algorithm to achieve above was developed, but this had a slow response [26].

The structure of the NARMA neuro controller to implement Eq. (11) is shown in Fig. 3 [28]. The NARMA controller design involves the following steps:

Fig. 3. Structure of the NARMA controller [28].

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

1765

Fig. 4. NARMA graphical interface.

• The collected data is normalized and processed. It is divided into three sets namely training data, validation data and testing data. • The number of input, hidden and output nodes and the number of hidden layers are chosen with an activation function. • The collected and processed data are used for training validation and testing.

were connected to current control input, current plant output, one delayed plant input and 3 delayed plant outputs. The training algorithm which gave the fastest and best result was Levenberg–Marquardt back-propagation algorithm [29]. The tangent sigmoidal function has been used as activation function in design of NARMA controller.

NARMA controller’s training data collection graphical interface in MATLAB is used for training, validation and testing data collection from the plant model as shown in Fig. 4. The maximum and minimum plant inputs are set to 100 and 0 respectively. These values determine the on and off period of the pump for allowing flow of cooling agent. The minimum and maximum sampling interval for collection of two samples, is set to 0.1 and 1 respectively. The data is collected and arranged automatically with above procedure. The optimum size of training, validation and testing data sets was found to be 20,000, 10,000 and 10,000 samples respectively after performing repeated trial and error experiments. The above data samples were obtained keeping in view both the under and over training data values. The structure of the neural network was chosen to have a compromise between network complexity and precision. In our design, there was one hidden layer with seven neurons, one output neuron and 6 input neurons. Six input neurons

3.2. ISA-PID controller ISA-PID controllers have been used for pH and dissolved oxygen control. A classical PID controller is described by the following equation:



U(s) = E(s) Kp +

Ki + Kd s s



(14)

Eq. (14) is known to have a single degree of freedom since the closed loop transfer function can be optimized only for a single task in which for step-response, the performance of controller becomes poor for load disturbance attenuation [30]. ISA-PID controller also known to have two degree of freedom PID (2DOF-PID) is defined by the following equation: U(s) = Kp (bR(s) − Y (s)) +

K  i

s

(R(s) − Y (s)) +

(sKd (cR(s) − Y (s))) (1 + sKd /Kp N) (15)

This equation can be rearranged as U(s) = R(s)F(s) − Y (s)H(s)

(16)

in which H(s) and F(s) are defined as H(s) = Kp + F(s) = Kp b + Fig. 5. Block diagram of two degree of freedom controller [31].

sKd Ki + s (1 + sKd /Kp N) Ki sKd c + s (1 + sKd /Kp N)

(17) (18)

1766

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

Fig. 6. Simulink diagram of the process implementation.

In Eq. (15), Kp is the proportional gain, Ki is the integral gain, Kd is the derivative gain, E(s) is the control error, U(s) is the control signal, Y(s) is the process output, N is the derivative filter, R(s) is the set point, and c and b are weights that influence the set-point response. The transfer functions F(s) and H(s) can be optimized for two different tasks. One can be used to maintain good output performance and the other could be used to maintain good regulator performance [31]. The block diagram of two degree of freedom PID controller, where G(s) is the transfer function and D(s) is the disturbance signal, is shown in Fig. 5. Both H(s) and F(s) have common terms namely Kp , Ki , Kd and N, therefore change in parameters, which affect regulator performance also affect the output performance. However the weights c and b in Eq. (18) can be selected appropriately to enhance one parameter without affecting the other parameter [31]. MATLAB’s auto-tuning feature was used for two degree of freedom PID controller and to ensure precision in following a desired profile for optimizing regulator and servo controller. The precision manual tuning has also been carried out to optimize the output. For DO controller 250, −1.6 and 123 were the values chosen for K. The values of Kp , Kd and Ki for the pH were chosen to be 25, −1 and 2 respectively.

4. Implementation in MATLAB Fig. 6 shows the implementation diagram of the process in Simulink. Block A represents reference profiles for temperature, pH and DO that the controller must follow with precision. The temperature, pH and DO controllers are shown in block B. By changing the flow of cooling agent (Fag ), the temperature (T) is controlled indirectly by using NARMA neuro controller. The pH level is controlled by changing the flow of base/acid (FB ). By changing the stirring speed of the impellers, the dissolved oxygen (DO) concentration is controlled. The chosen profiles are made to change randomly to verify the operation of controllers for any desired profiles in our simulation and the response and performance of the controllers is observed.

5. Stability analysis Bounded input bounded output stability tests were performed on the designed controllers over an operational range of input and output values. The operational range for each controller was determined by literature review. For the NARMA controller the temperature profile was varied randomly in temperature range of 26.2–30.3 ◦ C. This range has been selected because yeast cells acclimate to a temperature range of 27–30 ◦ C for optimum production of ethanol [32]. Fig. 7 shows the results of the experiment where it can be seen that the output of the controller remains bounded with the range of 0–100 l/h for an input range of 25.5–30.2 ◦ C which indicates that the controller is stable over this range. It was suggested by Singh and Trivedi [33] that pH range 5–6 is optimal for the production of ethanol. Hence the stability of the 2DOF-PID controller was tested over the range of 4–7.6 pH. Fig. 8 shows the input and output for the pH 2DOF-PID controller where it can be seen that the output of the pH controller remains bounded within the range of −6.5 to 8.2 l/h. The positive and negative outputs represent the flow of base and acid respectively. Hence the pH controller is deemed to be stable over this range. For optimum production of ethanol the minimum dissolved oxygen concentration is 3 g/l [34] but if the dissolved oxygen concentration is too high it may lead to hyperoxic condition resulting in reduction of cell growth. Therefore the dissolved oxygen concentration must be maintained just above 3 for an optimal cell growth. The stability was thus tested in the range of 3.3–3.9 g/l. Fig. 9 shows the input and output variations of 2DOF-PID DO controller. The controller remains bounded within the range of 19–50 rpm for the range of input values, which shows that the controller is stable over this range of values.

6. Controller implementation The controllers developed in Simulink were converted into machine language code and implemented onto the digital signal processing board (DSP) by means of Texas Code Composer Studio. This has been done to check the feasibility of implementation. The

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

1767

Fig. 7. NARMA controller input and output.

Fig. 8. pH controller input and output.

Simulink representation of the controllers that were implemented on the DSP is shown in Fig. 10. The parameters of the bioreactor model is given in Appendix II. In Fig. 10, blocks labeled 1, 2 and 3 are the inputs to the DSP. Block 1 is the input where temperature feedback from the bioreactor. Block 2 is the pH feedback and block 3 is the input where the concentration of dissolved oxygen is connected. The DSP board was tested against the bioreactor mathematical model in

Simulink via Arduino board functioning as an input/output device and not on an actual bioreactor. Block diagram representation of the implementation is shown in Fig. 11. The computer is interfaced with the Arduino board. Since the voltage levels between Arduino board and the DSP board are different, a voltage level converter has been used as interface between two. Fig. 12 shows the details of voltage interface block between the DSP board and Arduino.

Fig. 9. DO controller input and output.

1768

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

Fig. 10. Simulink diagram of implemented controllers.

As seen in the figures 12 pins from the DSP and Arduino are connected to each other via the interface voltage level converter. The maximum input and output voltage of DSP is +3 V, voltage divider circuits are used to reduce the peak voltage coming from the Arduino which is 5.0–2.7 V. The block representation of the voltage divider circuit is shown in Fig. 13 in which a 500  resistor and a Zener diode with breakdown voltage of 2.7 V is used. An open collector buffer integrated circuit with supply voltage of +5 V is used to change the voltage to 5 V for the Arduino input as shown in Fig. 14. The DSP which represents the set of controllers has 3 input pins and 3 output pins. The input pins are temperature, pH and dissolved oxygen levels of the reactor and the output pins are rate of cooling agent flow, flow of acid/base and stirring speed. Similarly out of the 6 pins used on the Arduino 3 input pins are rate of cooling agent flow, flow of acid/base and stirring speed and the 3 output pins are temperature, pH and dissolved oxygen levels of the reactor.

7. Simulation results and discussion This section discusses the simulation results obtained for different profile controllers. For the purpose of comparison an antiwindup PID (AWU-PID) described in [35] has also been simulated. 7.1. Temperature profile testing The profile given in Table 2 has been chosen to test profile following capability of the NARMA neuro controller. Set point change in the profile occurs after every 50 h and contains both small and large changes so as to evaluate not only the delay in the controller’s response and rise time but residual error and overshoot as well. The combined results of the anti-windup PID and the NARMA are shown in Fig. 15 and the corresponding changes in flow rate have been shown in Fig. 16. From Fig. 15 we note that the NARMA

Fig. 11. Implementation block diagram.

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

1769

Fig. 12. Voltage interface specifics.

controller’s response to change in the set point is faster than the anti-windup PID. Due to the frequent change in set point the antiwindup PID is unable to drive the system temperature very near to the set point and as such we observe a large rise time and residual error. The NARMA controller on the other hand is able to drive the

system temperature to the set point level with greater accuracy and significantly reduced rise time. The rise time for the NARMA controller is significantly less than the PID controller. The level of overshoot in the PID controller is proportional to the change in set point level as such at a large change we can observe substantial

1770

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777 Table 2 Temperature profile.

Fig. 13. Voltage divider circuit.

Time (h)

Temperature (◦ C)

0–50 50–100 100–150 150–200 200–250 250–300 300–350 350–400 400–450 450–500

31 28.5 29.5 28 29 27 30 28.5 30.5 28

overshoot. The overshoot is more evident at 0–50 h and 100 h. Overshoot for the NARMA controller however remains negligible throughout the entire simulation. These remarks can be explained by observing the flow rate of cooling agent shown in Fig. 16. It can be seen that the flow rate of the anti-windup PID changes slowly at all times including the instances where the temperature has large change. This shows that the anti-windup PID has undercompensated the large change and fails to differentiate between a small change and a large change. In contrast to this the NARMA has shown much better compensation for large changes as well as small changes. It can be concluded from these experiment that the NARMA controller shows improvement in overshoot, residual error, delay as well as rise time when compared to the anti-windup PID while following the desired profile. 7.2. pH profile testing

Fig. 14. Open collector buffer integrated circuit.

The designed 2DOF-PID was simulated along with the AWUPID for comparison. The results of the simulation are shown in Figs. 17 and 18. The target profile which is designed to have both large and small changes frequently to better evaluate the performance is given in Table 3. The system output versus the desired profile is shown in Fig. 17 and the flow of acid/base is shown in Fig. 18. The negative part in Fig. 18 represents flow of acid and the positive part represents flow of base. From the results shown in Fig. 17 it can be observed that the 2DOF-PID is able to accurately maintain the pH at the target level at all points while the AWU-PID not only shows great delay in rise time but also residual error. A slight overshoot is observed in the 2DOFPID’s control when a major change occurs in the desired profile towards increase in basicity at 50 h but the settling time however

Fig. 15. Anti-windup PID and NARMA combined results for temperature control.

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

1771

Fig. 16. Comparison of flow rates for temperature control.

Fig. 17. Comparison of 2DOF-PID and AWU-PID for pH control.

is seen to be very low. From Fig. 18 we note that for every change towards acidity in desired profile the 2DOF-PID increases the flow of acid which is followed by a slight flow of base and vice versa. This is done to reduce settling time as the pH of the medium can

be effected immediately by the flow acid/base since both acid and base are administered directly into the medium unlike temperature where cooling agent is regulated in an outer jacket and heat is reduced by conduction.

Fig. 18. Flow rate comparison of 2DOF-PID and AWU-PID for pH control.

1772

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

Fig. 19. Comparison of 2DOF-PID and AWU-PID control results for DO control.

Fig. 20. Comparison of stirring speed results of 2DOF-PID and AWU-PID for DO control.

Residual error and rise time can be explained by observing the manipulated variable, i.e. flow of acid/base in Fig. 18. Every change in the desired level of pH is responded by a rapid change in the flow rate of acid/base in case of 2DOF-PID whereas the

response by AWU-PID is noted to be slow. The flow rate change by 2DOF-PID itself can be seen to be higher than AWU-PID when a change occurs in the desired profile. This observation suggests that the AWU-PID has undercompensated the change in pH level

Fig. 21. Temperature implementation results.

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

1773

Fig. 22. pH implementation results.

Fig. 23. DO implementation results.

Table 3 pH profile. Time (h)

pH

0–50 50–100 100–150 150–200 200–250 250–300 300–350 350–400 400–450 450–500

4.5 6.5 6 7.5 3.5 2.5 3 5 4.5 3

resulting in increase in rise time and the presence of residual error. In conclusion the level of overshoot is directly proportional to the amount of change in the target profile but it lies within the tolerance level of ±0.5. Apart from this improvement has been seen in rise time, settling time and residual error. 7.3. DO simulation results Dissolved oxygen concentration plays an important role in the growth of cells. At high cell concentration the rate of oxygen

consumption may exceed the rate of oxygen supply, leading to oxygen limitations. In such conditions the growth rate varies directly with the concentration of dissolved oxygen and as such the dissolved oxygen concentration is maintained above a critical value at all times and the designed profile generally does not contain frequent changes. However to test the profile following capability of the 2DOF-PID controller the designed profile contains several changes in desired value with the lowest value being 2.7 g/l and the highest being 3.6 g/l. The complete profile is shown in Table 4.

Table 4 Dissolved oxygen. Time (h)

DO (g/l)

0–50 50–100 100–150 150–200 200–250 250–300 300–350 350–400 400–450 450–500

3.4 3.1 3.3 3.3 3.6 3.1 3.2 3.5 2.7 2.8

1774

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

Fig. 24. Magnified plot of temperature at set point changes.

Dissolved oxygen control by 2DOF-PID was simulated along with the anti-windup PID for comparison. The results for DO concentration of the plant and the corresponding stirring speed have been plotted in Figs. 19 and 20 respectively. It can be observed from the figures that the 2DOF-PID controller is able to accurately maintain the DO concentration to the set point with reduced rise time when compared with the anti-windup PID controller. This demonstrates the ability of the 2DOF-PID controller not only to assess large and small changes in the desired profile but its quick response in adjusting stirrer speed. In order to maintain the DO level the stirrer speed must counter the oxygen consumption rate and it can be observed in Fig. 20 that even though the desired oxygen level is same the stirred speed is slightly increased or decreased to maintain balance with the oxygen consumption rate. The residual error for the anti-windup PID controller seen in Fig. 19 arises from the inability of the controller to effectively counter the oxygen consumption rate which in the case of the 2DOF-PID is seen to be negligible. In conclusion it can be said that the AWU-PID fails to respond as quickly and as accurately as the 2DOF-PID controller and the designed controller shows improvement in response time and reduced residual error.

8. Implementation results and discussion The designed NARMA and 2DOF-PID controllers were implemented on a DSP board and the results for temperature, pH and DO control shown in Figs. 21–23 are for optimized profiles

for maximum product output and have been plotted along with the simulation results. Table 5 shows the optimized profile and Figs. 24–26 show magnified plots of temperature, pH and DO respectively at set point changes. It is observed that the control results for the parameter values show a ripple around the target values with an average error of ±0.01–0.05. In Fig. 21 it can be seen that the residual error is of the order of ±0.3 which in comparison to the simulation results shown in the same figure is quite high but still within the tolerance limit of ±0.5. The rise time remains the same as that of the simulation but the overshoot is higher at step times 70 h and 300 h. This can be clearly seen in Fig. 24 in which shows a magnified plot of set point changes. In Fig. 22 it can be noticed that the rise time has been effected at time 300 and 400 h and apart from this the rise time at the start remains the same as that of the simulation results. In Fig. 23 no overshoot is seen but a delay in rise is observed at all step changes. Apart from time 70 h to 120 h the residual error also follows the same pattern as the simulation.

Table 5 Optimized profile. Time

Temperature

pH

DO

0–75 75–200 200–300 300–400 400–500

29.5 28 29 28.5 30

5 5 5 5.6 6

4 4 3.8 3.7 3.4

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

1775

Fig. 25. Magnified plot of pH at set point changes.

Fig. 26. Magnified plot of DO at set point changes.

Changes in parameter values at step changes in target profile have been zoomed in and plotted in Figs. 25 and 26. In these figures it can be clearly seen that the residual error, rise time and overshoot of the implemented controllers is slightly higher when matched to the simulation results. This can be attributed to the noise and interference caused due to the connections, wires, etc. In spite of these errors the implementation results remain under the tolerance level of ±0.5 and as such can still be applied to a bioreactor since the time frame of the complete bioreactor run is estimated to be around 500 h. Hence it can be concluded that the implementation has not only successfully confirmed the feasibility of using DSP processor board to be capable of running both the NARMA neuro controller as well as the 2DOF-PID controllers but also proven that

the programmed DSP board can autonomously run an entire 500 h bioprocess efficiently.

9. Conclusion The biochemical reactor model for fermentation by Saccharomyces yeast was developed with temperature, pH and DO as critical output variables affecting the process dynamics. The kinetic model of the process, e.g. dependence of kinetic parameters of temperature, the mass transfer of oxygen and the influence of temperature and ionic strength on the mass transfer coefficient have been modified to cater for the nonlinear process dynamics.

1776

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

The biochemical reactor model was implemented in MATLAB and Simulink. The nonlinear auto regressive moving average (NARMA) neuro controller has been trained to act as controller for temperature control and the pH and DO parameters have been controlled by using two degree of freedom PID (2DOF-PID). The controller performances have been compared with the AWU-PID controllers. It has been observed that the AWU-PID fails to respond quickly and accurately compared to the NARMA and the 2DOF-PID controllers. The NARMA and the 2DOF-PID controllers show improvement in response time, reduced residual error and delay time. The behavior of the manipulated variable has also been discussed with respect to the changes in the desired profiles. The proposed controllers were tested for rapid set point changes in the desired profiles. The accuracy of these controllers was studied by predefined profile temperature. The performance of the controllers showed improvements in rise time, overshoot and residual error. The controllers have also been implemented on a digital signal processing board. Results of the implementation indicated the response of the controllers was satisfactory in real time which confirmed the applicability the DSP board for controller implementation and the feasibility of the overall design to run a bioreactor process autonomously.

Acknowledgement This research is financially supported by Universiti of Malaya, Ministry of Higher Education High Impact Research (UM.C/HIR/MOHE/ENG/20).

Appendix I. A.1. Levenberg–Marquardt algorithm Least squares is a standard approach to find the solution of a system of equations that minimizes the sum of squares of the errors made in the results of every single equation. In general the problem reduces to finding x that is a local minimizer to a function that is a sum of squares subjected to some constraints represented by f(x). The function f(x) is minimized that is a sum of squares in the least-squares problem



2

minx f (x) = F(x) = 2

 i

Fi2 (x)

(B1)

Such problems arise frequently when fitting model functions to data in most practical applications, e.g. estimation of nonlinear parameters. Such problems also arise in control systems when the output, y(x,t) must follow some predefined profile ϕ(t), where x is any vector and t is any scalar. Such a type of problem can be described as



where the weights of the quadrature formula are y¯ and ϕ. ¯ In this scenario F(x) can be described by Eq. (B4)

⎡

¯ t1 ) − ϕ(t ¯ 1) y(x,

⎢ ¯ t2 ) − ϕ(t ¯ 2) ⎢ y(x, ⎢ F(x) = ⎢ ¯ t3 ) − ϕ(t ¯ 3) ⎢ y(x, ⎢ ⎣ ...

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

¯ y(x, tm ) − ϕ(t ¯ m) The residual F(x) in such problems is generally very small near the optimum because generally realistic and achievable target profile are. The function in the left hand side (LS) can be minimized by means of any minimization strategy but to enhance the efficiency of the solving procedure a few characteristics of the problem can be explored. The Hessian matrix and gradient of LS have a special structure in which the main diagonal does not contain any mixed partial derivatives. Representing J(x) as the m-by-n Jacobian matrix of F(x) and G(x) as the gradient vector of f(x), H(x) as the Hessian matrix of f(x) and Hi (x) as the Hessian matrix of each Fi (x) we will get G(x) = 2J(x)T F(x)

(B5)

H(x) = 2J(x)T J(x) + 2Q (x)

(B6)

in which Q (x) =

m i=1

Fi (x) · Hi (x)

minx=Rn

(y(x, t) − ϕ(t))2 dt

(B2)

t1

here both y(x,t) and ϕ(t) are scalars. By means of a quadrature formula the integral is discretized and Eq. (B2) can be rewritten as

minx=Rn f (x) =

m i=1

¯ (y(x, ti ) − ϕ(t ¯ i ))2

(B3)

(B7)

When the residual F(x) moves to zero as x comes closer to the solution, Q(x) also moves toward zero. A very effective method will be to employ Gauss–Newton direction for the optimization procedure when F(x) is very small near solution. A search direction, dk is calculated at each major iteration k, in the Gauss–Newton method which comes from solving the linear least-squares problem:



2

minx=Rn J(xk ) − F(xk )

2

(B8)

When Q(x) can be ignored, the direction calculated from this procedure is equal to the Gauss–Newton direction. To ensure that the function f(x) in each iteration decreases, dk is used as part of a line search technique. When the second-order term Q(x) is significantly large, the Gauss–Newton method encounters problems. Levenberg–Marquardt method is a method that overcomes this problem. A solution of the linear set of equation (B9) is used as a search direction (J(xk )T J(xk ) + k I)dk = −J(xk )T F(xk )

(B9)

otherwise Eq. (B10) (J(xk )T J(xk ) + k diag(J(xk )T J(xk )))dk = −J(xk )T F(xk )

t2

(B4)

(B10)

where k controls both the direction and magnitude of dk . When k is zero, the direction dk is identical to that of the Gauss–Newton method. When k moves toward infinity, dk moves toward the sharpest direction of descent. Hence for a large k , the term F(xk + dk ) < F(xk ) is valid. k can hence be varied to make the direction toward descent even if second-order terms, which confine the effectiveness of the Gauss–Newton method, are present. Hence it can be concluded that the Levenberg–Marquardt method uses a search direction which is a mix of the steepest descent direction and the Gauss–Newton direction [28].

U. Imtiaz et al. / Journal of Process Control 24 (2014) 1761–1777

Appendix II. Parameters of bioreactor model Parameter

Value

Unit

A1 A2 AT Cheat,ag Cheat,r Ea1 Ea2 HCa HCO3 HCl HH HMg HNa HOH KO2 Kp KP1 Ks KSI KT KW mCaCO3 mMgCl2 mNaCl R RSP RSX Vj YO2 Hr O2 p ag r

9.5 × 108 2.55 × 1033 1 4.18 4.18 55,000 220,000 −0.303 0.485 0.844 −0.774 −0.314 −0.550 0.941 8.86 0.139 0.070 1.030 1.680 3 × 105 10−14 100 100 500 8.31 0.435 0.607 50 0.970 518 0.5 1.790 1000 1080

– – m2 J g−1 K−1 J g−1 K−1 J/mol J/mol – – – – – – – mg/l g/l g/l g/l g/l J h−1 m−2 K−1 (mol/l)2 g g g J mol−1 K−1 – – l mg/mg J/mol O2 h−1 h−1 g/l g/l

References [1] M. Henson, Biochemical reactor modelling and control: exploiting cellular biology to manufacture high-value products, IEEE Control Syst. 26 (2006) 54–62. [2] I. Grubecki, Analytical determination of the optimal temperature profiles for the reactions with parallel deactivation of enzyme encapsulated inside microorganism cells, Chem. Biochem. Eng. Q. 26 (1) (2012) 31–43. [3] I. Gruubecki, M. Wojcik, Analytical determination of the optimal temperature profiles for the reaction occurring in the presence of microorganisms cells, Biochem. Eng. J. 39 (2008) 362–368. [4] I. Grubecki, M. Wojcik, How much of enzyme can be saved in the process with the optimal temperature control? J. Food Eng. 116 (2013) 255–259. [5] J. Williams, Keys to bioreactor selections, Chem. Eng. Prog. March (2002) 34–41. [6] M. Shuler, F. Kargi, Bioprocess Engineering: Basic Concepts, Prentice Hall, Upper Saddle River, NJ, 2002. [7] A. Meszaros, A. Andrasik, P. Mizsey, Z. Fonyo, V. Illeova, Computer control of pH and DO in a laboratory fermenter using neural network technique, Bioproc. Biosyst. Eng. 26 (2004) 331–340. [8] W. Hochfeld, Producing Biomolecular Substances with Fermenters, Bioreactors and Biomolecular Substances, CRC Press, Florida, USA, 2006.

1777

[9] A. Costa, L. Meleiro, R. Maciel Filho, Non-linear predictive control of an extractive alcoholic fermentation process, Process Biochem. 38 (2002) 743–750. [10] J. Dowd, K. Kwok, J. Piret, Predictive modelling and loose-loop control for perfusion bioreactors, Biochem. Eng. J. 9 (2001) 1–9. [11] B. Foss, T. Johansen, A. Srenson, Nonlinear predictive control using local models applied to a batch fermentation process, Control Eng. Pract. 3 (1995) 389–396. [12] K. Preub, M. Lann, M. Cabassud, G. Anne-Archard, Implementation procedure of an advanced supervisory and control strategy in the pharmaceutical industry, Control Eng. Pract. 11 (2003) 1449–1458. [13] M. Guay, D. Dochain, M. Perrier, Adaptive extremum seeking control of a continuous stirred tank bioreactor with unknown growth kinetics, Automatica 40 (2004) 881–888. [14] S. Valentinnotti, B. Srinicasan, U. Holmberg, D. Bonvin, C. Cannizarro, M. Rhiel, U. Stockat, Optimal operation of a fed batch fermentations via adaptive control of overflow metabolite, Control Eng. Pract. 11 (2003) 665–674. [15] I. Smets, J. Claes, J. November, G. Bastin, J.F.V. Impe, Optimal adaptive control of (bio)chemical reactors: past, present and future, J. Process Control 14 (2004) 795–805. [16] K. Nakano, R. Katsu, K. Tada, M. Matsumura, Production of highly concentrated xylitol by Candida magnoliae under a micro aerobic condition maintained by simple fuzzy control, J. Biosci. Bioeng. 89 (2000) 372–376. [17] C. Karakuzu, M. Turker, S. Ozturk, Modelling, on line state estimation and fuzzy control of production scale fed-batch Baker’s yeast fermentation, Control Eng. Pract. 14 (2006) 959–974. [18] U. Imtiaz, A. Assadzadeh, S. Jamuar, J. Sahu, Bioreactor temperature profile controller using inverse neural network (INN) for production of ethanol, J. Process Control 23 (2013) 731–742. [19] Z.K. Nagy, Model based control of yeast fermentation bioreactor using optimally designed artificial neural networks, Chem. Eng. J. 127 (2007) 95–109. [20] E. Sivaraman, S. Aruselvi, Neuro modelling and control strategies for a pH process, Int. J. Eng. Sci. Technol. 3 (2011) 194–203. [21] J. Alford, Bioprocess control: advances and challenges, Comput. Chem. Eng. 30 (2006) 1464–1475. [22] A. Scragg, Bioreactors in Biotechnology: A Practical Approach, Horwood, New York, 1991. [23] F. Mjalli, S. Al-Asheh, Neural-networks-based feedback linearization versus model predictive control of continuous alcoholic fermentation process, Chem. Eng. Technol. 28 (2005) 1191–1200. [24] K. Levenberg, A method for the solution of certain non-linear problems in least squares, Q. J. Appl. Math. 2 (1944) 164–168. [25] D. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, SIAM J. Appl. Math 11 (1963) 431–441. [26] K.S. Narendra, K. Parthasarathy, Learning automata approach to hierarchical multiobjective analysis, IEEE Trans. Syst. Man Cybern. 21 (1991) 263–272. [27] K. Narendra, S. Mukhopadhyay, Adaptive control using neural networks and approximate models, IEEE Trans. Neural Netw. 8 (1997) 475–485. [28] M.H. Beale, M. Hagan, H.B. Demuth, Neural Network ToolboxTM User’s Guide, The Math Works Inc., Natick, 2012. [29] J. Nocedal, S. Wright, Numerical Optimization, Springer, New York, 1999. [30] R. Vilanova, V. Alfaro, O. Arrieta, Analytical robust tuning approach for twodegree-of-freedom PI/PID controllers, Eng. Lett. 19 (3) (2011) 204–214. [31] D. Czarkowski, T. Mahony, Intelligent controller design based on gain and phase margin specifications, in: Proceedings of Irish Signals and Systems Conference, IET Digital Libraries, Belfast, 2004, pp. 381–386. [32] S.L. Tai, P. Daran-Lapujade, M.C. Walsh, J.T. Pronk, J.-M. Daran, Acclimation of Saccharomyces cerevisiae to low temperature: a chemostat-based transcriptome analysis, Mol. Biol. Cell 18 (2007) 5100–5112. [33] D. Singh, R. Trivedi, Acid and alkaline pretreatment of lignocellulosic biomass to produce ethanol as a biofuel, Int. J. ChemTech Res. 5 (2013) 727–734. [34] A. Gharahvaran, Design of profile controller for biochemical reactor, Universiti Putra Malaysia, Serdang, Selangor, Malaysia, 2010. [35] K. Astrom, T. Hagglund, Advanced PID Control, ISA, North Carolina, 2005.