Acta Polytechnica Hungarica Vol. 12, No. 1, 2015 Biomass-fired Boiler Control Using Simulated Annealing Optimized Improved Varela Immune Controller Alexander Hošovský Department of Mathematics, Informatics and Cybernetics, Faculty of Manufacturing Technologies with seat in Prešov, Technical University of Košice, Bayerova 1, 08001 Prešov, Slovakia, e-mail: [email protected] Abstract: Water temperature control of biomass-fired boilers represents a process with extreme delay which makes it quite difficult to stabilize it using conventional PID controller tuned according to commonly used PID tuning rules. It is proposed here to use a controller based on Varela-Countinho second generation immune network which was shown to have very good anti-delay capabilities. Since there are currently no tuning rules for this type of controller, simulated annealing algorithm is used for optimizing the set of controller parameters to achieve good performance according to IAE criterion. The resulting controller is shown to offer stable performance even for such a long time-delay which is also robust for up to 30% variations in system time constant compared to a nominal case. Keywords: B-cells; antibody; temperature; optimization; controller 1 Introduction An effective control of biomass combustion and biomass-fired boilers in general in regard to meeting ever stricter requirements on emission levels remains a challenging task [1], [2], [3]. In works [2] and [3] the aspects of biomass-fired boiler operation in terms of carbon monoxide and oxygen concentration relationship are analyzed. Despite the complexities of processes associated with biomass combustion, the commonest controller type applied to this task is still a conventional PID controller [4], [5], [6]. In work [5] the ways of improving the performance of current control techniques are inspected while in [6] Smith predictor is applied to compensate for the performance deterioration due to long time delay of a boiler. The prevalence of PID control can almost certainly be attributed not that much to PID’s satisfying performance in that task but more probably to its ease of use, multitudes of tuning rules and good comprehensibility. Nevertheless, ambitions for improving the performance of biomass-fired boiler control using more advanced techniques are still present, e.g. in works [5], [7], [8]. – 23 – A. Hošovský Biomass-fired Boiler Control Using SA Optimized Varela Immune Controller The works [7] and [8] present a preliminary analysis of possible intelligent techniques applied to the control of biomass-fired boilers. In [6] a control design model for medium scale biomass-fired boiler was developed. Precisely, this kind of plant should be treated as multivariable system with several inputs and several outputs that are cross-linked. However, based on the measured data it was shown that water temperature control could be treated as single input single output (SISO) system with fuel feed as an input variable and water temperature as an output variable. An identified time-delay for this system is 480 s (for increasing fuel feed), which makes it difficult to be stabilized using standard tuning rules for PID controllers. A method of time-delay compensation using Smith predictor (e.g. [9]) is proposed – this, however, requires the model to be equivalent to the plant, which might not always be fulfiled in practice and this may lead to deterioration of controller’s performance. It is therefore proposed to use Varela immune controller which was shown [10], [11] to have very good anti-delay capabilities. Since there are no tuning rules for this kind of controller and the effect of its parameters on controller’s performance is less intuitive than in case of PID controller, it is also proposed to optimize its parameters using simulated annealing. It is shown that this combination provides a stable and robust performance even for a system with such long time-delay. 2 Used Methods 2.1 Improved Varela Immune Controller Varela immune controller belongs to the group of nonlinear controllers which make use of some of natural immune system paradigms. The theory of immune controllers was extensively studied by Chinese researchers [10], [11], [12]. Varela immune controller is based on the idea of second generation immune networks proposed by Varela and Countinho in [13]. This model takes into account the interactions between B-cells and antibodies present as free soluble molecules. The central idea lies in the relationship between B-cell maturation (only fully mature B-cell can produce antibodies) and proliferation probability and cumulative receptor occupancy (or sensitivity). This characteristic has the form of two mutually shifted bell functions from which the basic fact that B-cells are activated only at intermediate receptor occupancy can be observed (Fig. 1). – 24 – Acta Polytechnica Hungarica Vol. 12, No. 1, 2015 Figure 1 Probability functions for B-cell maturation and proliferation in relation to network sensitivity The dynamics of second generation immune network can be described by two differential equations [13], [14]: df i f (k k )k bMat() (1) i 1 2 i 3 i i dt db i b (k Prolk )Meta[i] (2) i 4 i 5 dt where k – antibody natural death rate, k – antibody interaction death rate, k – 1 2 3 antibody production rate by fully matured B-cells, k – B-cell proliferation rate, k 4 5 – B-cell death rate, f – free soluble antibody with i idiotype, b – B-cell with i i i idiotype, Meta[i] – additional cells from the resting pool, Mat() – maturation i probability function, Prol() – proliferation probability function. The symbol  is i i the sensitivity of network to i idiotype defined as: N (t)m f (3) i ij j j1 where m is affinity between i and j idiotype. ij It is clear from (1) and (2) that the dynamics of free antibodies as well as B-cells with specific idiotype are related to the network sensitivity to this idiotype (i.e. receptor occupancy). It is possible to construct a controller which uses the foundations of second generation immune networks on the condition that some special assumptions relevant to control theory are made. The theory of second generation immune networks describes the activity of immune system in the absence of antigens [10]. The adaptation of immune system theory to automatic control, however, uses the notion of antigen (as a foreign element that needs to be eliminated) to represent a control error. It is therefore necessary to incorporate the effects of antigen presence into the immune controller model. According to [15], the rate of antigens can be given by the following equation: a k a k f a (4) i ag i e i i – 25 – A. Hošovský Biomass-fired Boiler Control Using SA Optimized Varela Immune Controller where k – multiplication rate of antigens, k – elimination rate of antigens, a – ag e i antigen with i idiotype. Moreover, the Meta term in equation (2) might be modified in order to account for B-cells supplied from bone marrow modeled by constant k and B-cells proliferation due to antigen presence modeled by factor bm k a. The final model used for constructing Varela controller can be summarized bag i as follows: da i k a k fa ag i e i i dt df i f (k k)k bMat() i 1 2 i 3 i i dt (5) db i b(k Prolk )k k a i 4 i 5 bm bag i dt N (t)m f i ij j j1 The four equations (5) constitute the biological model of immune network dynamics according to Varela-Countinho with addition of antigen dynamics equation (first equation in (5)). This biological model has to be transformed into a model that would use control variables needed for its function as a controller. The interactions in (5) are shown in Fig. 2 (with + sign denoting stimulation effect and – sign denoting suppression effect) together with certain analogies between biological variables and control variables. Figure 2 Interactions in biological model of Varela immune network with inclusion of antigen effect (left) and some analogies between biological and control model (right) According to [12], two assumptions are made in regard to the implementation of Varela controller when compared to its biological model shown in (5) – 1) Antigen is capable of replication in human body, which is modeled by the first term in the first equation of biological model. As was mentioned above, antigen is treated as a control error in immune controllers but there is no analogy for its replication in control systems and thus it is not considered in the development of Varela immune controller. – 26 – Acta Polytechnica Hungarica Vol. 12, No. 1, 2015 2) Also, in contrast to B cell and antibody concentration, controller effort and control error can assume negative values. Maturation and proliferation functions were thus modified in the following way: Muk (e1u e2u )sign(u) M (6) Puk (e1u e2u )sign(u) P where ,  – negative constants and  < , and k , k – positive constants. 1 2 2 1 M P Also the network sensitivity was replaced with u in maturation function and with uin proliferation function. The shape of this function is depicted in Fig.3 for normalized values of respective variables. Figure 3 Maturation and proliferation function in Varela immune controller Using the abovementioned assumptions and analogies shown in Fig. 2 we may rewrite the equations (5) to the following form: de k f(t)e(t) e dt df k f(t)k M(u)u(t) (7) f 3 dt du u(t)(k Prol(u(t))k )k k e(t) 4 5 bm bag dt where in the first equation the self-replication term of antigen (now control error) has been removed and the term (k + k) has been replaced with k. It can also be 1 2 f seen that the antibody concentration f has been preserved also in Varela immune controller model as there is no analogy for this variable in controller design – 27 – A. Hošovský Biomass-fired Boiler Control Using SA Optimized Varela Immune Controller (actually f is one of the three state variables of Varela controller, other two being u and du/dt). In order to implement a practical controller, it was proposed in [11] to further modify equations in (7) by taking time derivative of the last equation, neglecting the second term in derivative of u(t)k Prol(u(t))and substituting the 4 control error rate equation into this equation to get: df k f(t)k M(u)u(t) f 3 dt (8) d2u u(t)k Prol(u(t))k u(t)k f(t)e(t) dt2 4 5 I where k = kk . This is the final model of Varela immune controller used for I e bAg experimentation. 2.2 Simulated Annealing Simulated annealing is an optimization technique based on the idea of annealing process in metalurgy. This process is described as a way of obtaining low energy state of a solid in heat bath [16]. First a solid is heated to a high temperature at which the atoms in lattice are arranged randomly and then it is cooled in a controlled manner so that a regular crystal lattice is formed – this state corresponds to a state with minimal energy. If the temperature is decreased to quickly, a state with higher than minimal energy (corresponding to a polycrystalline structure) is attained. Inventors of simulated annealing algorithm used an analogy of this process for optimization technique where the state of a solid corresponds to a problem solution, energy associated with a current state corresponds to a fitness function value and the temperature represents a control parameter. The points within the vicinity of current solution are generated using the following function [17]: x2 y(x) 1 gX ,TX T ;y(x) exp 2 (9) C c c y(x) 2 2 where X – n-dimensional current point vector, T – current temperature vector, c c y(x) – pseudorandom generator function, x – a number taken from the default random stream and . - 2-norm. In this implementation, the perturbation of 2 current solution depends on current temperature (represented as a vector, i.e. can be different for every element of solution vector). The acceptance function decides whether the generated solution is accepted or not (if the new solution decreases the fitness function value it is certainly accepted while if it increases the fitness function value it can still be accepted with a probability given below): – 28 – Acta Polytechnica Hungarica Vol. 12, No. 1, 2015 f  f  f 0,P(X X )1 n c n c 1 f  f  f 0,P(X X ) (10) n c n c f 1expmax(Tc) where f – cost function for a new point, f – cost function for a current point, P(.) n c – probability function. The performance of simulated annealing is strongly dependent on the selection of cooling schedule which was chosen to be of exponential type. The schedule could be expressed in this form [17]: T T 0.95k (11) n c where T – new temperature vector, k – annealing parameter vector. The n implementation of simulated annealing algorithm had also capability of reannealing (repeated increase in temperature vector values according to the sensitivity of fitness function to the changes in respective parameters) expressed in the form of sensitivity function [17]: f(X ) T max(S) SR c ,k ln 0 (12) X T S c c where S – parameter sensitivity vector, R – parameter range vector, T – initial 0 temperature vector and max(.) is maximum element in vector. 2.3 Plant The model of plant used for experimentation was based on the identification carried out for a biomass-fired boiler in [6]. This boiler was a medium-power unit (1 MW) fed by various types of biomass fuels: woodchips, bark, wood waste, agricultural plants and straw. From the control point of view it is a process with very long time delay owing to the fact that it takes approximately 60 seconds to feed the fuel by feeder and some additional time to deliver the fuel by grate into the furnace. Thus, the overall time delay was found out to be around 7 to 8 minutes [6]. The system itself was considered to be of single-input single-output type where the input variable was fuel feed [kg.h-1] and the output variable was water temperature [oC]. It was found out that the system responses to increase and decrease in fuel feed were qualitatively different and therefore two distinct linear transfer functions were identified – one for the case when the fuel feed rate was positive and other one for the case when it was negative (Fig. 4). These two transfer functions (G (s) and G (s)) had the following form: B+ B- 0,0615 0,06 G (s) e480s G (s) e150s (13) B (200s1)4 B (165s1)2 – 29 – A. Hošovský Biomass-fired Boiler Control Using SA Optimized Varela Immune Controller Figure 4 Nonlinear plant transfer function represented by two different linear transfer functions for increasing and decreasing fuel feed (i.e.input variable) (according to [6]) The main disturbances that can be considered for biomass-fired boilers are: thermal losses, fluctuations in heat demand and variable fuel capacity. Of these, the second one is considered also in the presented model. The fluctuations in heat demand will eventually be reflected in the changes of water temperature at the output of boiler. 3 Results and Discussion 3.1 Problem Description and Initialization From the control point of view, the water temperature control system of biomass- fired boilers is specific especially due to its excessive time-delay which significantly obstructs the process of obtaining good performance. As was shown in [6], standardly it was difficult to stabilize the plant using the experience-based tuning of a controller. Thus, some kind of a controller that could not only stabilize the plant but possibly also achieve robust performance for this kind of plant was needed. In Fig. 5 the block diagram of water temperature control using simulated annealing optimized Varela immune controller is shown. The following nomenclature is used : w(t) – desired water temperature, e(t) – water temperature control error, u(t) – fuel feed (control action), v(t) – heat demand fluctuations (disturbance), y(t) – actual water temperature. The control action was upper limited with a value of 350 kg.h-1 (represented as saturation block at the controller output). According to equation 8, Varela immune controller contains several free parameters that affect its performance: – 30 – Acta Polytechnica Hungarica Vol. 12, No. 1, 2015 k – combined parameter of antibody death rate (due to natural death of antibodies f as well as due to their interaction) k – antibody production rate by mature B-cells 3 k - B-cell proliferation rate 4 k – B-cell death rate 5 k – combined parameter of antigen elimination rate and B-cell proliferation rate I due to antigen presence. In addition to these parameters, another three parameters could also be taken into account: initial value of antibody concentration (represented as initial setting of antibody integrator) and scaling coefficients for proliferation and maturation function. These three parameters were not incorporated into optimization and were set in advance based on the trial-and-error results. Furthermore, it was found out that the best results were obtained when k parameter (i.e. B-cell death rate) was 5 set to zero. Even though the number of modifiable parameters was thus reduced to four, it would have been difficult to achieve optimal settings using manual tuning as the interactions between parameters were quite complex. Figure 5 Block diagram of water temperature control using Varela immune controller with simulated annealing optimization To achieve good Varela controller settings (good in the sense of near optimal solution to the optimization problem), simulated annealing algorithm was used. The fitness function for simulated annealing algorithm was in the form of integral absolute error criterion (IAE): n  IAE K T T (14) S rk k k1 where K – scaling constant set to 0.001, T – reference temperature for k-th S rk sample [oC], T – actual temperature for k-th sample [oC], n – number of samples. k – 31 – A. Hošovský Biomass-fired Boiler Control Using SA Optimized Varela Immune Controller Several parameters and functions of the simulated annealing algorithm needed to be set before the optimization. The main three functions (generation and acceptance functions and temperature schedule) forming the basis of the algorithm are expressed in equations (9)-(11). The main control parameters were initial temperature, reannealing interval and number of stall generations (stopping criterion) which were respectively set to the following values: T = 150,  = 100 0 R it and  = 2000 it. SG 3.2 Optimization and Control Test Since the transfer functions representing the system for increasing or decreasing fuel feed were quite different, it was necessary to optimize the controller for both situations separately (it would have been difficult to obtain satisfactory performance with only one set of controller parameter values). Nevertheless, the situations for which different gain values should have been used were considered easily discernible (increasing or decreasing fuel feed) and were not therefore viewed as obstruction to the use of Varela controller. Even when treated separately for increasing or decreasing fuel feed, the whole control system had to be considered nonlinear due to the Varela controller’s nature causing the system performance to be dependent also on excitation signal magnitude. Moreover, in order to avoid extreme immune response of the Varela controller it was necessary to optimize the performance for several consecutive steps of desired water temperature (using only one step would result in very fast yet undesirable response due to the controller’s insensitivity to future changes in reference water temperature). In Fig. 6 the reference temperature step pattern used for optimization is shown. This pattern consists of four consecutive steps in reference temperature change (denoted as T) in the range of 0-8 oC with the time duration of 15000 seconds r for each step. The same step pattern was used for optimization of both increasing and decreasing fuel feed with using transfer functions given in (13). The results of both optimization runs are found in Table 1 with the graphs of best and current fitness values as a function of iteration count shown in Fig. 7. The results obtained using simulated annealing algorithm were improved using a local optimizer (LO). An exponential cooling schedule was preferred over logarithmic or linear due to being much faster than other two schedules. On the other hand, excessive speed of the temperature decrease might be detrimental to optimality of the obtained result. It was thus helpful to use reannealing interval after which the temperature was increased again so that a probability of accepting worse point was higher (and algorithm could jump out of a local minimum). The reannealing points can be seen as spikes in current function value on fitness-iteration graphs in Fig. 7. In both cases the results could be further improved using Nelder-Mead algorithm. – 32 –