APPLICATION OF MULTI-MODEL ADAPTIVE CONThOL TO PULPING PROCESSES by Navid Mohsenian B.A.Sc., Electrical Engineering, University of British Columbia, Vancouver B.C., 1990 A THESIS SUBMITIED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA November 1993 © Navid Mohsenian, 1993 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of &j,CTi(Ci4L The University of British Columbia Vancouver, Canada Date DE-6 (2/88) M1/U/99’ 1’ Abstract In this thesis, the theory of Multi-Model Adaptive Contro l is applied to pulping processes, particularly for wood species compensation and identif ication, and dead-time compensation. Wood species identification is performed on a kraft pulping process whose kinetics are described by the Hatton’ s equations where the ratios of the different wood species in the incoming feed stream are identified on-line by a recursive least-squares based estimation algorithm. Upon wood species identification, a multi-model adaptive controller takes action to compensate for the variability in the wood chips supply content and to control the kappa numbe r. Dead-time compensation is applied to a chlorination tower in a bleaching process with chlorine and chlorine dioxide. The corresponding mass balance equations are developed. The process is simula ted, sampled, and identified in openloop. A multi-model adaptive controller is utilized in closed -loop to compensate for the variations in the pulp quality and transport time delay and to control the bleaching kappa number. The controller’s performance is compared against an adaptive controller which is coupled with a recursive least-squares based algorithm and the Fixed Model Variable Regression Estimation algorithm to compensate for the same variations. U Table of Contents Abstract ii List of Tables v List of Figures vi Acknowledgments xi Introduction 11 I Background 1 1.2 Motivation for this work 2 1.3 Thesis outline 2 2 MMAC as a wood species compensator and identifier 2.1 3 MMAC structure for a kraft pulping process 3 2.1.1 Process kinetics 4 2.1.2 Process kinetics and dynamics 5 2.1.3 The controller 8 2.1.3.a Feedback control 9 2.1.3.b Feedforward control 10 2.1.3.c Feedbacklfeedforward control 12 2.1.4 The weight adapter 13 2.1.4.a Recursive least-squares with exponential forgett ing (RLS) 2.1.4.b Recursive least-squares with exponential forgetting and resetting (EFRA) 2.1.4.c A different mathematical model based on a new perform ance index. 2.1.4.d Identification of more than two species 2.2 MMAC versus conventional adaptive control (CAC) 2.3 Conclusions 15 20 . . . 25 34 39 45 lii 3 MMAC as a dead-time compensator 3.1 46 Dead-time compensation by MMAC 46 3.1.1 The controller 47 3.1.2 The weight adapter 3.2 48 Dead-time compensation by Fixed Model Variable Regression Estimation (FMVRE) 3.2.1 Formulation of FMVRE 3.2.2 Estimation of the process parameters 3.2.3 Implementation 3.3.2 Second-order systems 3.4.1 Kinetics and dynamics of chlorine and chlo rine dioxide prebleaching 3.4.2 Process identification 3.4.3 Noise model 51 53 54 54 64 65 68 73 3.4.4 Closed-loop control issues 74 3.4.5 Simulation results 76 Conclusions 90 General Conclusions and future work 92 4.1 Thesis summary 4.2 50 60 3.4 Adaptive control of a chlorination tower by MMAC and FMVREIRLS (an industrial example) 4 . 54 3.3 MMAC vs FMVRE in adaptive control (simulati on results) 3.3.1 First-order systems 3.5 . 92 Future research 92 References 94 iv List of Tables 2.1 Values of species dependent constants for softwoods in the Hatton’s equation 2.2 Values of species dependent constants for hardwoods in the Hatton’s equation 2.4 Variations of the wood chips supply composition with time 2.3 Values of species dependent constants for western hemlo ck and western red cedar in the Hatton’s equation 5 17 17 2.5 Values of the tuning constants for the exponential forgetting and resetting algorithm 22 2.6 Sinusoidal variations of the wood chips supply compo sition with time 24 2.7 Values of the tuning constants for the exponential forgett ing and resetting algorithm 24 2.8 Values of the tuning constants for the exponential forgett ing and resetting algorithm 27 2.9 Simulation table of a mixture of different wood species . (S/S) softwood/softwood, (H/H) hardwood/hardwood, (S/H) softwood/hardwood 29 2.10 The EFRA constants for the conventional adaptive controller 41 3.11 Time-varying parameters of the simulated process 55 3.12 A comparison between the approximated and the estimated models 73 V List of Figures 1.1 Multi-model adaptive control 1 2.2 Sampling of a first-order continuous-time system 6 2.3 Open-loop block diagram of the single-species proc ess for kraft pulping 7 2.4 Open-loop block diagram of the multi-species proc ess for kraft pulping 8 2.5 Dablin controller in closed-loop 9 2.6 Feedforward control 10 2.7 Combined feedbacklfeedforward control system of the kraft pulping plant 2.8 Multi-model adapter 12 14 2.9 Estimates of w 1 (western hemlock), and its trace of covarianc e for the RLS with exponential forgetting 2.10 Estimates of w2 (western red cedar), and its trace of covariance for the RLS with exponential forgetting 2.11 The control signal (H-factor) 18 19 19 2.12 Estimates of w 1 (western hemlock), and its trace of covarianc e for the EFRA 23 2.13 Estimates of w2 (western red cedar), and its trace of covariance for the EFRA 23 2.14 Estimates of 1 w (western hemlock) with the previous EFRA cons tants 25 2.15 Estimates of w 1 (western hemlock) with the new EFRA constants 25 2.16 Estimates of w 1 (western hemlock) and w 2 (western red cedar) with the EFRA based on the new prediction model 2.17 Process output (kappa number) and the control signal (H-factor) vi 28 28 2.18 Estimates of w 1 (balsam fir) and w 2 (douglas fir) with the EFRA based on the new prediction model 30 2.19 Process output (kappa number) and the control signal (H-fac tor) for a combination of two softwoods 2.20 Estimates of w 1 (trembling aspen) and w2 • . • . • . (beech) with the EFRA based on the new prediction model 2.21 Process output (kappa number) and the control signal (H-fac tor) for a combination of two hardwoods 2.22 Estimates of w 1 (douglas fir) and w 2 (trembling aspen) with the EFRA based on the new prediction model • . • . 2.23 Process output (kappa number) and the control signal (H-fac tor) for a combination of one softwoods and one hardwood 2.24 Estimates of wi at different values of the effective alkali charge when a mixture of H/S is used 2.25 Estimates of w (western hemlock) 2.26 Estimates of W2 (western red cedar) 31 32 32 33 35 35 2.28 Estimates of w 1 (balsam fir) W2 31 35 2.27 Estimates of 3 w (jack pine) 2.29 Estimates of 30 37 (trembling aspen) 38 2.30 Estimates of w 3 (red alder) 38 2.31 Estimates of zv (balsam fir) at [OH}=2 38 2.32 Estimates of w2 (trembling aspen) [OH]=2 39 2.33 Estimates of w 3 (red alder) [OH]=2 39 2.34 Estimates of the process gain (MMAC) 41 vu 2.35 Estimates of the process gain (CAC) 42 2.36 Estimates of the disturbance (MMAC) 42 2.37 Estimates of the disturbance (CAC) 42 2.38 Process output (kappa number) and the control signal (H-factor), (MMAC). 43 2.39 Process output (kappa number) and the control signal (H-factor), (CAC). 44 2.40 Output variance for CAC and MMAC at different weight fractions 44 3.41 The predicted output of individual models 49 3.42 FMVRE-based adaptive control 50 3.43 Process response, control signal, and gain variations (no adaptation) 56 3.44 Process response and control signal (adaptation by MMAC) 56 3.45 Process time delay and gain estimated by MMAC 57 3.46 Process response and control signal (adaptation by FMVREIRLS) 57 3.47 Process time delay and gain estimated by FMVREIRLS 58 3.48 Process response, control signal, and time delay variati ons (no adaptation). 58 3.49 Process response and control signal (adaptation by MMAC ) 59 3.50 Process time delay and gain estimated by MMAC 59 3.51 Process response and control signal (adaptation by FMVR E/RLS) 60 3.52 Process time delay and gain estimated by FMVREIRLS 60 3.53 Open-loop identification of a first-order model parame ters 61 3.54 Control of a second-order system (adaptation by MMAC ) 62 3.55 Process time delay and gain estimated by MMAC 62 3.56 Control of a second-order system (adaptation by FMVR E/RLS) 63 3.57 Process time delay and parameters estimated by FMVR EIRLS. 64 vm 3.58 Prebleaching stage with a mixture of chlorine and chlorin e diode 65 . 3.59 Static input/output relationship of the process 69 3.60 Prebleaching response to a change in chlorine and chlorin e dioxide flow 70 3.61 Bode magnitude plot of the true process and its approx imated 1st-order model 71 3.62 Process identification based on a first-order model approx imation 71 3.63 Estimated parameters of a first-order model before filter installation 72 3.64 Estimated parameters of a first-order model after filter installa tion 72 3.65 Process open-loop block diagram with a noise model 73 3.66 Prebleaching kappa number (process output) and bleach ing agents flow (control signal). Adaptation by FMVREIRLS 3.67 Estimated time delay and gain of the prebleaching proces s by FMVRE/RLS algorithm. . 79 . . 3.68 Prebleaching kappa number (process output) and bleach ing agents flow (control signal). Adaptation by MMAC 80 3.69 Estimated time delay and gain of the prebleaching proces s by MMAC algorithm 80 3.70 Prebleaching kappa number (process output) and bleach ing agents flow (control signal). No adaptation 81 3.71 Prebleaching kappa number (process output) and bleach ing agents flow (control signal). Adaptation by FMVREIRLS 81 3.72 Prebleaching kappa number (process output) and bleach ing agents flow (control signal). Adaptation by MMAC 3.73 Estimated time delay and gain of the prebleaching proces s by FMVREIRLS algorithm. . 3.74 Estimated time delay and gain of the prebleaching proces s by MIvIAC algorithm 3.75 Prebleaching kappa number (process output) and bleach ing agents flow (control signal). No adaptation ix 79 82 . . 82 83 83 3.76 Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by FMVRE/RLS 84 3.77 Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by MMAC 3.78 Estimated time delay and gain of the prebleach ing process by FMVREIRLS algorithm. . 84 . . 85 3.79 Estimated time delay and gain of the prebleach ing process by MMAC algorithm 85 3.80 Prebleaching kappa number (process output) at a closed-loop pole of Q=-0.77 86 3.81 Bleaching agents flow (control signal) at a close d-loop pole of Q=-0.77 86 3.82 Estimated process time delay at a closed-loop pole of Q=-0.77 87 3.83 Estimated process gain at a closed-loop pole of Q=-0.77 87 3.84 Prebleaching kappa number (process output) at a closed-loop pole of Q=-0.4 88 3.85 Bleaching agents flow (control signal) at a close d-loop pole of Q=-0.4 88 3.86 Estimated process time delay at a closed-loo p pole of Q=-0.4 89 3.87 Estimated process gain at a closed-loop pole of Q=-0.4 89 x Acknowledgments I would like to express my thanks and gratitude to Dr. Guy Dumont for giving me the opportunity to pursue my post graduate studies under his guidance and supervision. I wish to thank him for providing me with his invaluable insight and assistance throughout my research. I would also like to thank Dr. Ashraf Elnaggar who was always there to patiently answer my questions and discuss my ideas. My colleagues and good friends Dr. Shiva Shar areh and Mr. Jahan Ghofraniha for their input and support. Mr. Rick Morrison and Mr. Kristinn Kris tinsson for having the computer facilities up and running around the clock. I would also like to than k all my friends and family who have been supportive throughout my study. Finally I am most grate ful to my better half, Susanne, who believed in me and encouraged me, and for her unlimited supp ort and patience during my post graduate studies. xi Chapter 1: Introduction Chapter 1 Introduction 1.1 Background The goal of robust control design is to obtain controllers that maintain desirable performance in the presence of modeffing uncertainties. If modelling unce rtainties can be adequately represented by means of a bank of a finite number of models bounded by the lowest extreme and the highest extreme uncertainties, then the unknown process can be assumed to lie within the models, or more accurately, be characterized by a weighted sum [6] of all those mod els, which is the essence of multi-model adaptive control (MMAC). Figure 1.1 illustrates an over view of such a scheme. The weights reflect n(t) 61 Modell Modem I AdaprT VJ1,W2,...,Wn Figure 1.1: Multi-model adaptive control the probability of a particular model to closely represent the unknown process, and are normalized such that their value stays between zero and one and their sum is one. This is particularly important since the weighted sum of the models should preserve the gain of the unknown process. Weight tuning is achieved by the on-line weight adapter. At every samp le, the weight adapter measures the process output, calculates the outputs of every model, and by computing some modelling error to represent 1 Chapter 1: Introduction a performance index, estimates the weights recursive ly through an estimation algorithm based on minimization of that performance index. The controlle r will then utilize these weights in order to calculate the control action since it will be designed based on the weighted sum of all the models. 1.2 Motivation for this work In pulping processes such as chemical pulping and bleaching, the variety in chemical composi tions of different wood species can often influence control to a degree which would be related to the relative rate of reaction of the species in the presence of chemical agents in order to achieve a desired final product. Most often these variations can caus e process upsets and therefore affect the quality of the final product. Therefore it is desirable to acco unt for the wood species content of chip supply and take compensating actions. This can be achie ved by a multi-model adaptive controller which will not only identify the weight ratios of different wood species in the chip supply by constructing a model for each species, but also compensate for the wood species variations adaptively. 1.3 Thesis outline This thesis is presented in two parts. Chapter 2 form s the first part which includes the formulation of the kraft pulping process kinetics and dynamics , the multi-model adaptive controller design, and the wood species identification algorithm. Finally the performance of the multi-model adaptive controlle r is compared against a conventional adaptive cont roller. Chapter 3 covers the remainder of the thesis. It includes the introduction of two different algorithms for dead-time compensation namely Fixe d Model Variable Regression Estimation and the multi-model approach. The two algorithms are utiliz ed in turn to control generic processes and their performance are compared. Later the equations representing kinetics and dynamics of prebleach ing with chlorine and chlorine dioxide in a chlorinat ion tower are formulated. Again both adaptive algorithms are used to control the chlorination towe r and their performance are compared. Finally chapter 4 gives a summary of this work where concluding remarks are made and future research directions are given. 2 Chapter 2: MMAC as a wood species compensator and identjf ier Chapter 2 MMAC as a wood species compensator and identifier In this chapter MMAC is applied to a kraft pulping process. In kraft pulping, wood chips are exposed to cooking liquor in a vessel where heat is applied. The objective in kraft pulping is to cook wood chips to a target range of kappa number whic h represents the lignin content of the produced pulp. The applied heat will speed the chemical reaction between the cooking liquor and ligni n contained in wood chips. This process is called delignification and the end result is productio n of pulp [28]. The two driving forces for kraft pulp ing reactions are alkali concentration (as measured by either effective alkali or active alkali) and temperatu re. By arbitrarily assigning a relative reaction rate of 1 for 100° C, a method has been developed [29] for expressing the cooking time and temperatu re as one single variable. This variable called the H-fa ctor is applied in cooking control to control the kappa number. Here, by identifying ratios of diffe rent wood species contained in the wood chip s feed stream, the H-factor will be manipulated usin g a MMAC scheme in order to compensate for wood chips variations while the kappa number is kept constant at a desired set point. 2.1 MMAC structure for a kraft pulping process The wood chips supply to a mill consists of one or many species. The variety in chem ical compositions of each individual species can influ ence control to a degree dependent upon the relat ive rates of cooking of the species in order to reach a target kappa number. Therefore it is desirable to know an accurate measure of the ratio of each species pulped in order to achieve a high qual ity final product. In [19] and [16] it has been reco mmended that in order to overcome control prob lems arising from pulping of mixed species, wood chip s from separate storage piles have to be uniformly blended. A method based on multi-model adap tive control has been proposed by Dumont [6], when the raw material consists of two species, that not only will compensate for species variations , but most importantly will provide an on-line estim ation of the wood species content of the raw mate rial. Here, this method is refined, and its applicati on to more than two species is investigated. 3 Chapter 2: MMAC as a wood species compensator and identif ier 2.1.1 Process kinetics In the past 25 years many studies of the kinetics and transport behavior of the kraft pulping process have been performed [15]. With the data from these studies, kraft pulping models of varying complexity have been developed for control purposes . These kinetic models are generally divided into two categories, theoretical and empirical. Theo retical models have been derived based on work by Kerr [20] and Smith [26]. In Kerr’s kinetic mod el, Equation 2.1, the form of the delignification rate equation is a simplistic one: (2.1) L represents the pulp lignin content, k is the rate cons tant, [OH] is the effective alkali concentration, and t is the time. An example of a more complex kinetic model is Equation 2.2 which has been developed by Smith: = (k [ 1 OH] + 2 [k OHJ[S])L (2.2) Comparing Equations 2.1 and 2.2, it can be seen that Smith has embedded [5], the sulfide concentra tion, along with an extra rate constant into his model. Theoretical pulping models have the advantage of being able, at least in principle, to accurately desc ribe the complex behavior of pulping systems and to predict non-uniformity. A disadvantage however , is that it is difficult to obtain sufficiently reliable parameter estimates to justify the effort involved in deriving and applying the models[2l1. Another drawback is that such models are relatively diffi cult to apply and require substantial computing re sources to perform the lengthy numerical integratio n of simultaneous partial differential equations. For these reasons empirical models are more popular. They explicitly relate such dependent variables as yield and kappa number to control variables such as effective alkali and H-factor. Empirical models have their own disadvantages. Because they are entirely without theoretical basis they lack generality. The most popular empirical kinetic mod el is the Hatton’s equation [17], Equation 2.3: K = — /3 (log 10 H) ([OHJ)z . (2.3) It represents a simple relationship that relates kapp a number to H-factor and effective alkali charge under the assumption that sulfidity remains cons tant. Here, K is the kappa number, H is the H-factor, 4 Chapter 2: MMAC as a wood species compensator and iden4fier and [OH] is the effective alkali charge. a, 3, and 0 < n < 1 are species dependent constants that are obtained by curve fitting of experimental data. Hatton has done extensive research to determine the numerical values of these constants for a number of wood species [17, 181. The results have been summarized in Tables 2.1 and 2.2. The Hatton’s equation provides a very compact description SOFTOODS . ,.:... . Wood Species a /3 n Western hemlock 259.3 22.57 0.41 Western red cedar 304.9 33.56 0.35 Jack pine 279.3 30.18 0.35 Alpine fir 340.2 36.81 0.35 Balsam fir 315.0 38.98 0.30 Douglas fir 246.5 23.24 0.40 Table 2.1: Values of species dependent constants for softwoods in the Hatton’s equation. . ....: ... . HAWOODS ..: . .. ..,. ,.. Wood Species a /3 ii Trembling aspen 124.7 5.03 0.76 Beech 106.3 3.55 0.80 Hard maple 86.9 4.19 0.65 Red alder 87.6 4.83 0.60 Yellow birch 75.2 2.65 0.75 .,., Table 2.2: Values of species dependent constants for hardwoods in the Hatton’ s equation. of the pulping behavior of several species, and therefore it will be utilized here to represent a kinetic model for kraft pulping when the wood chip supply consists of more than one species. 2.1.2 Process kinetics and dynamics Assuming kraft pulping takes place in a hypothetical reactor that can be approximated by a first-order dynamics plus a time delay, the sampled system can be described (as shown in Figure 2.2) by the open-loop transfer function, HOL (s): e_Tds HOL(S) +1 5 (2.4) Chapter 2: MMAC as a wood species compensator and identifier I 1_:3T8 Figure 2.2: Sampling of a first-order continuous-time system preceded by a zero-order hold HZOH(s), where: 1 HZOH(s) — (2.5) = Here, g is the process gain, Td the time delay, T the time constant, and T 8 is the sampling time. The equivalent pulse transfer function of the system in Figure 2.2, in discrete-time domain, can easily be calculated as[27}: Tr1()Jq —3T ( 1 j—e )= s ge —Ts rs+l (2.6) g(1 + a)q’ —k q — 1+aq’ — where g, as before, is the process gain, q’ the backward shift operator, (i.e. 1 q y (t) = y(t — 1)), and a and k are the process pole and time delay respectively defined as: a = _e_TsIT , k (2.7) Now considering the kraft pulping of n different species, as discussed in section 2.1.1, the kinetics of each species can be described by the Hatton’s equation, relating the final kappa number K, to the H-factor H, and effective alkali [OH], where i denotes the species index: = — . 10 H) ([OH])’’ (log for i = 1,2,..., n . (2.8) Then by relating the input and output of the process, Equation 2.8, to the process dynamics, Equation 2.6, and assuming a unit gain, the following equality can be reached: — — . (log 10 H) ([OH}) . 6 — (1+a)q —k q 1 1+ aq 29 (.) Chapter 2: MMAC as a wood species compensator and identifier The equivalent difference equation for every species can then be obtained by cross-multiplying Equation 2.9 to obtain: K(t) —aKi(t — 1) + (1 + a)aj(t — 1 — k) — (1 + a){f3( [OH])’’ (log 10 H)}t_l_k (2.10) for i=1,2,...,n Equation 2.10 can then be written in the standard form of: K(t) = —aK(t — 1) + bu(t — 1 — k) + cv(t — 1 — k) (2.11) for i=1,2,...,n or ‘ K(t) = q_ku(t) 1 + q_kvi(t) 1 aq— (2.12) where u(t — 1 — k) = 10 H)t_l_k (log v(t—1—k) =a(t—1—k) (2.13) = —(1 + ( 3 a){ [OH])}_l_k / c=1+a Equation 2.12 represents a single-species model for the kraft pulping process. An open-loop block diagram of this model is depicted in Figure 2.3. vj(t) u(t) K(t) Figure 2.3: Open-loop block diagram of the single-species process for kraft pulping. Equation 2.11 can be expanded to account for a mixture of species by taking the weighted sum of both sides of the equation. Let 0 w < 1 be the weight fraction of species i into the reactor, and K be the final kappa number when more than one species are present. Then the final kappa 7 Chapter 2: MMAC as a wood species compensator and identifier number out of the reactor can be represented as the weighted sum of the kappa number contributions of each individual species: K(t) = wK(t) (2.14) By substituting Equation 2.12 into Equation 2.14 we obtain: K(t) w{—aK(t — 1) + bu(t — 1 — k) + cvj(t — 1 — k)} (2.15) = Expansion of Equation 2.15 leads to: K(t) = —aK(t — 1) + wbu(t — 1 — k) + cwv(t — 1 — k) (2.16) or E wbq’ K(t) = 1 aq —1 q_ku(t) + q_k wv(t) 1 +aq_’ (2.17) which represents the multi-species model for the kraft pulping process. An open-loop block diagram of the model of Equation 2.17 is depicted in Figure 2.4. u(t) K(t) Figure 2.4: Open-loop block diagram of the multi-species process for kraft pulping. 2.1.3 The controller In section 2.1.2, a first-order open-loop model for delignification of pulp in kraft pulping was derived. In this section we turn our attention to closing the loop and designing a controller for this model. Looking at Figures 2.3 and 2.4, it is evident that the control signal u(t) and the “load disturbance” v(t), both affect the final kappa number K(t). Therefore, in order to keep the kappa number at a desired set point and compensate for the load disturbance, a feedforward controller along with a feedback controller have to be designed to control the process. 8 Chapter 2: MMAC as a wood species compensator and identifier 2.L3.a Feedback control One of the most widely used feedback controllers in pulp and paper industry, for processes with a significant time delay [5], is the Dahlin controller. The Dahlin controller was first introduced [4] in 1968, and since then it has proven to be very useful for dead-time compensation provided an accurate model of the process under control exists [5]. ) 1 HCL(q Figure 2.5: Dahlin controller in closed-loop A Dahuin controller is designed such that the closed-loop response required is that of a first-order system with a gain of one and with a time delay equal to the time delay of the open-loop system. Suppose the open-loop system is described by: HOL = B( _1) q_k A(q’) (2.18) where q 1 is the backward shift operator, B and A are polynomials in terms of q , and k is the 1 open-loop time delay. Let the desired first-order closed-loop response of gain one and time delay of k be described by: = HCL where Q (1+ 1+Qq (2.19) denotes the closed-loop pole. Then according to the closed-loop configuration depicted in Figure 2.5 the following relationship holds: HFB HCL TT . HOL 1 + FB TT LOL (2.20) where HFB is the feedback Dahlin controller. Equation 2.20 can be easily solved for HFB to yield: = HFB 1 HOL 9 HCL l—HCL (2.21) Chapter 2: MMAC a. a wood species compensator and identifier Substitution of Equations 2.18 and 2.19 into Equation 2.21 yields Equation 2.22, which is the Dahlin controller: 1 A B(q’) FB 1 1+ Qq’ Consequently, by substituting for A(q’) —1 — (1+ Q)qk_l 1 + aq 1 and B(q ) 1 = bq’ from Equation 2.12 into Equation 2.22, the single-species kappa number feedback controller based on a Dalilin controller is obtained: HFBjq —1\j=. . 1 ,, u(t) ref 1+aq’ (223 — 1+Q 1 + Qq’ (1 + — Q)q_k_l — where Kr€f (t) is the reference signal. Similarly the multi-species kappa number feedback Dablin controller can be obtained by substituting for A(q’) = 1 + aq 1 and B(q’) = D wbq’ from Equation 2.17 into Equation 2.22: HFBq , ‘ = — n(t) Kre(t) K(t) 1 + aq 1 1+Q 1 (1fQ)q_k_I 1+Qq Ewjbj — (2.24) — — 2.1.3.b Feedforward control If disturbances can be detected and measured as they enter a system, suitable compensating control actions can be introduced in order to keep the output variables constant. Such adjustment in control actions is called feedforward control [22], and is depicted in Figure 2.6. The load disturbance v(t) enters the process through the Hv(q ) process transfer function. The load disturbance is also 1 fed into a feedforward controller that has a transfer function HFF (q’). The feedforward controller v(t) —‘ H,,(q’) HFF (q’) + I I u(t) ‘] + ) 1 Hj(q I Process Figure 2.6: Feedforward 10 control y (t) Chapter 2: MMAC as a wood species compensator and identjfier detects changes in the load disturbance v(t), and makes compensating changes in the control action u(t) so that contributions to the process output due to disturbances are zero [221: 0 = v(t) Hu(q’) + v(t) HFF(q’) Hoi,(q’) (2.25) Solving for HFF(q’), the feedforward controller relating the control action, u(t) to the load disturbance, v(t) is: 7:11-1 HFF(q —1 ) = — ) HOL(q —1 (2.26) However, the designer has to bear in mind that Equation 2.26 is only feasible if two conditions are satisfied. The first condition regards the stability of the feedforward controller where the open-loop transfer function HOL(S) has to have a stable inverse (i.e. be minimum phase). The second condition regards maintaining the causality of the feedforward controller where the time delay of HOL(S) should be less than that of H (q’). Having assumed these two conditions, and substituting for clq_k Hv(q’) (2.27) and ‘jq HOL (q’) q_k = (2.28) from Equation 2.12 into Equation 2.26, the single-species kappa number feedforward controller is obtained: , —1\ HFFq j u(t) = vzi) c = —- (2.29) 2 Similarly the multi-species kappa number feedforward controller can be obtained by substituting for Hv(q’) q_k 1 1 + aq (2.30) and HOL (q ) 1 E wbq q_k aq (2.31) from Equation 2.17 into Equation 2.26 to obtain: HFF(q) C = — wb i=1 11 (2.32) Chapter 2: MMAC as a wood species compensator and identifier 2.1.3.c Feedbacklfeedforward control In sections 2. l.3.a and 2. l.3.b the transfer function of a feedback Dablin controller and a feedforward controller for the case of a single-species and multi-species wood chips in a kraft pulping process were derived. Figure 2.7 shows the complete control systems for the process under study where HFB (q’) and HFF(q’) are Equations 2.23 and 2.29 for the case of a single-species, and Equations 2.24 and 2.32 for the case of multi-species respectively. The control signal is influenced Figure 2.7: Combined feedbacklfeedforward control system of the kraft pulping plant. by both feedback and feedforward controllers. The process and its control algorithm can then be completely described by: bq’ K(t) = u(t) 1 + aq 1 1+aq’ = 1 + Qq’ q_ku(t) + 1+Q (1 + — 1 cq q_kvi(t) (2.33) 1 1 + aq Q)q_k_I (Kref(t) — K(t)) c — —vi(t) IN (2.34) for the single-species case and n wbq’ K(t) — — u(t)= 1 1 + aq 1 + aq 1 wb 1 + Qq’ 1+Q (1 + — q_ku(t) + Q)q_C_l 1 qk cq 1 + aq 1 (Kref(t) — wv(t) C K(t)) — E i=1 (2.35) i=1 1 wb > wv(t) (2.36) i=1 for the multi-species case respectively. In [6] it has been shown that the multi-species controller, Equation 2.36, can be expressed as a linear combination of the single-species controller, Equation 12 Chapter 2: MMAC as a wood species compensator and identifier 2.34, i.e.: u(t) (2.37) = where wb = (2.38) wb i=1 The multi-species controller, Equation 2.36, assumes knowledge of the weight ratios of the wood supply composition at every sample. If these ratios are unknown, then in order to calculate the correct adaptive control action, the correct values of either w or \ have to be recursively identified. Here, we choose direct identification of the w which is the subject of the proceeding section. 2.1.4 The weight adapter In section 2.1.2 a model representing kinetics and dynamics of kraft pulping when raw material consisted of more than one species was developed. Section 2.1.3 was dedicated to the design of a multi-species controller associated with the developed kraft pulping model based on a Dahlin controller and feedforward control. The focus of this section is the estimation of the weights associated with every species. Figure 2.8 shows a block diagram of the proposed multi-model adapter. The raw material containing wood chips of different species are fed into the digester where they are cooked to a final kappa number K(t). A bank of n models represents kinetics and dynamics of the n species present in the raw material. Each model, expressed by the Hatton’s equation and a first-order dynamics delayed by the process time delay, represents the kraft cooking of each individual species. The models are driven by the control signal, Equation 2.36, and the output of every model, K(t), is calculated according to Equation 2.11. The objective is to minimize a performance index that will lead to the recursive identification of the weights. One of the most popular and basic methods for parameter estimation is least-squares and its variants, based on the principle formulated by Gauss. According to this principle if a mathematical model representing a process has the property of being linear in parameters, then the unknown model 13 Chapter 2: MMAC as a wood species compensator and identifier n(t) u(t) Figure 2.8: Multi-model adapter parameters should be chosen such that the sum of the squares of the difference between the actual measured output of the process and its predicted value computed by the mathematical model is a minimum [2]. Suppose the predicted output is represented by the regression model: (t) = (t) 1 8 (t) + 2 1 + 02q + Omçbm(t) (2.39) = with (t) [ ( 1 t) ( 2 t) = (2.40) 8= [8i 82 where (t) is the predicted output, functions, and 01,62,. , i, q,. . . , •.. Om]T qS, are called regression variables and are known m are the unknown model parameters to be estimated. If the measured 8 process output is denoted by y(t), the least-squares error between the measured output and the predicted output computed by Equation 2.39 leads to the following performance index V(8) = ‘N {y(t) t = 1 - (2.41) > — 14 Chapter 2: MMAC as a wood species compensator and identifier where N is the number of measurements and 0 < A 1 called the forgetting factor, is a scalar that gives us the flexibility to base the identification on the most recent data rather than on the old one in case the parameters are time-varying. In the following sections a mathematical model based on Equation 2.39, relating the predicted final kappa number to the weights associated with every species model will be constructed. This model will then be utilized in Equation 2.41, and the perfonnance index will be minimized with respect to the model parameters to obtain a recursive least-squares based algorithm for the identification of the weights. 2.1.4.a Recursive least-squares with exponential forgetting (RLS) A mathematical model in the form of Equation 2.39 can be directly constructed from Equation 2.14 to represent the predicted kappa number as: I(t) = zii(t —1). 1< (t) + ti 1 (t —1). K 2 (t) + 2 +thn(t 1). K(t) (2.42) = ( T t)8(t— 1) with q(t) = [K ( 1 t) 2 K ( t) • . T K(t)J . (2.43) 0(t) = [ii(t) 2 z1 ( t) T th(t)] •.. where (t) is the predicted final kappa number, 1 K ( t), 2 K ( t), ..., K,(t) are the regression variables calculated from the family of n models, and th 1 (t), r 2 (t), ..., zI (t) are the model parameters representing the estimates of the unknown weights associated with every model, normalized due to the constraint Ew 1 as: zi(t) w(t)= for j = 1,...,n (2.44) Ezij(t) i=1 Then the minimization of Equation 2.41 which is now written as: V(0) 1 = A_t{K(t) - 2 f(t)} (2.45) with respect to the model parameters will lead to the following recursive least-squares algorithm [14]: e(t) = K(t) — 4 T (t)Ô(t 15 — 1) (2.46) Chapter 2: 0(t) MMAC as a wood species compensator and identifier P(t—1)q5(t) + A + ,T(t)p(t l)(t)0(t) P(t 1)(t)T(t)P(t 1) (t- )(t)P(t-1)(t) T A+ 8(t — 1) (2.47a) — (t)_ — — • 1 2 47b f where e(t) is the computed prediction error between the measured and the predicted kappa number which is fed into the weight tuning block, and Ô(t) is an estimate of 0(t), calculated recursively based on its old estimate, and inside the weight tuning block as shown in figure 2.8. P(t) is the covariance matrix whose value is proportional to the variance of the estimate ö(t). The forgetting factor A prevents the covariance matrix from decaying, and therefore maintains the algorithm’s adaptivity for time-varying parameters. In summary estimation of the weights is implemented in the following steps: 1. Choose an initial guess for the parameters and the covariance matrix: ö(O) 2. Calculate the kappa number K (t) for every model, (Equation 2.11). 3. Construct the parameter vector 4. Calculate the prediction error e(t), (Equation 2.46). 5. Estimate a new 9(t) and update the covariance matrix P(t), (Equation 2.47). 6. Normalize the estimated weights, (Equation 2.44). 7. New sample; go to step 2. Ô(t) = , P(O) = 0 P and the regression vector q(t), (Equation 2.43). Simulation results For simplicity, kraft pulping of a mixture of only two softwoods namely western hemlock and western red cedar is first simulated. The corresponding constants for the Hatton’s equations describing kraft pulping of theses softwoods are shown again in Table 2.3 for convenience. The last column shows the range of kappa number for each species where the constants a, 3, and n fit the experimental data [17]. Therefore the target kappa number Kref has to lie within these ranges. According to [3] the SOFTWOODS Wood Species a Range of K Table 2.3: Values of species dependent constants for western hemlock and western red cedar in the Hatton’ s equation. (Continued) 16 Chapter 2: ¶::::.:. .:::... MMAC as a wood species compensator and identifier ::. Thner 0t<500 wj(t) w ( 2 t) 0.1 0.9 500 t < 1000 0.2 0.8 1000 t < 1500 0.3 0.7 Table 2.4: Variations of the wood chips supply composition with time. [ [ Western hemlock 259.3 22.57 0.41 20.6-49.0 Western red cedar 304.9 33.56 0.35 20.0-51.0 Table 2.3: Values of species dependent constants for western hemlock and western red cedar in the Hatton’s equation. operating range for kraft pulping of softwoods is at 28—35 and for hardwoods at 13—15. Here, Kref is chosen to be 3 Krf= O . The corresponding effective alkali charge for the chosen kappa number target can be found in [17] to be [OH]=18. Since kraft pulping process like most pulping processes has a dominant time delay relative to its time constant, a slow sampling rate is usually chosen to account for the long time delay, and therefore to achieve a sensible number of delay samples in the discrete-time domain. A slow sampling rate will in turn push the process pole closer to the centre of the unit disk according to Equation 2.7. For this reason the open-loop process time delay and pole are chosen as k=5 and a=—0.6 respectively. Assuming a faster closed-loop response is desired, the closed-loop pole is set to Q=—0.5. The closedloop time delay is unchanged at k=5. Furthermore, an independent white noise sequence of zero-mean and variance 2 a = 0.25 was superimposed on the final kappa number in Equation 2.17, to make the simulation more realistic: K(t) (biw Q 1 )+ w2(t))q 5 + 1 0’ 1 (w,(t)v,(t) + 2 5 q 1 (t)v w ( t)) 0. q q + n(t) 1 0.6q’ (2.48) — where b,, 2 b 1 , v and v , 2 are calculated from Equation 2.13 for each species, n(t) is the white noise sequence, and the wood chips supply composition is varied with time according to Table 2.4. where the weights are normalized. Then assuming the controller has prior knowledge of the process pole 17 Chapter 2: MMAC as a wood species compensator and identifier and time delay, the control signal is calculated according to Equations 2.36 as: u(t) 1—06—’ 05 b,zi(t 1) + b (t —1) 1— 0.5q’ z 2 0.4(th,(t 1)v,(t) + th (t 1)v2(t)) 2 1) + 2 I) b ( t 1) — — — — 6 0.5q ref(t) — K(t)) (2.49) — — — 2 are the estimates of wi and where th 1 and ti) W2 respectively. Figures 2.9 and 2.10 show the true and the estimated weights and the corresponding trajectories of the covariance matrix P(t), when the recursive least-squares algorithm with a forgetting factor of X=O.98 is used. Although for the wi and its estimate 1 0.8 0.6 0.4 0.2 200 400 200 400 600 800 Time in samples Trace of the covariance matrix 2 1 1’ “0 600 800 Time in samples 1000 1200 1400 Figure 2.9: Estimates of w 1 (western hemlock), and its trace of covanance for the RLS with exponential forgetting. case of only two wood species it is sufficient to estimate only one weight and calculate the other as 1 — ti(t), here we have decided to estimate both weights independently. The dotted and the solid lines represent the true and the estimated weights respectively. Evidently the estimator does not have a satisfactory performance as the weights have large variances. This can be explained by analyzing the trace of the covariance matrix for both estimates. Under low excitation, since q(t)q5T(t) 0, the bracketed term in the covariance matrix update, Equation 2.47b, decays slower than it is inflated by the multiplier 1/A. This leads to the exponential growth of P(t) as P(t) 18 P(t — 1)/A and Chapter 2: w2 MMAC as a wood species compensator and identifier and its estimate Trace of the covariance matrix 600 800 Time in samples Figure 2.10: Estimates of w 2 (western red cedar), and its trace of covariance for the RLS with exponential forgetting. therefore large variations occur in the estimated parameters as shown in Figures 2.9 and 2.10. It is H-factor 600 800 Time in samples Figure 2.11: The control signal (11-factor). interesting to note that the rapid decay of the covariance matrix at t=500 and 1000 samples is due to changes in the variation of the wood chips supply composition. The control signal, shown in Figure 2.11, will adjust itself to compensate for these variations. This compensation will lead to a change in the calculated kappa number out of each model namely the regression variables. Consequently rich excitation will be introduced into the regression vector q(t). Now the bracketed term in Equation 2.47b will decay faster than it is grown by the forgetting factor, which will lead to the decay of 19 Chapter 2: MMAC as a wood species compensator and identifier the trace of the covariance matrix. Subsequently the estimates produced by the estimator are closer to their true values. As it was seen the recursive least-squares with exponential forgetting algorithm prevents the covariance matrix from going to zero and therefore keeps the estimator alert. However, as it was observed, one difficulty with this algorithm is that with poor excitation the covariance matrix can grow without bound which contributes to undesirable variance of the estimates. One alternative is to allow the estimator to smoothly reset F(t) to some sensible value whenever the excitation is poor. This algorithm called the exponential forgetting and resetting algorithm (EFRA) [23] is another variant of the recursive least-squares, and will be utilized here to improve the performance of the weight estimator. 2.1.4.b Recursive least-squares with exponential forgetting and resetting (EFRA) As the name of the EFRA suggests, not only does the algorithm incorporate exponential forgetting of older data, but also features upper and lower bounding of the covariance matrix under low excitation conditions, and resetting it to its upper bound. Choosing the same mathematical model representation of the process as in Equations 2.42 and 2.43, the corresponding exponential forgetting and resetting algorithm is [23]: e(t) 0(t) P(t) where i, = P(t = K(t) — (cb T t)Ô(t — 1) l)4(t) l)()e(t) + 4T(t)p(t + 1 iiP(t—i)c5(t) ( T t)P(t— b 1) (t 2 + p1— 6P 0(t — — = — (2.50) — 1.) (2.51a) — — 1) (2.51b) p, and 6 are constants that tune the estimator, and I is the identity matrix. If the estimator is to live up to its expectations and perform well whenever the excitation is poor, the four tuning constants have to be chosen under some constraints. It is shown in [23] that if the four constants satisfy the following constraints: (2.52a) 20 Chapter 2: ( MMAC a.s’ a wood species compensator and identifier + 4iz6 < (1 — p>0, oI — 6>0 P(0) (2.52b) (2.52c) VI (2.52d) where 1A ( (2.53a) + (2.53b) + (‘) (_ + + (q7)2) (2.53c) and I is the identity matrix, then for an arbitrary regression vector c(t) and all t, the covariance matrix P(t) is bounded from above by VI and from below by ãI, and under low excitation conditions, i.e., (t)=O, the covariance matrix will approach 171 asymptotically which is the exponential resettling property of the EFRA. Although the algorithm contains four timing constants, the selection of these constants in practice is relatively easy and is done according to the following general guidelines [231: 1. q adjusts the estimator’s gain; typically q E [0.1, 2. i is a small constant which is directly related to the minimum eigenvalue of P; typically ft E [0, 0.01]. 0.5]. 3. A is the forgetting factor; typically A E [0.9, 0.99]. 4. 6 is a small constant which is inversely related to the maximum eigenvalue of F; typically 5 E [0, 0.01]. With the above choices of parameter ranges a and P are approximated as [23]: 67 a 17—7 21 (2.54a) (2.54b) Chapter 2: MMAC as a wood species compensator and identjfier Simulation results For this simulation the process and the controller’s parameters remain unchanged as before. The only difference is the implementation of the new estimator, the EFRA, in the place of the recursive least-squares estimator. The wood chips supply composition to the plant remains unchanged as shown in Table 2.4 for the purpose of performance comparison between the two estimation algorithms. The estimator’s four constants are chosen according to the above guidelines and are displayed in Table 2.5. With these choices of the estimator’s constants, the constraints in Equations 2.52a and 2.52c EFRA CONSTANTS 0.5 0.005 A 6 0.98 0.005 Table 2.5 Values of the tuning constants for the exponential forgetting and resetting algorithm. are obviously satisfied, and for 2.52b we have: (,, — + 4i6 = 0.23 < 0.25 = (1 — 77)2 (2.55) To find an initial value P(0), for the covariance matrix we have to satisfy Equation 2.52d. Numerical values of approximate lower and the upper bounds of P, namely ö and V are calculated from Equation 2.54 to be 0.01 and 4.3 respectively. For this we choose P(0)=41. Figures 2.12 and 2.13 show the true and the estimated weights and their corresponding trajectories of the covariance matrix, P(t), when the exponential forgetting and resetting algorithm is used. Again the dotted and the solid line represent the true and the estimated weights. It is obvious that the estimated weights follow their true values closely. The trace of the covariance matrix for both estimates does not decay to zero; therefore the algorithm retains its alertness and adaptivity in order to estimate the time-varying weights. Nor does the P matrix grow due to its boundness unlike the recursive least-squares algorithm. As a 22 Chapter 2: MMAC as a wood species compensator and identifier wi and its estimate (western hemlock) I I 0.8 0.6 200 400 600 800 1000 1200 1400 I I I 1000 1200 1400 Time in samples Trace of the covariance matrix 3 I I 2 1 200 400 600 800 Time in samples Figure 2.12: Estimates of w 1 (western hemlock), and its trace of covariance for the EFRA. w2 and its estimate (western xd cedar) 1 I 0.8 0.6 0.4 0.2 “0 200 400 600 800 Time in samples 1000 1200 1400 1200 1400 Trace of the covariance matrix 3 I 2 200 400 600 800 Time in samples Figure 2.13: Estimates of w 2 (western red cedar), and 23 1000 its trace of covariance for the EFRA. Chapter 2: MMAC as a wood species compensator and identifier result the estimated weights stay close to their true values. However, further simulations revealed that although the EFRA provides superior performance over the RLS algorithm in terms of weight estimation, the same choice of estimator’s constants, shown in Table 2.5, does not generalize the problem of weight estimation in the case of sinusoidal weight variations, and the estimator’s constants have to be slightly modified in order to account for these types of variations. Assume now that the composition of the wood chips supply varies with time according to Table 2.6. In order to make 0t<250 2 8 250 t<1000 2+cos(0.Olt+O.5) 2 1000t<1500 8 2 Table 2.6: Sinusoidal variations of the wood chips supply composition with time. the estimator more active, the estimator’s gain is increased to one and the forgetting factor .X is decreased to 0.95. The new EFRA constants are displayed in Table 2.7. The lower and the upper bounds of the covariance matrix are again calculated from Equation 2.54 to be a=0.005 and ii=10.6 respectively. For this we choose P(0)=10I. EFRA CONSTANTS ‘1 II 1 0.005 (S 0.95 0.005 Table 2.7 Values of the tuning constants for the exponential forgetting and resetting algorithm. Figures 2.14 and 2.15 show the difference in the estimation of the weights when different sets of the EFRA constants are used. For the purpose of comparison only the estimates of w 1 are shown here. Clearly the new EFRA constants keep the estimator alive where the old EFRA constants fail. However, there is a price to be paid here and that is an increase in the variance of the estimated weights. This was to be expected according to Equations 2.53a and 2.54a since a 0.3 decrease in ). from 0.98 to 0.95 means more than a double increase in y from 0.02 to 0.05. Therefore the upper bound of P which is a function of ‘y is pushed up twice as much as before. This leads to a looser upper bound and, consequently an increase in the variance of the estimates. The search for finding 24 Chapter 2: MMAC as a wood species compensator and identifier a solution to compensate for the large variance of the estimates led to a different approach to the problem which is the topic of the next section. wi and its estimate (western hemlock) 600 800 Time in samples Figure 2.14: Estimates of w 1 (western hemlock) with the previous EFRA constants. wi and its estimate (western hemlock) Figure 2.15: Estimates of w 1 (western hemlock) with the new EFRA constants. 2.1.4.c A different mathematical model based on a new performance index In sections 2.l.4.a and 2.l.4.b the performance index to be minimized with respect to the model parameters was chosen as the sum of the squares of the prediction error between the measured final kappa number out of the process and the predicted one calculated by the mathematical model of Equation 2.42. Suppose that now we set our goal to minimize the sum of the squares of the prediction error between the calculated kappa number for each individual species from each individual model, and the same kappa number ‘extracted’ from the measured final kappa number out of the process. Then in this case the least-squares error for only two species leads to the following performance index: V(8) = > ({it - 2 ( 1 t)} + {K (t) 2 25 - 2 ( t)}2) (2.56) Chapter 2: MMAC as a wood species compensator and identifier 1 (t) and K where K 2 (t) are kappa numbers calculated from the model of each individual species as: K(t) —aK ( 1 t 1) + bu(t — — k) + cvj(t — 1 — k) (2.57) for and 1 — i=1,2 1 (t) and R R 2 (t) are the same kappa numbers ‘extracted’ from the measured final kappa number in the following fashion: K(t) i(t) 2 t) K ti ( 2 t — wi(t K(t) — — — — — — w ( 2 t 1)K ( 2 t) 1) (t) 1 1)K 1) — — (2.58a) (2.58b) where K(t) is the measured final kappa number. Expansion of Equation 2.58 leads to the following mathematical model representation of K 1 (t) and K 2 (t): 1 — 1) = KQ) th ( 2 t — ‘K ( 2 t) 1) — — (2.59) 1) — with ç’i(t) o and t ( 2 t) = = I = — 1 (t 2 w — [K(t) FI — zi, ( 2 t)]T 1 c’(t)O2(t K(t) — (2.60) wi(t) LW1(t) 1) T ( 2 K t)] wi(t — — t ( 2 t — (t) K 1 1) (2.61) 1) with 2 ( t) o [K(t) ‘ 2 I — FI — (2.62) T 1 thi(t) 1 LW2(t) w ( 2 t) It should be mentioned that here, identification of the parameters is performed by estimating the reciprocal of the desired weights, therefore in order to avoid possible division by zero a small positive constant can be added to the weight estimates. Then the minimization of Equation 2.56 with respect to the model parameters O and 02, unlike the previous performance index, will lead instead of one, to two recursive estimators based on the exponential forgetting and resetting algorithm as shown here: eQ) = K(t) — 1T(t)0(t 26 — 1) (2.63) Chapter 2: O(t) P ( 2 t) = — = 1) 1) (2.64a) + 1+T(t)l(t_l)I(t)et(t) 1)4ij(t)4(t)P(t 1) jPj(t + Pz’ — — MMAC as a wood species compensator and identifier — — 1) (2.64b) where index i=1,2 represents the species number. Clearly the difference between this approach to weight identification and the one discussed in sections 2.1 .4.a and 2.1 .4.b is the choice of the minimization of a different prediction error which leads to construction of a new mathematical model. Evidently such a selection allows each species to be identified separately by an estimator of its own which differs from the previous case where the identification was achieved by only one estimator. Furthermore, the weight estimates for each species are calculated by taking the reciprocal of the first element where the second element of the parameter vector 8(t) is disregarded in parameter estimation. Simulation results Again for this simulation, the process and the controller’s parameters stay unchanged. The only difference is the implementation of the new prediction model. The wood chips supply composition is again varied sinusoidally according to table 2.6 and the EFRA constants for both estimators are shown in Table 2.8. For this choice of the EFRA constants, the lower and the upper bounds of the . FRA CONSTA T TS (S 1i 1 i’. 1 1 0.01 0.95 0.001 2 1 0.01 0.95 0.001 Table 2.8 Values of the tuning constants for the exponential forgetting and resetting algorithm. covariance matrix for both estimators can again be calculated from Equation 2.54. Figure 2.16 shows the estimated weights when the new performance index with respect to the new model parameters is minimized. Figure 2.17 shows the corresponding process output (kappa number), and the control action (H-factor). As can be seen not only does the multi-model adaptive controller provide good control of the kappa number despite the wide variation of the wood chips supply composition, but it also gives an accurate measure of the composition of the incoming wood chips. 27 Chapter 2: MMAC as a wood species compensator and iden4fier wi and its estimate (western hemlock) 0 G 200 400 600 800 1000 1200 1400 1200 1400 Time in samples w2 and its estimate (western red cedar) 200 400 600 800 Time in samples 1000 Figure 2.16: Estimates of w 1 (western hemlock) and w 2 (western red cedar) with the EFRA based on the new prediction model. Final kappa number (softwoods) 2 I 200 30C 400 I 600 800 Time in samples H-factor I 1000 1200 1400 I I I 1200 Ii 80% 200 400 600 800 Time in samples 1000 1200 1 1400 Figure 2.17: Process output (kappa number) and the control signal (H-factor). 28 Chapter 2: MMAC as a wood species compensator and iden4fier The new scheme shows a great improvement with respect to the variance of the estimates over the previous one showed in Figure 2.15. It is believed that the reduction in the variance of the estimates is due to the method by which the estimates are calculated; namely by the taking the reciprocal of the elements of the parameter vector. Furthermore, although the estimator provides fairly accurate estimates of the weights, the estimates do not follow their true values very tightly. The reason for that is the designer’s choice to select the estimator’s tuning constants ‘q, )., .t, and 6 such that the estimator is robust enough to accommodate different species of different characteristics. Figures 2.18 through 2.23 show a series of simulation of kraft pulping for such species of different natures. Throughout all these simulations the EFRA constants remained unchanged at their values shown in Table 2.8. However, the kappa number set point Kref and the effective alkali [OH] were varied according to the type of the wood species used in the simulation. Table 2.9 shows the wood species type, their corresponding kappa number set point and effective alkali along with their corresponding figures. Comparing the estimated weights for different wood mixtures in Figures 2.18, 2.20, and 2.22 Mixture Wood Specie.s Krej [OH] Figure S/S balsam fir / douglas fir 30 18 2.18-2.19 HIH trembling aspen I beech 15 14 2.20-2.21 S/H douglas fir / trembling aspen 25 15 2.22-2.23 Table 2.9 Simulation table of a mixture of different wood species. (S/S) softwood/softwood, (HJH) hardwood/hardwood, (S/H) softwood/hardwood. it is observed that the weights have the tendency to follow their true values closer when the wood mixture is composed of a hardwood/softwood (HIS) than when the composition is made of S/S or HJH. The reason for that can be tracked back to the structure of the Hatton’s equations describing the kinetics of the wood species: 1 K 2 K = — — . /2 10 H) ([OH])1 (log . 10 H) . ([OH])n2 (log (2.65) As can be seen the Hatton’s equations are distinguished by the two constants a and /3g. Although the equations are not an exact multiple of each other, a value for the effective alkali, [OH], can be 29 Chapter 2: MMAC as a wood species compensator and identifier wi and its estimate (balsam fir) 600 800 Time in samples wl and its estimate (douglas fir) I I I 200 400 -- 0.8 0.6 0.4 0.2 600 800 Time in samples 1000 1200 1400 Figure 2.18: Estimates of w 1 (balsam fir) and w 2 (douglas fir) with the EFRA based on the new prediction model. Final kappa number (softwoods) H-factor Figure 2.19: Process output (kappa number) and the control signal (H-factor) for a combination of two softwoods. 30 Chapter 2: MMAC as a wood species compensator and identifier wi and its estimate (trembling aspen) 200 400 200 400 1000 600 800 Time in samples w2 and its estimate (beech) 600 800 1000 1200 1400 1200 1400 Time in samples Figure 2.20: Estimates of w 2 1 (trembling aspen) and w (beech) with the EFRA based on the new prediction model. Final kappa number (hardwoods) ‘)1’ 200 400 600 800 1000 1200 1400 1000 1200 1400 Time in samples H-factor 0 200 400 600 800 Time in samples Figure 2.21: Process output (kappa number) and the control signal (H-factor) for a combination of two hardwoods. 31 Chapter 2: MMAC as a wood species compensator and identifier wi and its estimate (douglas fir) 1 0.8 0.6 0.4 0.2 C 200 400 200 400 600 1000 800 Time in samples w2 and its estimate (trembling aspen) 1200 1400 1200 1400 1 0.8 0.6 0.4 0.2 “0 600 800 Time in samples 1000 Figure 2.22: Estimates of w 1 (douglas fir) and w 2 (trembling aspen) with the EFRA based on the new prediction model. Final kappa number (softwood/hardwood) 28 26 i.hLiIàiaiI IIiIIii1ILIiIJkLhi b L 1 1i1 , L 1 ,. 1 ia .i 24 22 200 400 600 800 Time in samples H-factor 1000 1200 1400 200 400 600 800 Time in samples 1000 1200 1400 200C 1500 1000 500 “0 Figure 2.23: Process output (kappa number) and the control signal (H-factor) for a combination of one softwoods and one hardwood. 32 MMAC as a wood species compensator and identifier Chapter 2: chosen such that through the nonlinear relation, [OH] , one model can be made close to a multiple 75 of the other as: (2.66) cr This will in turn create some dependencies between the two models which can bring the excitation matrix, close to singularity, which as a consequence can affect the parameter convergence[21. In fact, this is exactly what is happening. Figure 2.24 shows the estimated w 1 at an effective alkali Estimates of Wi at an effective alkali of 15 400 600 800 Time in samples 1000 Estimates of Wi at an effective alkali of 10 0.8 0.6 0.4 0.2 200 400 600 800 Time in samples - 1000 1200 1400 Figure 2.24: Estimates of wl at different values of the effective alkali charge when a mixture of HIS is used. value of 15 and 10 respectively when the same mixture of hardwood and softwood namely trembling aspen and douglas fir is used. Clearly the first case produces closer estimates since according to the ratios in Equation 2.66 the models are less dependent: -20 (/3 2 [OHJ)n2 ’ t 11 ([OH] f 1.7 at [OH]=15 0 . 2 j at [OH] = Furthermore, looking at the condition number of the excitation matrices, 10 it is observed that indeed an effective alkali of 10 produces excitation matrices that are closer to singularities: 33 Chapter 2: Effective MMAC as a wood species compensator and iden4fier Condition number Alkak ‘i 15 23 100 10 354 225 2.1.4.d Identification of more than two species Assuming that now the incoming wood chips into the process are composed of three wood species, the corresponding mathematical models representing the extracted kappa numbers are: f ( 1 t) 1 = — 1) = K(t) — —_K ( 2 t) wi(t 1) —_1)K(t) 1) — — — (2.68) 1) — with (t) 1 4 oi R ( 2 t) = = with 1 (t 2 w — [K(t) = I 1) ‘(t)0 ( 2 t (t) 2 41 — — W3(t = — 1) t)O ( 3 t — T ( 3 K t)] (t) 2 w (t)]T 3 th wi(t) wi(t) —_‘K ( 1 t) (t 2 t 1) [K(t) [ K(t) — — z ( 3 t — — ti ( 2 t — K ( 3 t) 1) (2.70) T ( 3 K t)J 1 tii(t) (t)] 3 w (t) 2 w (t) 2 w W2(t)J - — thi(t — W3(t — (2.69) 1) — = 1 = 1 LW1(t) — — = I ( 3 t) K ( 2 t) [I K(t) (t) 2 0 and — — ‘Ki(t) 1) — (2.71) z ( 2 t — z ( 3 t — 1)K(t) 1) (2.72) 1) with (t) 3 q = [K(t) f (t) 3 6 = 1 L ( 3 t) — Ki(t) — T ( 2 K t)j (t) 1 t1 1 ( 2 zi, T t) (t) 3 th (t)] 3 w (2.73) Figures 2.25, 2.26, and 2.27 show the estimated weights when a mixture of three softwoods namely western hemlock, western red cedar, and jack pine are pulped at an effective alkali of 18. As can be seen not only do the weights exhibit a very slow convergence, but also the estimates are biased. In fact, as it wifi become clear the convergence of the weights is not guaranteed when the wood chips supply is composed of more than two species. We will show that by incorporating a third species 34 Chapter 2: MMAC as a wood species compensator and iden4fier the minimum of the performance index is not unique and therefore the condition for uniqueness of the estimates is not satisfied. This can be explained by looking at the Hatton’ s equation which is Wi and its estimate (western hemlock) - o 200 400 1000 600 800 Time in samples 1200 1400 Figure 2.25: Estimates of w 1 (western hemlock). w2 and its estimate (western red cedar) 0.8 200 400 600 800 Time In samples 1000 1200 1400 Figure 2.26: Estimates of w 2 (western red cedar). w3 and its estimate (jack pine) I I 200 400 I I I I 1000 1200 0.8 0.6 0 V 600 800 Time in samples 1400 Figure 2.27: Estimates of w 3 (jack pine). rewritten here for convenience: K = — 10 H) ([OH])’’ I3 (log . 35 (2.74) Chapter 2: MMAC as a wood species compensator and identifier Writing the Hatton’s equation in the following form: K where = A + BX 2 A 2 B = (2.75) c!j ([OH])’ = (2.76) 1 H) X = (log 2 and B . Now 2 it is evident that the Hatton’s equation is defined by only two constants namely A suppose that the Hatton’ s equations describing the chemical pulping of three individual wood species are written in accordance to Equation 2.75 as: 1 K = 1+B A X 1 (2.77a) 2 K = 2+B A X 2 (2.77b) 3 K = 3+B A X 3 (2.77c) Then we can always find two scalars a and b such that the chemical pulping kinetics of the third wood species can be written as a linear combination of the first two species as follows: =aKi+bK 3 K 2 (2.78) (aAi + bA ) + (aBi + bB 2 )X 2 where a and b can be calculated by equating the corresponding coefficients of Equations 2.77c and 2.78. With that in mind we turn our attention to the regression vectors çb for i=l,...,3 defined for three wood species in Equations 2.69, 2.71, and 2.73. Initially the analysis will be done for q and the others will follow accordingly. By substituting for K in Equation 2.69, q can be written as: c5i(t) = Q) + 2 K 1 [w Jt) w ( {+w (t) K 3 — (t) 2 K — (t)]T 3 K (2.79) or after taking N measurements as: (1)+waK K 2 wiKi(1)+w ( 3 1) (1) 2 -K —1<3(1) (2.80) wiKi(N) + w (N) + w K 2 (N) K 3 36 (N) 2 -K (N) 3 —K Chapter 2: , MMAC as a wood species compeator and idennfier by column manipulation, can be reduced to: (1) 2 —K wiKi(1) -K ( 3 1) (2.81) wiKi(N) 2 -K ( N) -K (N) 3 K as we saw earlier, can be written as a linear combination of K , Now since 3 1 and K 2 or a multiple thereof, whose third column will be linearly dependent on the other two will either not have a full rank or due the disturbances in the system be close to not having a full rank. Consequently the excitation matrix fi will be close to being singular which will affect the uniqueness of the estimated parameters[2]. The same analogy can be applied to show that 2 and fall into the same category. Although the singularity can be to some degree reduced by choosing the effective alkali and therefore the parameter B 2 in Equation 2.77 to be time-varying so that a and b in Equation 2.78 are different at every sample, it cannot be completely overcome. The reason for that arises from the strong dependencies of individual Hatton equations on one another as explained in the end of Section 2.1 .4.c. By choosing kinetic models that are more distinct based on their constants, and by choosing a value for the effective alkali that does not alter this distinction very much, the singularity of the excitation matrix can be further reduced. However, in reality, choosing a value for the effective alkali is not at user’s disposal and its value has to be chosen according to some guidelines. But here, in order to show the effect of different effective alkali values on better weight identification, hypothetical cases are simulated. Figures 2.28, 2.29, and 2.30 show the estimated weights where now the wood supply consists of one softwood, balsam fir, and two hardwoods, trembling aspen and red alder, respectively and the constant effective alkali is as before 18. As it can be seen although Wi and its estimate (balsam fir) 600 800 Time in samples Figure 2.28: Estimates of w 1 (balsam fir). 37 Chapter 2: MMAC as a wood species compensator and identjfier w2 and its estimate (trembling aspen) 0.8 0.6 — — — — — — — — — 0.4 0.2 200 400 600 800 Time in samples 1000 1200 1400 Figure 2.29: Estimates of w 2 (trembling aspen). w3 and its estimate (red alder) Time in samples Figure 2.30: Estimates of w 3 (red alder). the estimates still show bias, the parameter estimation has improved a great deal compared to the last case. Further simulations revealed that a variable effective alkali around a mean value of 18 had a very small effect on the improvement of weight identification. However, the convergence of the parameters showed an even better improvement when in an unrealistic case the variable effective alkali was drastically reduced to 2 with a standard deviation of 0.5. Figures 2.31, 2.32, and 2.33 wi and its estimate (balsam fir) Figure 2.31: Estimates of w 1 (balsam fir) at [OHJ=2. 38 Chapter 2: MMAC as a wood species compensator and identjfier w2 and its estimate (trembling aspen) Time in samples Figure 2.32: Estimates of w 2 (trembling aspen) [OHj=2. w3 and its estimate (red alder) 600 800 Time in samples Figure 2.33: Estimates of w 3 (red alder) [OHJ=2. show the estimated weights. The reason for the improvement is believed to be the contribution of the effective alkali through the relation [OH]’ to more model separation between the wood species. 2.2 MMAC versus conventional adaptive control (CAC) In the last few sections an overview of the proposed multi-model adaptive control was presented and its application to the kraft pulping particularly for wood species identification and compensation was investigated. We showed that MMAC scheme estimates the weights of the species, and the controller’s gain is calculated based on the estimated weights. In this section however, our focus will be to investigate and assess the control performance of the MIvIAC. To this end, we compare it against a conventional adaptive controller (CAC) where the controller’s gain is calculated directly based on the estimated overall process gain. The gain estimation in CAC will be accomplished by a recursive least-squares based estimator without having any knowledge of the weights associated with every species. 39 Chapter 2: MMAC as a wood species compensator and identifier To formulate our conventional adaptive controller we once again have to minimize the square of the error between the measured and the predicted kappa numbers. We start with a mathematical model describing the predicted kappa number out of the process based on kraft pulping of two species. Equation 2.82 represents such a model: (t) = —aK(t — = 1) + wbu(t — 1 — v(t 1 k) + c w — 1 — k) (2.82) 1) — where = [—K(t — 1) u(t — 1 — k) jT 1 (2.83) and ê(t) = [a (t) (2.84) describe the regression and the parameter vector respectively. The first element of the parameter vector is the process pole which need not be estimated as before assuming we have previous knowledge of it. The second term, (t) is the estimate of the overall process gain, the estimate of the overall ‘disturbance’ c , and the third term, J(t) is 1 wb Clearly the latter two terms will be time dependent since they are a function of the weights w which are varied with time. Then the minimization of the performance index will lead to a least-squares based algorithm to recursively estimate the elements of the parameter vector. The estimation algorithm used here will be once again the EFRA in order to give the CAC the same advantages as MMAC. e(t) Ô(t P(t) = P(t — 1) — = K(t) — — 1) (2.85) 1) (2.86a) + 1 + cbT(:)P(t_i)b(t)e(t) P(t—i)qS(t)b 7 , ( T t)P(t— 1) + — — 40 6P ( 2 t) (2.86b) Chapter 2: MMAC as a wood species compensator and identifier Shnulation results Again for this simulation, the only difference with the previous cases is the implementation of the new mathematical model (Equation 2.82). The wood species used for this simulation are western hemlock and western red cedar. Table 2.10 shows the EFRA constants for CAC that give the closest estimates of the process gain and disturbance. As it was mentioned before, the process pole is assumed known and need not be estimated. This is easily accomplished by setting the first row and column of the initial covariance matrix P(0) and the identity matrix I in Equation 2.86 to zero. EFRA CONSTANTS 11 11 A 1 0.01 0.90 0.001 Table 2.10 The EFRA constants for the conventional adaptive controller Figures 2.34 and 2.36 show the process gain (t) and the disturbance d(t) when they are estimated indirectly based on estimation of tI’ and th 2 via MMAC. Figures 2.35 and 2.37 show the same estimates when CAC is used to directly estimate the process gain and the disturbance. The difference Proces gain (MMAC) 200 400 600 800 Time in samples 1000 1200 Figure 2.34: Estimates of the process gain (MMAC) 41 1400 Chapter 2: MMAC as a wood species compensator and identifier Pmces gain (conventional adaptive control) 0 200 400 800 600 Time in samples 1000 1200 1400 Figure 2.35: Estimates of the process gain (CAC) Disturbance (MMAC) 125 0 lu., 200 400 800 600 Time in samples 1000 1200 1400 Figure 2.36: Estimates of the disturbance (MMAC) Disturbance (conventional adaptive control) Time in sairiples Figure 2.37: Estimates of the disturbance (CAC) is visually obvious. Clearly MMAC provides much better estimates of the gain and the disturbance. Not only is the convergence quicker but the estimates also exhibit a lower variance. The reason for this difference can be tracked back to the choice of the mathematical model. In section 2.1 .4.b we showed that when the prediction error between the measured and the predicted kappa numbers is minimized, the parameters estimated via MMAC show a large variance. However, in section 2.1.4.c 42 Chapter 2: MMAC as a wood species compensator and identifier we showed that the variance of the estimates was decreased when the prediction error minimized was between the calculated kappa number from each model and the same kappa number extracted from the measured kappa number out of the process. Now since the prediction error minimized via CAC is also that of the difference between the predicted and the measured kappa numbers, parameters estimated via CAC are expected to show higher variance than those estimated via MIvIAC. The reason for the quick convergence of the parameters in the case of MMAC is believed to be the normalization of the weights at each sampling interval; something CAC does not have the capability of doing since it has no knowledge of the weights. Figures 2.38 and 2.39 show the process output namely the kappa number, and the control signal namely the H-factor respectively. Although CAC provides good control, the mismatch between the true and the estimated gain causes the control signal and consequently the process output to have a higher variance. MMAC on the other hand exhibits lower output variance since its scheme provides good gain estimates. The difference in the output variance for both algorithms, for a number of simulations, is shown in Figure 2.40 where w 1 and 2 are varied from 0 to 1 w Final kappa number (MMAC) H-factor (MMAC) Figure 2.38: Process output (kappa number) and the control signal (H-factor), (MMAC). 43 Chapter 2: MMAC as a wood species compensator and identifier Final kappa number (conventional adaptive conftvil) 200 ‘O 400 600 800 1000 1200 1400 1000 1200 1400 Time in samples H-factor (conventional adaptive control) 13CC 1200 1100 1000 900 OUJØ 200 400 600 800 Time in samples Figure 2.39: Process output (kappa number) and the control signal (H-factor), (CAC). Output variance for CAC and MMAC — 0.7 GAG MMAC 0.6 0.5 0.4 0.3 0.2 0.1 0.2 0.3 0.4 0.5 0.6 wl(t) = 1 w2(t) 0.7 0.8 0.9 - Figure 2.40: Output variance for CAC and MMAC at different weight fractions. 44 Chapter 2: MMAC as a wood species compensator and identifier 2.3 Conclusions The structure of a multi-model adaptive controller applied to chemical pulping was presented in this chapter. It was demonstrated that the algorithm can successfully identify the wood species composition of the chip supply and take control actions to compensate for the variations resulting from a change in the composition. The structure of the proposed multi-model adaptive controller was built on the assumption that the reactor accommodating the chemical cooking of the wood chips behaves as a first order process plus a transport time delay. Furthermore, it was assumed that the designer has a prior knowledge of the process pole and the time delay. An empirical model, the Hatton’s equation, was utilized to represent the kinetics of chemical pulping of different wood species. This model, unlike the theoretical models, is easy to implement and does not require substantial computing resources to perform lengthy numerical integrations of simultaneous differential equations. Furthermore, extensive research has been done in identifying kinetic models of this type for different wood species. Because of the structure of the Hatton’s equations, the kinetic models show dependencies on one another which can affect the rate of parameter convergence. These dependencies can be suppressed by a choice of the effective alkali through the nonlinear relation [OH] in order to make the models more distinct. The algorithm’s performance is satisfactory when identifying a mixture of two species. However, due to the simplicity of the Hatton’s equations, the kinetics of a third species can always be linearly dependent of the other two or an approximate thereof. As a consequence, model dependencies are further increased and the algorithm looses its identifiability considerably when identifying three species. Although the proposed multi-model adaptive controller (MMAC) and a conventional adaptive control (CAC) both provide satisfactory control of the process, MMAC has a slight edge over CAC since it can provide a more accurate gain estimation to the controller. As a result the output variance is lower when the process control is provided by MMAC. 45 Chapter 3: MMAC as a dead-time compensator Chapter 3 MMAC as a dead-time compensator In chapter 2 a multi-model adaptive controller was introduced and its applications to wood species identification and compensation were studied. lii this chapter, a different application of multi-model adaptive control, namely dead-time compensation will be investigated. The design of controllers for processes with uncertain or time-varying parameters is a challenging problem. A fixed-parameter controller will take wrong control actions and consequently show bad performance, or even instability, if tuned based on the wrong approximation of process gain or time delay. To avoid this problem, the process time delay and dynamics should continually be tracked, parameters estimated, and the controller adjusted accordingly in order to achieve the desired performance. Two algorithms that provide such on-line information to the controller, are the Fixed Model Variable Regression Estimation in [8] proposed by Elnaggar, and the multi-model approach in [91 proposed by Gendron. Here, the two proposed algorithms will be introduced and simulated, to compare their performances in servo and regulation type problems against a fixed parameter controller. 3.1 Dead-time compensation by MMAC In [9], Gendron has proposed a method that improves the robustness of Dahlin controllers, when controlling low-order systems, in the presence of time-varying process parameters such as process time delay and gain. Assuming that the range of uncertainty in the process gain and time delay is known, the interval of uncertainty of these parameters can then be discretized in a reasonable fashion. Then assuming that the designer has prior knowledge of the process time constant, i.e. the process pole, models representing the process can be constructed. Each model assumes a discrete gain and a discrete time delay from the intervals of uncertainty where the total number of constructed models equals the product of the numbers of the discrete elements in both intervals. The true process can then be represented by a weighted sum of all the models where the weight associated with each model is estimated through an adaptation mechanism by measuring the calculated output of that model against the measured process output. In order to maintain the overall steady-state gain, the 46 Chapter 3: MMAC as a dead-time compensator weights are normalized so that their sum is one. Once there is a model available to represent the process, an adaptive Dahlin controller can be designed based on knowledge of the weights. 3.1.1 The controller Since most industrial processes can be approximated by a first-order model plus a time delay, the following model is chosen to represent the process: eTd8 HOL(S) + = where g, Td, (3.87) . and r are gain, time delay, and time constant of the process respectively. In discrete-time domain, Equation 3.87 is written in terms of its pulse transfer function as: HOL(q’) g(1 + a)1’q_k 1 + aq = Since we assume a prior knowledge of the process time constant (3.88) (r), the process pole (a) will be known. The gain and time delay of the process however, are assumed to be varying with time. Given a prior knowledge of their ranges, their interval of uncertainty can be discretized in the following manner [9]: g e [gz g] ‘ = 1, 2, ..., (3.89) ,j=1,2,...,m k] where subscripts 1 and u indicate lower and upper bounds of the uncertainties respectively, and n kE[kz and m represent the number of discrete elements in each interval. Then one possible model for the process transfer function will be: i(1+tt)q_kj H(q’) (3.90) However, since the model of the process has to account for all uncertainties the process transfer function is written as a weighted sum of Equation 3.90: n HQL(q’) m = i=1 j=1 n in _1 = z=1 )1 (1+a)q —1 1+aq’ 47 g1+a)q 1 + aq 1 fl q —k m ---.. wjjgjq i=1 j=1 (3.91) Chapter 3: MMAC as a dead-time compensator where the weights wj are normalized. Now a feedback Dahlin controller based on the process model of Equation 3.91 can be designed. However, since the process model now carries zeros in its transfer function, the design should be such that the zeros of the open-loop transfer function are not cancelled and reappear in the closed-loop transfer function as they may be poorly damped or even unstable. Q, Assuming that the desired closed-loop pole is denoted by the closed-loop transfer function of a unity gain is written as: n ,-- (1 + Q)q / —1\ 1 1Jr—1 m E E wjjgjq i=lj=I n m EZWijgj i=lj=1 Having found an expression for the process model (Equation 3.91) and the closed-loop response (Equation 3.92), the Dahlin controller for this process model is derived as: / HFBq ,j 1 HOL(q —1 ) = HcL(q) (3.93) 1—HCL(q —1 ) which reduces to: HFB(q —1 ) 1+Q n = 1+aq’ m n • m (3.94) (1 + a) 1 + Qq’ — (1 + ’ 1 Q)’ E£ i=1 2=1 3.1.2 The weight adapter Since calculation of the control signal (Equation 3.94) requires knowledge of the weights (w) some adaptation mechanism has to be introduced on-line in order to calculate the weights recursively. According to [10] a reasonable choice would be to calculate the weights based on the minimization of some norm of a prediction error. Assuming that all stationary or long-term disturbances are denoted by a slowly time-varying parameter 1’, the predicted process output is: n (t) = m wjJH(q’)u(t) + 1’ (3.95) i=1 3=1 The parameter 1’ need not be estimated if the incremental changes of input and output are considered: n m = i1 (3.96) = i=1 j=1 48 Chapter 3: 1 [k (t) 11 , MMAC as a dead-time compensator k] (t) (t) 12 j2m(t) _j m 1 (t) Figure 3.41: The predicted output of individual models where i = 1— q’ and t4jj (t) is the predicted output of each individual model generated according to Figure 3.41. Each individual prediction error is then calculated by taking the difference between the true process output and the output of each model as: ejj(t) = y(t) (3.97) — The norm of the prediction error is chosen as: v(t) (3.98) = where F(q’) is a linear filter defined as: ) 1 F(q = where 1 — q_1 (3.99) f sets the tracking speed of the algorithm, and is calculated based on the time delay uncertainty as: N,—i (3.100) where Nk is the size of the time delay uncertainty interval. Then the weight associated with every model is calculated such that it is inversely proportional to some power of its corresponding error 49 Chapter 3: MMAC as a dead-time compensator norm, and normalized as: wjj(t) (3.101) = 1Jy where a small positive value is added to v (t) to avoid division by zero. The power of the error norm sets the tracking ability of the estimates. As p approaches iiiflnity the weight that corresponds to the smallest prediction error approaches one and all other weights approach zero. 3.2 Dead-time compensation by Fixed Model Variable Regression Estimation (FMVRE) Another method for dead-time compensation, proposed by Elnaggar [8], is based on the direct on-line estimation of the process time delay provided the designer has prior knowledge of the time delay uncertainty interval. The so called Fixed-Model Variable Regression Estimator (FMVRE) provides fast time delay estimation since it does not require any knowledge of the estimated process parameters such as process pole and gain. The FMVRE makes use of an auxiliary model for time delay estimation. This auxiliary model is chosen as a first-order model plus an integer time delay since many industrial processes are over-damped and can be well approximated by a first-order model plus a time delay. Elnaggar has shown [7] that a first-order auxiliary model is enough to estimate Figure 3.42: FMVRE-based adaptive control 50 Chapter 3: MMAC as a dead-time compensator the time delay of processes even if their dynamics are represented by a higher-oder model. Once an estimate of the process time delay is available by FMVRE, a recursive least-squares based estimator can be used on-line to estimate the process parameters. This information along with the estimated time delay is then fed to the controller so that the control action can be calculated. Figure 3.42 shows a pictorial representation of a FMVRE-based adaptive controller. 3.2.1 Formulation of FMVRE Let the process under study be represented by the following discrete-time transfer function of order n: HOL(q —1 B(q’) y(t) ) = q’ + 1 b q’’ —k 1 +b q 1 + aiq 1 + + aq ... = = (3.102) ... where b÷i represents the additional zero in the presence of fractional time delay, and k = integer() (3.103) is the discrete-time integer time delay with Td being the continuous-time process time delay and T 3 the sampling period respectively. Then the process can be written in terms of its difference equation, and consequently its regression model as: y(t) where 4’ = T(t, k)8 (3.104) and 0 defined as [—y(t — 1),. . . , n —y(t — n), u(t — k — r 1), t . . . , u(t — k — n — iT (3.105) are the regression and the parameter vectors respectively. Let the predicted output be calculated by a first-order plus time delay auxiliary model of regression vector cbm and parameter vector — [_y(t — m where k E [kmin, — as: (t,A)0m = = °m 1),u(t— )]T 1 (3.106) [am, m] kmaxl the estimated time delay, is known to lie within a delay interval of uncertainty , bounded by kmin and kmax, and is estimated based on minimization of some performance index. 51 Chapter 3: MMAC as a dead4ime compensator The performance index is chosen as the expected value of the square of the prediction error between the true process output and the predicted process output: v(i) {y(t) = — q(t, A)Om} 2 N = 1 (t,k) 2 e = Efr ( 2 t,k)] (3.107) To minimize Equation 3.107 it first has to be expanded by substituting Equation 3.106 into the expression for the prediction error. (t)] + aE[y 2 (t 2 E[e ( 2 t,k)] =E[y — 2amE[y(t)y(t — — 1)] + bE [u2 (t — — 1)] 1)] (3.108) — 2bmE[y(t)u(t + 2bmE[amy(t — — 1)1 k 1)u(t — — 1)] Then by definition we can write: E[e ( 2 t,k)] = r(0) + ar(0) + br(0) — 2amry(1) (3.109) — 2bm(ruy(k + 1) — amruy(k)) where r denotes the cross-correlation function between two signals. Then the performance index can be written as: V(k) = — 2bmEi (3.110) with 0 E 1 E r(0) + ar(0) + br(0) = — + 1)— amruy(&) Evidently only E 1 is dependent on the estimated time delay 2amry(1) (3.llla) (3.lllb) k. Therefore minimization of the performance index with respect to the time delay will only affect E . According to Equation 3.110 1 minimization of the performance index means that E 1 has to be maximized if the process gain is 52 MMAC as a dead-time compensator Chapter 3: positive, and minimized for a negative gain. However, before E 1 can be optimized, the value of the auxiliary model pole am has to be known. lii [7] it has been shown that if the auxiliary model is 1 becomes: chosen to be an integrator, i.e. am = 1, the best result is achieved. Thus E 1 F + 1) = = E[y(t)u(t — E[(y(t) = E[y(t)u(t 1)] — — = — r(Ic) y(t — — k — 1))u(t — E[y(t — k — — 1)u(t — — 1)1 (3.112) 1)] 1)] = r(k+1) where z = 1 — q. Clearly step disturbances and set-point changes will not affect the time delay simulation since they will be eliminated in zy. Now optimization of E 1 can easily be performed on-line at every sample by computing the cross-correlation function between the process input and incremental changes of the process output for every time delay in the range of interest through a linear filter with a pole )‘E as: Ei(t, k ) 1 = )EE1(t — 1, k) + (y(t) (maxEi (t, k) (t,k(t)) 1 E — 1))u(t y(t — V k E [kmin, kmax] , — 1),V k 1 E [kmin, kmaxl (3.113a) for a positive gain (3.113b) = (minE ( 1 t,k) V k E [kmin, kmazl , for a negative gain Then depending on the sign of the process gain, the time delay that optimizes E 1 is chosen as the estimated time delay. 3.2.2 Estimation of the process parameters Once an estimate of the process time delay is available by FMVRE, estimating the process parameters is a straightforward problem. The estimated time delay will be used in the regression vector to change its structure at a sample. Then the resulting regression vector, (t, ) — = [_y(t — 1),..., —y(t — n),u(t — (t) — 1),... ,u(t — (t) — n )JT (3.114) and the estimated parameter vector, ê= (3.115) 53 Chapter 3: MMAC as a dead-time compensator can be used in any recursive least-squares based estimator to estimate the process parameters at that sample. 3.2.3 Implementation Steps required at every sample to implement the FMVRE algorithm are simple and are outlined below: 1. Choose initial estimates O(O) and kmin 2. Generate the corresponding regression vector (t, k), Equation 3.114. 3. Update the old parameter vector 8 using an estimation algorithm of your choice. A(O) kmax. 4. Update E 1 V k E [kmjn, kmaz], Equation 3.113a. 5. Estimate k by choosing the time delay that optimizes E , Equation 3.113b. 1 6. New sample: go to 2. 3.3 MMAC vs FMVRE in adaptive control (simulation results) In sections 3.1 and 3.2 two different methods for dead-time compensation based on the FMVRE and the MMAC were introduced. In this section these algorithms will be implemented and simulated in order to control generic processes of different orders with a time-varying time delay and gain in closed-loop. The purpose however, is to compare the performance of each adaptive algorithm against a fixed parameter Dahlin controller in servo type problems. Since the MMAC algorithm is designed to control low-order systems, the comparison will be done only for first-order and second-order systems. 3.3.1 First-order systems The process under control is a first-order plus a time delay given by Equation 3.116: y(t) ay(t — 1) + g(t)(l + a)u(t — k(t) — 1) + n(t) (3.116) where g(t) and k(t) are the time-varying gain and time delay respectively, and n(t) is a normally distributed random sequence with zero mean and variance a 2 =0.05. The reference signal is chosen as a square wave with a period of 100 samples, a mean of zero, and a peak-to-peak amplitude of 54 Chapter 3: MMAC as a dead-time compensator 2. Since the MMAC assumes a prior knowledge of the process open-loop pole, this parameter will be assumed a known constant and will not be estimated throughout the simulation. However, the remaining process parameters such as time delay and gain will be varied with time according to Table 3.11. The process is controlled by a Dablin controller that sets the closed-loop pole to 0.5 from its open-loop value of 0.75. The time delay interval of uncertainty for both algorithms is assumed to be Time ft) Dead-time, kft Gain, g(t) Pole (a) 0<1<300 7 I 0.75 300t<600 5 2 0.75 t 2 5 0.75 600 1000 Table 3.11 Time-varying parameters of the simulated process. [1,8] with increments of 1. The gain interval of uncertainty for the MMAC algorithm is assumed to be [0.5,5.5] with increments of 0.5. The power p of the error norm in equation 3.101 was set to 10 since it gave the best simulation results. The gain estimator used along with FMVRE is a recursive least-squares estimator (RLS) with exponential forgetting. In the first simulation the process time delay was kept constant and the gain was varied. Figure 3.43 shows the process response when a fixed-parameter Dahlin controller is used. Clearly in the absence of adaptation the controller performs poorly. For the first 600 samples, as it is shown, the fixed gain is higher than the true gain of the process. Therefore since the controller’s gain will be smaller than what it should be, a sluggish behavior is expected from the controller. After 600 samples however, the fixed gain is lower than the true gain. As a result the controller will have a higher gain than expected which will result in oscillations in the process response. Figure 3.44 shows the same response when adaptation is provided by MMAC. Evidently MMAC frees the process from the severe oscillations observed before, and provides a much better performance. However, the controller again shows a sluggish performance before 600 samples and some minor oscillations after 600 samples. This can be explained by looking at Figure 3.45 which displays the time delay and gain estimated by MMAC. It should be mentioned that although MMAC does not directly estimate the time delay, a measure of the time delay estimates can be plotted by plotting the time delays that correspond to the largest computed weights. Looking at the estimates of the gain it is apparent that, 55 Chapter 3: MMAC as a dead-time compensator Process output (fixed-parameter controller) 400 600 500 Time in samples Control signal (fixed-parameter controller) 300 400 600 500 700 Time in samples Fixed gain (fixed-parameter controller) 1000 S Fixedgain True gain 6 4 2 I I 100 200 300 400 600 500 Thne in samples 700 800 900 1000 Figure 3.43: Process response, control signal, and gain variations (no adaptation). Process output (MMAC) 100 200 300 400 600 500 Time in samples Control signal (MMAC) 700 800 900 10 00 100 200 300 400 500 600 Time in samples 700 800 900 1000 Figure 3.44: Process response and control signal (adaptation by MMAC). at low excitations, the estimates tend to converge to some median that lies in the centre of the gain uncertainty interval. This provides an advantage, analogous to parameter estimation with ‘leakage’, since in the steady-state MIvIAC will prevent the estimated gain from drifting. However, since the median lies above the tme gain before 600 samples, the controller as expected shows sluggish performance due to overestimation of the gain. The minor oscillations in the response observed 56 Chapter 3: MMAC as a dead-time compensator Estimated delay (MMAC) 8 LI 6 }v 100 6 6 200 ! ‘1 ‘‘ 300 400 600 500 Time In samples Estimated gain (MMAC) 700 800 900 1000 300 400 500 600 Time in samples 700 800 900 10 00 I I Estimated gain 100 200 Figure 3.45: Process time delay and gain estimated by MMAC. Process output (FMVRE/RLS) ‘ml 100 200 300 400 500 600 Time in samples Control signal (FMVREIRLS) 700 800 900 10 00 Time in samples Figure 3.46: Process response and control signal (adaptation by FMVRE/RLS). at set-point changes are due to the gain underestimation. At the steady-state the control action is calculated based on an underestimated gain. Therefore at the instant when the set-point changes the response will shoot beyond its reference. Due to high excitations at this instant, the MMAC is able to track the time-varying parameter correctly and feed the correct gain estimates to the controller. However, the correct gain is fed to the controller after the oscillations has already occurred and all the controller can do now is to stabilize the response. Figure 3.46 shows the same response when a combination of FMVRE and RLS are used to estimate the time delay and the gain respectively. As seen in this figure, the FMVREIRLS provides a better closed-loop performance than that of the MMAC. The reason for that is the direct and correct gain estimation shown in Figure 3.47. The 57 Chapter 3: MMAC as a dead-time compensator Estimated delay (FMVRE) I I I I 6 4 2 100 200 300 400 500 600 Time in samples Estimated gain (RLS) 700 800 900 10 DO 400 500 600 Time in samples 700 800 900 1000 Estimated gain 6 100 200 300 Figure 3.47: Process time delay and gain estimated by FMVRE/RLS. Process output (fixed-parameter controller) 100 5 200 300 400 500 600 700 Time in samples Control signal (fixed-parameter controllar) I I I I 800 I I 700 800 900 1000 900 1000 6 4 Fixed delay Truedelay 2 100 200 300 400 500 600 Time in samples Figure 3.48: Process response, control signal, and time delay variations (no adaptation). overshoots observed in the response are again due to the gain mismatch between the true and the estimated gain. However, unlike the MMAC response, they only occur at the time of abrupt changes of the true gain. Next, the process gain is kept constant and the process time delay is varied with time according to Table 3.11. Again in the absence of adaptation the fixed-parameter Dablin controller performs very poorly in the presence of time-varying time delay as shown in Figure 3.48. However, 58 Chapter 3: MMAC as a dead-time compensator both MMAC and FMVREIRLS show good performance, handling the time delay variations as shown in Figures 3.49 to 3.52. The MMAC shows good robustness properties and closed-loop control performance although it does not estimate the time delay explicitly. FMVREIRLS on the other hand achieves the same task by providing exact estimates of the process time delay and gain. Process output (MMAC) 0 —1 nnnnnuu 100 200 300 400 500 600 Time in samples Control signal (MMAC) 700 800 900 l( Time in samples Figure 3.49: Process response and control signal (adaptation by MMAC). 100 200 300 100 200 300 400 500 600 Time in samples Estimated gain (MMAC) 700 800 900 1000 700 800 900 1000 6 400 500 600 Time in samples Figure 3.50: Process time delay and gain estimated by MMAC. 59 Chapter 3: MMAC as a dead-time compensator Process output (FMVRE/RLS) 100 200 300 400 600 500 Time in samples Control signal (FMVRE/RLS) 700 800 900 1C 400 600 500 Time in samples 700 800 900 1000 Figure 3.51: Process response and control signal (adaptation by FMVRE/RLS). Estimated delay (FMVRE) I I 100 200 8 300 I 400 600 500 Time in samples Estimated gain (RLS) I I I I I 700 800 900 I I I 1000 6 IIIlIIl 100 200 300 400 600 500 Time in samples 700 800 900 1000 Figure 3.52: Process time delay and gain estimated by FMVRE/RLS. 3.3.2 Second-order systems As it was shown in section 3.1, the MMAC algorithm assumes a first-order model of the process under control regardless of the process order. If the actual process has an order greater than one, the Dahlin controller will be designed with mismatches between the actual process and its first-order model. This model uncertainty will in turn degrade the performance of MMAC since Dahlin controllers will perform well when an accurate model of the actual process exists [5]. The FMVRE/RLS algorithm on the other hand, does not suffer from the same limitations in choosing a 60 Chapter 3: MMAC as a dead-time compensator model for the actual process. Since it has an on-line recursive least-squares estimator, it can always estimate the parameters of a higher order model so that the controller can be designed based on a more accurate model of the process. The difference between the two algorithms will become more apparent by simulating a second-order process in closed-loop. Assuming the process under control has a second-order continuous-time transfer function given by Equation 3.117 where g and Td are the time-varying gain and time delay respectively, off-line HOL(S) = g e_Td8 2 (3.117) (s+1) identification was performed at a sampling period of T =1 second, in order to identify a first-order 8 discrete model for the process transfer function so that it can be used in the MIvIAC algorithm. Figure 3.53 shows the estimated open-loop parameters when the reference signal is a square wave with a period of 50 samples, and the process time delay and gain are 5 and 1 respectively. The reason for Estimated parameters of a first-order model Gain= 0.26 0.5 0 100 Pole=- 0.77 if11ll11 lit! I,Irl1rtt 1.1, 200 400 300 500 600 700 Time in samples Estimated time delay of a first-order model 100 200 -0.5 rlj ‘0 ,i rrrl 800 r., 900 800 900 , ‘C DO 5 6 4 2 n 300 400 500 600 Time in samples 700 1000 Figure 3.53: Open-loop identification of a first-order model parameters. the drift in the estimated pole at the down-going or up-going edge of the reference signal is due to the high frequency content of the reference signal at these instants which will excite the high order dynamics of the plant when in fact the estimated model is only a first-order model. Further simulation revealed that variations in the gain and time delay do not affect the estimated pole of the 61 MMAC as a dead-time compensator Chapter 3: first-order model in anyway. Therefore since the MMAC requires a known constant pole before it can be implemented, the pole was chosen as estimated, leading to the following discrete-time model: H0L(q’) q_k 1 O26 1 g (3.118) Having approximated a first-order model for the second-order process, the MMAC algorithm is easily implemented as before. The process gain and time delay are again varied according to Table 3.11. Figures 3.54 and 3.55 show the results when MMAC is used to control the second-order process of Equation 3.117 based on the approximated first-order model of Equation 3.118. Clearly MMAC Process output (MMAC) 50 100 150 200 250 300 Time in samples Control signal (MMAC) 350 400 450 500 50 100 150 200 250 300 Time In samples 350 400 450 500 Figure 3.54: Control of a second-order system (adaptation by MMAC). Estimated delay (MMAC) 1 A 50 A: I 50 100 I 100 150 I 150 200 250 300 Thne in samples Estimated gain (MMAC) I I I 200 250 300 Time in samples 400 450 I I I 350 400 450 350 Figure 3.55: Process time delay and gain estimated by MMAC. 62 5( 500 Chapter 3: MMAC as a dead-time compensator still stabilizes the process as the algorithm tracks variations of time delay and gain. However, the mismatch between the second-order actual process and the approximated first-order model will affect the performance of the Dahlin controller which is evident in Figure 3.54. Further simulations reveal that the FMVREIRLS algorithm produces the same results when estimating a first-order model of the process. As it was mentioned earlier Dablin controllers require an accurate model of the process under control. For processes described by transfer functions of higher orders than one, the MMAC’s performance is limited since it assumes a first-order model of the process. The FMVRE/RLS algorithm however, does not suffer from these limitations. The versatility of the on-line RLS estimator allows the algorithm to estimate parameters of a higher order model by means of simple computations. The MMAC algorithm however, requires an enormous number of computations at every sample which grows as the number of models are increased. Figures 3.56 and 3.57 show the results when Process output (FMVRE/RLS) Time in samples Conirol signal (FMVRE/RLS) Figure 3.56: Control of a second-order system (adaptation by FMVRE/RLS). FMVRE/RLS is utilized to control the same process when time delay and gain are varied. The on-line RLS estimator assumes a discrete-time second-order model of the form displayed in Equation 3.119 once it receives a time delay estimate, k from the FMVRE algorithm. b , b 1 , &, and 2 — ) 1 H0L(q — q_ b 1 q 2 +b q_2 2 1+aiq_1+a a2 are (3.119) the parameters of the second-order model estimated by the RLS estimator. It can be seen that for 63 Chapter 3: MMAC as a dead-time compensator controlling processes of higher order the FMVRE /RLS algorithm is superior. The adaptive Dahlin controller shows good performance since it is designed based on a second-order discrete-time model whose parameters are estimated recursively at every sample as shown in Figure 3.57. Estimated delay (FMVRE) 50 100 150 200 250 300 Time in samples Estimated parameters (RL5) 350 a, 50 100 150 200 250 300 Tin in samples 350 400 450 400 450 - 500 Figure 3.57: Process time delay and parameters estimated by FMVRE/RLS. Up until now, in Section 3.3, the two discussed methods of dead-time compensation have been compared within the frame of servo problems where the objective is to have the process output respond to changes in the set-point in a specified manner. In the following section however, we turn our attention to the other control task called regulation where the objective is to regulate the process output around a fixed set-point despite the presence of disturbances. Since the MMAC algorithm has been successfully implemented on a chlorination tower[lO], in order to compare the control performance in regulation, the prebleaching stage in a bleach plant will be simulated, identified, and then controlled in closed-loop by a Dahlin controller where the adaptation will be provided in turn by the two algorithms. 3.4 Adaptive control of a chlorination tower by MMAC and FMVREIRLS (an industrial example) Figure 3.58 shows a chlorination tower which is the first stage in a multi-stage bleaching sequence. Here, prebleaching of pulp is performed with a mixture of chlorine and chlorine dioxide substitution since the presence of chlorine dioxide influences the viscosity of pulp and is less harmful to the 64 Chapter 3: MMAC as a dead-time compensator Bleached Pulp Sensor 2 Unbleached J Pulp EEEME Mixer I Tower Figure 3.58: Prebleaching stage with a mixture of chlorine and chlorine dioxide. environment [13]. The percentage of chlorine dioxide substitution is varied from time to time but is usually set to 50%. In what comes next the equations relating kinetics and dynamics of the kappa number out of a chlorination tower to the flow of chlorine and chlorine dioxide will be derived. Subsequently the process will be simulated, and identified in open-loop. Finally, assuming that a kappa number sensor is available, adaptive Dahlin controllers will be implemented based on the MMAC and FMVREIRLS algorithms to regulate the prebleaching kappa number in a closed-loop. 3.4.1 Kinetics and dynamics of chlorine and chlorine dioxide prebleaching Assuming that prebleaching can be modelled by a combination of ideal mixing represented by a continuous stirred tank reactor (CSTR), and a plug flow represented by a pure time delay, we start by writing the material balance equations for the CSTR when only chlorine is used in pulp treatment, and extend it to bleaching with chlorine and chlorine dioxide. The material balance equations are [1]: dM = d(MK) dt F v(Mo — F = — 0 (C — C) M) MK) — (3.120a) — pMsr Mr (3.120b) (3.120c) where M=pulp consistency/100, K=kappa number, C=chlorine concentration (gmIL), Fpulp-water slurry flow rate (L/s), V=volume (L), r=reaction rate, p=slurry density (gmIL), s=stoichiometric 65 Chapter 3: MMAC as a dead-time compensator factor, t=time (s), and the subscript o indicates initial values at the inlet valves. For Equation 3.120 to accommodate the addition of chlorine dioxide, the reaction rate, r, and the stoichiometric factor, s, will have to be different, and the chemical concentration, C, will have to represent the sum of the concentrations of the two bleaching agents. The reaction rate for a mixture of chlorine and chlorine dioxide is [13]: = = [ClO 1 —k [ ] 2 K 5 Cl ° where K=kappa number, k =rate constant, and 1 (3.121) [ ] denotes the concentrations of chlorine and chlorine dioxide in gmIL. However, Equation 3.121 requires the individual concentrations of both bleaching agents at time t which will have to be extracted from Equation 3. 120c which now represents the sum of the concentrations of both chlorine and chlorine dioxide. In order to derive such expressions, we take advantage of the stoichiometry and the relative rates of consumption of chlorine and chlorine dioxide during bleaching given as [11, 12]: = 2 + Cl —(zClO ) 2 i.(LCl ) 2 ) 2 QClO = [H]°° 2 k K 8 ’ / ] 2 [Cl aV 2 [do ] — — (3.122a) (3 122b) where s is the new stoichiometric factor, [H+] is the concentration of the hydrogen ions, a and k 2 are constants, is the differential operator, and LCl 2 and i.ClO 2 are the total chlorine and chlorine dioxide consumption in gmlgm pulp during a given bleaching time t. By rewriting Equations 3. 122a and 3.122b we have: = 2 + zCl —(zClO ) 2 a/[] (1?ClO ) 2 (3.123a) (3.123b) Solving for i..Cl 2 and zClO 2 we arrive at the following expressions for the total chlorine and chlorine dioxide consumptions at any given time t: ) 2 (zCl (3.124a) = + 66 aV ] 2 [do Chapter 3: MMAC as a dead-time compensator __(sr) = 1 (3.124b) +aV [do ] 2 Expressions for the individual concentrations of each bleaching agent at time t can be easily found by first rewriting Equation 3.120c in terms of chlorine and chlorine dioxide concentrations as: 2 + dO (Cl ) 2 = 20 + C10 (Cl 20 2 Cl — — ) + pMsr 2 dO (3.125) , and C10 2 20 are the concentrations at the inlet valves. Then by , Cl 0 where now, instead of C substituting Equation 3. 123a into Equation 3.125 we have: ) 2 + dO = , C10 , j(Cl2, + 2 2 Cl — — ) 2 dO — ) 2 pM(i.ClO2 + zSCl (3.126) Expressions for the individual change in the concentration rate of each bleaching agent can be written by separating Equation 3.126 in the following fashion: ) 2 (Cl ) 2 ClO = = 20 (Cl 20 (ClO — — ) 2 Cl — ) 2 C10 ) 2 pM(zCl — ) 2 pM(LiClO (3.127a) (3.127b) Finally by substituting Equations 3.124 into 3.127, rates of change in chlorine and chlorine dioxide concentrations can be represented individually in terms of the concentrations of both agents at time t: ) 2 -(Cl ) 2 ClO = = 20 -(Cl 20 ClO — 21(pMsr) )+ 2 Cl — )+ 2 dO 1 (pMsr) 1 _ (3.128a) (3.128b) + The new material balance equations describing the prebleaching process with chlorine and chlorine dioxide are: = 0 (M 67 — M) (3.129a) Chapter 3: d(MK) = ) 2 (Cl (ClO 4 ) 2 = = K 0 ç(M 20 (Cl — 20 (ClO — MK) )+ 2 Cl — — MMAC as a dead-time compensator Mr (pMsr) 1 (pMsr) 1 , )+ 2 dO (3.129b) (3.129c) (3.129d) + = [ClO 1 —k ] [ J 2 K 5 Cl ° ° (3.129e) K’ 8 2 [H+] °° k (3.1290 = Finally in order to represent the plug flow in the chlorination tower, the measured process variables such as consistency, kappa number, and bleaching agents concentrations are delayed by a time delay equal to their residence time in the chlorination tower until they reach a kappa number sensor. From Equation 3.129 it is clear that the effect of chlorine and chlorine dioxide on the outcoming kappa number is a non-linear relationship. However, as it will be shown in the next section, the prebleaching dynamics are well damped, stable, and resemble a first-order process plus a time delay. 3.4.2 Process identification In this section the prebleaching process governed by the system of differential equations derived previously will be simulated. The attempt will then be to utilize an identification method in order to estimate a transfer function relating the kappa number to the flow of the bleaching agents by a first-order system plus a time delay. The pole of the estimated first-order model will later be used as a fixed parameter in closed-loop control. Here, the method of open-loop identification chosen will be FMVRE/RLS algorithm. Figure 3.59 shows the static input/output relationship of the process found under steady state conditions. As expected the relation is a nonlinear one, and as a result the process gain can vary drastically at different operating points; therefore before the process can be identified a desired operating region has to be chosen in order to ensure correct parameter identification. This region is depicted in Figure 3.59 as the shaded area. According to [25] the kappa number out of a chlorination tower for kraft pulp is around 10, which unless mentioned otherwise, will be used here as the 68 Chapter 3: MMAC as a dead-time compensator Static inputloutput relationthip I 15 a2+c1o2 (g/L) Figure 3.59: Static inputloutput relationship of the process operating point in all simulations. Figure 3.60 shows the process response to a change in the flow of chlorine and chlorine dioxide assuming that a kappa number sensor is located 6 minutes away from the chlorine and chlorine dioxide valves. Clearly the response of the process, being over-damped, resembles a first-order dynamics plus a time delay. The parameters of such a model can easily be calculated from the response as shown in Figure 3.60, leading to the following continuous-time first-order approximation: = HOL(S) OOS3e where the time delay and time constant are given in minutes. (3.130) However, before the open-loop ‘approximated model’ of Equation 3.130 can be utilized in order to design an adaptive kappa number controller, we have to ensure that the approximated model is consistent with the ‘estimated model’ calculated by the estimator recursively. To have the estimator estimate a first-order model of the process, the existing model uncertainty between the actual process and the estimated model has to be taken into account. Since the actual process might contain higher order dynamics, we have to ensure that these dynamics are not perturbed by the high frequency content of the incoming control signal since in this case the mismatch between process and the model will be severe. To get a feel for the model mismatch it is convenient to compare the Bode magnitude diagrams of the process and its approximated model. Since no generic transfer function of the actual process was available, the Bode 69 Chapter 3: MMAC as a dead-time compensator Step increase in chlorine and chlorine dioxide flow 3( I I I I 25 20 15 1(:90 100 110 11 120 I --I 130 140 Time (mm) Kappa number I 150 160 170 18 0 160 170 180 I De1a°” io.: 0083 5 Z 1 I l\\G . Constant 10 tfllfl nç I 100 110 120 130 140 Time (mm) 150 Figure 3.60: Prebleaching response to a change in chlorine and chlorine dioxide flow. diagram in this case was produced by injecting a sinusoid of peak-to-peak amplitude of 10 and a mean value of 20 at various frequencies. Figure 3.61 shows the results. Evidently the approximated model of equation 3.130 provides a very good first-order match at frequencies lower than 0.002 radls. Therefore for good parameter estimation of a first-order model, any mismatches at higher frequencies have to be avoided by lowpass filtering of the signals required for process identification at some cutoff frequency of w<O.0O2 radls. However, the lower bound of the excitation signal frequency required to correctly estimate the process time delay by FMVRE is w raclls where T 3 and r are the sampling period and the fastest time constant of the process respectively [7]. At a sampling period of =3 minutes and a time constant of r=10 minutes the lower bound of the excitation signal frequency 8 T is calculated to be 0.004 rad/s. This value falls beyond the cutoff frequency of the filter. As a result parameter estimation will be conducted utilizing the filtered signals, and time delay estimation using the raw signals. Figure 3.62 shows a pictorial representation of the proposed identification scheme where LPF blocks denote lowpass filtering. In order to simulate the proposed scheme and identify the process, white noise superimposed on a step function of amplitude 20 was injected into the system, and the process was sampled at a sampling period of T =3 minutes. The lowpass filters chosen 3 70 Chapter 3: MMAC as a dead-time compensator here were digital Butterworth filters of order 2 and cutoff frequencies of w=O.OOl rad/s. Figures 3.63 and 3.64 show the estimated parameters of a first-order model before and after filter installation respectively. Table 3.12 shows how the estimated models compare against the approximated model. Bode magnitude plot -25 -30 -35 — — • . . • . . • . . Trmcess Apprmomated model — . — -50 -55 10 1(12 1(i 10-I Frequency (radls) Figure 3.61: Bode magnitude plot of the true process and its approximated 1st-order model. Process Figure 3.62: Process identification based on a first-order model approximation. The first column displays the approximated model of Equation 3.130 in the discrete-time domain when 71 Chapter 3: Estimated pole Estimated gain -0.02 I I I I I I — ill. -0.2 -0.04 —- -0.4 +— —— — -+. ——— —+— — —-+--————+ I-————4———— -0.6 -0.06 -0.08 MMAC as a dead-time compensator --1 -0.8 I I I -0. 100 I .1. 2K) 4 ) 300 Time in samples Estimated time delay 51 K) 100 400 200 300 Time in samples 500 Approximated model .L. — —— - 4 —0.022’ 1— O74q” HOL(q ) 1 — 3 — —I— Estimated model: H0L(q ) 1 = 200 300 400 Time in samples 100 —0. 016q q 500 Figure 3.63: Estimated parameters of a first-order model before filter installation. Estimated gain Estimated pole U : ?4 -0.06 -0.08 •-—-i I I —01 •100 -0.2 I I I -I I I I I I t . +I I I I I I I I I I +----I I I 200 300 400 Time in samples Estimated time delay 500 — 100 200 400 300 Time in samples 500 Approximated model: £ 4 — L(q ‘ 0 IJ 3 = —0.022q-’ 1— -2 — Estimated model: ) I = I Olq_2 2 O.O 1— O.77r’ — 100 200 400 300 Time In samples 500 Figure 3.64: Estimated parameters of a first-order model after filter installation. sampled at T =3 minutes. The second and third columns show the estimated models at the same 3 sampling period produced by open-loop identification with and without filtering respectively. It is clear that correct estimates of a first-order discrete model parameters are only achieved when filtering is in effect since the parameters of this model are closer to those of the approximated model. 72 Chapter 3: MMAC as a dead-time compensator Table 3.12 A comparison between the approximated and the estimated models. It is worth mentioning that although in open-loop, FMVRE/RLS is capable of estimating a more accurate discrete model of a higher order for the process in the absence of filtering and under high excitation conditions, it was chosen not to do so since here the closed-loop control objective is to regulate the prebleaching process around a fixed reference point. In a regulation problem, due to constant value of the reference signal, the excitations in the signals required for parameter identification are only provided randomly by the presence of the disturbances in the system. In bleaching processes the disturbances acting on the systems vary with time due to the variations in the pulp quality such as pulp grade changes and variations in pulp residual lignin content. These random disturbances which are nonstationary due to their nature do not have a strong enough frequency content in order to keep the estimator active at all times. Therefore in closed-loop, by choosing to estimate the parameters of a higher order model, the chances of estimating a higher number of parameters correctly are further reduced. 3.4.3 Noise model Figure 3.65 shows a block diagram representation of the process where now a noise model has been included. n(t) is an independent normally distributed white noise with a zero mean and a variance a and V disturbances. = 1/(1 — q’) is an integrator which will simulate the effects of nonstationary Starting with the first-order discrete model estimated in the previous section, the n(t) y(t) u(t) Figure 3.65: Process open-loop block diagram with a noise model. 73 Chapter 3: MMAC as a dead-time compensator predicted model of the process in terms of its difference equation model is: (t) = 1 q u(t) + Vn(t) = —ay(t (3.131) where and k — 1) + (1 + a)u(t — k — i) + V(n(t) + an(t — 1)) are the process gain and time delay to be estimated, and a is the open-loop pole of the process always fixed at its estimated open-loop value. Note that although n(t) in Figure 3.65 is a white noise, its overall effect on the process output according to Equation 3.131 is colored since n(t) and n(t — 1) are correlated. Consequently the least-squares estimate will be biased when applied to the predicted model of Equation 3.131. However, one way to avoid this difficulty is to sample the process at widely spaced time intervals[2]. For this reason the sampling period of 3 minutes used earlier seems sensible. 3.4.4 Closed-loop control issues In section 3.4.2 a method was proposed to correctly and recursively estimate a first-order discrete model of a higher-order prebleaching process in open-loop. In this section the attempt will be to close the loop and to regulate the process around a fixed set-point in an adaptive mode via a Dahlin controller. To achieve this task the proposed model approximation method will be used in the FMVRE/RLS algorithm to achieve correct model estimation and therefore better closed-loop control. The results will then be compared against a Dahlin controller whose adaptation is provided by the MMAC algorithm. The proposed on-line method of identification which will be utilized in the FMVREIRLS algorithm will guarantee that a first-order model can be estimated by filtering the high frequency contents of the signals required for identification before they enter the estimator. However, if at any instant in time the control signal entering the process has higher frequencies than the cutoff frequency of the filter, the higher order dynamics of the system will be perturbed, and as a result the mismatch between the actual process and its estimate grows which can lead to instabilities. One method to overcome this problem and make the controller more robust is to choose the closed-loop time constant of the system such that its bandwidth is equal or less than the filter’s cutoff frequency. This way we guarantee that the process input and output will not contain frequencies that will excite higher 74 Chapter 3: MMAC as a dead-time compensator order dynamics of the process. As a result the closed-loop pole of the system will have to be greater than its open-loop value which can be justified by considering the fact that the control objective for this particular problem is not to have the kappa number track a specific moving reference signal at a certain speed but to have it stay close to a specific constant value. As it was mentioned earlier the only excitation introduced to the system in closed-loop is random and provided by nonstationary disturbances since the set-point will be constant. As a result the estimator will not be very active, and because of the reasons explained in chapter 2, the growth of the covariance matrix trace of the estimated parameters in the presence of very low excitations will lead to parameter estimate’s drift from their true values. Therefore the on-line estimator used here in the FMVREIRLS algorithm will not only be of the EFRA type but also will feature ‘leakage’ [24]. The term leakage is refeffed to the amount of update that is ‘leaked’ away from the new estimates in the absence of excitations. Ultimately, in the presence of no excitations, the overall leakage amount will be such that the parameter estimates will converge to a user specified set of parameters. The equation representing the EFRA parameter update with leakage is: Ô(t) where ( = (1- Ô(t -1) + (1- + (0* 1+ (3.132) sets the amount of leakage and 9* is a set of parameters specified by the user that will give reasonable control. It is clear from the parameter update equation that when there is no new information coming into the estimator, i.e., e(t)=O and ö(t)=ê(t Ô(t) = — 1), Equation 3.132 will reduce to 0’, and as a result the parameters are bounded. To simulate and control the process in closed-loop, the value of the leakage parameter (, in the FMVRE/RLS algorithm, is set to 0.01. The open-loop pole of the process is assumed to be known and for a sampling period of 3 minutes is kept at a constant value of -0.77. Therefore the only parameter to be estimated is the process gain. The time delay interval of uncertainty for both algorithms is assumed to be [1,8] with increments of 1. For the MMAC algorithm the gain interval of uncertainty is assumed to be [-0.01 ,-0. 1] with decrements of -0.01, and the parameter p is set to 10. The white noise sequence, n(t), has a zero mean and a variance 0.003. The Dablin controller 75 Chapter 3: MMAC as a dead-time compensator sets the closed-loop pole equal to the open-loop pole in order to keep the closed-loop bandwidth near the cut-off frequency of the lowpass filters for the reasons explained in section 3.4.2. 3.4.5 Simulation results Figures 3.66 through 3.69 show the prebleaching kappa number, the flow of bleaching agents, and the estimated process parameters under constant time delay and constant process gain conditions when control is provided by an adaptive Dahlin controller, and adaptation is provided by the FMVREIRLS and the MIvIAC schemes. The kappa number target was decreased after 300 samples and remained unchanged for a period of 100 samples before it was brought back to its original value. Both schemes showed a good closed-loop response to set-point changes. The overall nonstationarity behavior observed in the flow of the bleaching agents is due to the compensation of the nonstationary disturbances entering the process in order to keep the kappa number at its target value. In the next simulation the process gain was increased by 10% after 800 samples and remained unchanged for a period of 100 samples before it was brought back to its original value. Figure 3.70 shows the performance of a fixed-parameter Dahlin controller. Clearly under no adaptation the controller is not able to stabilize the kappa number due to the mismatch between the true process gain and the gain controller has assumed. A fixed-parameter Dahlin controller will assume a gain that is an underestimation of the true gain. As a result the controller’s gain will be higher than what it should be leading to an oscillatory behavior as shown in Figure 3.70. Figures 3.71 and 3.72 show the results under the same circumstances when the Dahlin controllers are adaptive, and the adaptations are provided by the FMVREIRLS and the MMAC algorithms respectively. Evidently both adaptive algorithms overcome the instabilities encountered in the absence of adaptation since now the control signals respond to a change in the process gain in order to compensate for it. Figures 3.73 and 3.74 show the parameters estimated by the two algorithms. The process gain estimated by the FMVRE/RLS algorithm deviates from its average value at 800th sample in order to track the new process gain and returns to its original value after the 900th sample. This tracking of the new process gain is more obvious for the FMVREIRLS algorithm since the leakage parameter ( chosen in the estimator has a low value, and as a result the estimates have more flexibility to follow new changes. 76 Chapter 3: MMAC as a dead-time compensator As for the MMAC algorithm the tracking of the new process gain between the 800th and the 900th samples is not as clear although it deviates from its average value towards a higher estimate. This is believed to be a result of the strong leakage which is a characteristic of the MMAC algorithm. Figure 3.75 shows the performance of a fixed-parameter Dahlin controller when there is a mismatch, for the first 500 samples, between the true time delay of 5 samples and the time delay of 2 samples that the controller has assumed. Although in reality it is very unlikely that the time delay would have a large variation between 5 and 2 samples at a sampling period of 3 minutes, this variation has been chosen to make the comparison more visual. The oscillations observed are due to the fixed-parameter controller’s underestimation of the true time delay. Figures 3.76 and 3.77 show the same results when the Dahlin controllers are adaptive via FMVREIRLS and MMAC algorithms. Evidently both algorithms stabilize the process since they are capable of tracking the variations of the time delay. Figures 3.78 and 3.79 show the estimated time delay associated with the two algorithms. Clearly the FMVRE/RLS algorithm provides the control with the correct and direct estimation of the process time delay. For the MMAC algorithm in Figure 3.79 a graphical measure of the estimated time delay is achieved by plotting the time delay corresponding to the highest estimated weight associated with a predicted model. It can be seen that although this is not an exact estimate of the time delay since in reality no weight estimates are zero, the time delay estimates tend to average around the true time delay which is a sign of adaptation. In the final simulation, performances of the two adaptive algorithms are compared under tighter closed-loop control requirements against one another. Figures 3.80 through 3.83 show the process output, the control signal, the estimated time delay, and the estimated process gain when the two algorithms are in turn employed for adaptation. The closed-loop pole is set to its previous value of -0.77 implying a loose control. Under these conditions both adaptation algorithms have acceptable performances. However, a difference is observed when the value of the closed-loop pole is gradually reduced. Although in the process industry, the trend in regulating a process output in the presence of uncertainties and disturbances is to choose a process closed-loop pole that is larger than its open-loop value in order to increase robustness and therefore reach acceptable control, here we have chosen to reduce the closed-loop pole for the sake of a comparison in order to find the restrictions on 77 Chapter 3: MMAC as a dead-time compensator the achievable performance of both adaptive algorithms. Figures 3.84 through 3.87 show the same simulation results when the value of the closed-loop pole is set to -0.4 implying a tight control. It is observed that closed-loop control based on FMVRE/RLS algorithm is more susceptible to model uncertainties. The reason for this sensitivity is believed to be due to an increase in the closed-loop bandwidth of the system. This allows a higher frequency control signal to enter the process and to excite higher process dynamics when the FMVRE/RLS algorithm only estimates, through lowpass filtering, a first-order model. This creates a mismatch between the model and the true process which leads to instabilities. 78 Chapter 3: MMAC as a dead-time compensator Kappa number 10 V-’- Th I 100 200 I 300 400 500 600 Time in samples I I I 700 800 900 10( C12+C102 (g/L) bC I I I 80 60 00 200 300 400 500 600 Time in samples 700 800 900 1000 Figure 3.66: Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by FMVRE/RLS. Estimated time delay E p p p I 6 4 100 200 300 400 500 600 Time in samples 700 800 900 10( 700 800 900 1000 Estimated gain -0.08 200 300 400 600 500 Time in samples Figure 3.67: Estimated time delay and gain of the prebleaching process by FMVRE/RLS algorithm. 79 Chapter 3: MMAC as a dead-time compensator Kappa number —-v. 10 0 I I 200 300 I 400 500 600 700 800 900 1000 700 800 900 1000 Time in samples C12+C102 (g/L) 80 60 40 20 I, TOO 200 300 400 500 600 Time in samples Figure 3.68: Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by MMAC. Estimated time delay 6 4 2 F DO 200 300 400 500 600 Time in samples 700 800 900 10 )0 700 800 900 1000 Estimated gain -0.02 -0.04 -0.06 -0.08 1O0 200 300 400 500 600 Time in samples Figure 3.69: Estimated time delay and gain of the prebleaching process by MMAC algorithm. 80 Chapter 3: MMAC as a dead-time compensator Kappa number Normal gain 10% increase in gain ic wiWA - 100 _r- I I 200 300 A- - - -w Jr- I 400 - I 500 600 Time in samples 700 800 900 1000 C12-i-C102 (g/L) i0 I I 80 00 200 300 400 500 600 Time in samples 700 800 900 1000 Figure 3.70: Prebleaching kappa number (process output) and bleaching agents flow (control signal). No adaptation. Kappa number Norrnaigain 10% increase in gain ic 00 -- — 200 300 ---j 400 V • 500 600 Time in samples 700 800 900 1000 C12+C102 (gIL) :OC I I 80 200 300 400 500 600 Time in samples 700 800 900 1000 Figure 3.71: Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by FMVRE/RLS. 81 Chapter 3: MMAC as a dead-time compensator Kappa number -• Normalgain 10% Increaseln gain lc C roo — — I 200 — --- I 400 300 500 600 Time in samples _.j.__ - — I I 700 800 900 lOC 700 800 900 1000 C12--C102 (g/L) oc p 80 60 00 200 400 300 500 600 Time in samples Figure 3.72: Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by MMAC. Estimated time delay I I I 6 4 2 100 200 300 400 500 600 Time in samples 700 800 900 1000 700 800 900 1000 Estimated gain foo Figure 3.73: Estimated 200 400 300 time delay and 500 600 Time in samples gain of the prebleaching process by FMVRE/RLS algorithm. 82 MMAC as a dead-time compensator Chapter 3: Estimated time delay I. — 100 200 1 300 I 400 I — 500 600 Time in samples 700 800 l0( 900 Estimated gain -0.02 :j4v -0.08 200 300 400 500 600 Time in samples 700 900 800 1000 Figure 3.74: Estimated time delay and gain of the prebleaching process by MMAC algorithm. Kappa number 15 I :: I I I 700 800 900 1000 700 800 900 1000 10 l00 200 300 400 500 600 Tune in samples C12+C102 (g/L) 100 I I 80 60 40 20 00 200 300 400 500 600 Time in samples Figure 3.75: Prebleaching kappa number (process output) and bleaching agents flow (control signal). No adaptation. 83 Chapter 3: MMAC as a dead-time compensator Kappa number d=5 d=2 K roo —-- 200 300 400 - 500 600 Time in samples - - ---- - 700 800 900 1000 700 800 900 1000 C12i-C102 (g/L) 80 60 00 200 300 400 500 600 Time in samples Figure 3.76: Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by FMVRE/RLS. Kappa number 15 I I I d=5 ic 00 d=2 - 200 300 I I 200 300 400 y_ - 500 600 Time in samples ..- --- 700 800 900 10( 700 800 900 1000 C12+C102 (g/L) bc 80 60 00 400 500 600 Time in samples Figure 3.77: Prebleaching kappa number (process output) and bleaching agents flow (control signal). Adaptation by MIvIAC. 84 Chapter 3: MMAC as a dead-time compensator Estimated time delay 100 200 300 400 500 600 Time in samples 700 800 900 1000 Estimated gain 0 I I I -0.02 EE 200 300 400 500 600 Time In samples 700 800 900 1000 Figure 3.78: Estimated time delay and gain of the prebleaching process by FMVRE/RLS algorithm. Estimated time delay 100 200 300 400 500 600 Time in samples 700 800 900 1000 Estimated gain C I I -0.02 -0.08 100200 3004005006007008009001000 Time in samples Figure 3.79: Estimated time delay and gain of the prebleaching process by MMAC algorithm. 85 Chapter 3: MMAC as a dead-time compensator Kappa number (Adaptation by FMVREIRLS) 00 200 300 400 500 600 Time in samples 700 800 900 1000 800 900 1000 Kappa number (Adaptation by MMAC) 11 W 9 00 200 300 400 500 600 Time in samples 700 Figure 3.80: Prebleaching kappa number (process output) at a closed-loop pole of Q=-O.77. C12-i-C102 (g/L) (Adaptation by FMVRE/RLS) DO 200 300 400 500 600 Time in samples 700 800 900 1O( 800 900 1000 C12-i-Q02 (g/L) (Adaptation by MMAC Ioo 200 300 400 500 600 Time in samples 700 Figure 3.81: Bleaching agents flow (control signal) at a closed-loop pole of Q=-O.77. 86 Chapter 3: MMAC as a dead-time compensator Estimated time delay by FMVREIRLS 2 .11 0 200 300 400 500 600 Time in samples 700 800 900 1000 Estimated time delay by MMAC 6 4 ñJ 4 L 1 2 I 10 200 300 400 500 600 Time in samples 700 800 900 1000 Figure 3.82: Estimated process time delay at a closed-loop pole of Q-O.77. Estimated process gain by FMVREIRLS 0 200 300 400 500 600 Time in samples 700 800 900 10 800 900 1000 Estimated process gain by MMAC -0.02 -0.04 -0.06 -0.08 200 300 400 500 600 Time in samples 700 Figure 3.83: Estimated process gain at a closed-loop pole of Q=-O.77. 87 Chapter 3: MMAC as a dead-time compensator Kappa number (Adaptation by FMVRE/RLS) 2003004005006007008O0900l000 Time in samples Kappa number (Adaptation by MMAC) 12 I 100 200 300 400 500 600 Time in samples I I I 700 800 900 1000 Figure 3.84: Prebleaching kappa number (process output) at a closed-loop pole of Q=-O.4. C12+C102 (g/L) (Adaptation by FMVRE/RLS) 40 II 30 20 10 0 200 300 400 600 500 Time in samples 700 800 900 1000 800 900 1000 C12+Q02 (gIL) (Adaptation by MMAC) 40 30 20 10 00 200 300 400 600 500 Time in samples 700 Figure 3.85: Bleaching agents flow (control signal) at a closed-loop pole of Q=-O.4. 88 Chapter 3: MMAC as a dead-time compensator Estimated time delay by FMVRE/RLS 8 I I I 800 900 I I 800 900 6 4 100 200 300 I I I 200 300 400 400 500 600 Time In samples 700 1000 Estimated time delay by MMAC I I 6 100 500 600 Time in samples 700 1000 Figure 3.86: Estimated process time delay at a closed-loop pole of Q=-O.4. Estimated process gain by FMVRE/RLS C -0.02 -0.04 -0.06 -0.08 0 200 300 400 500 600 700 800 900 1( 00 800 900 1000 Estimated process gain by MMAC I.. -0.02 -0.04 -0.06 -0.08 ioo 20 300 40 500 600 Time in samples 700 Figure 3.87: Estimated process gain at a closed-loop pole of Q=-O.4. 89 Chapter 3: MMAC as a dead-time compensator 3.5 Conclusions Two different algorithms for on-line gain estimation and dead-time compensation in adaptive control were introduced. The performances of the algorithms were first studied within the frame of servo problems utilizing generic models. Later the same study was done in a regulation problem in order to control a chlorination tower in a prebleaching process in closed-loop. Assuming that the unknown delay is known to lie within a discretized delay interval of uncertainty, the FMVRE/RLS algorithm compensates for the gain and time delay variations by initially estimating the time delay directly. This is achieved by computing the cross-correlation function at every sample between the control signal and the incremental changes of the process output for every delay in the range of interest. Time delay is then estimated by finding the time delay that maximizes (for positive gain) or minimizes (for negative gain) the cross-corelation function in that range. Once an estimate of the time delay is available the information is fed into a RLS based estimator so that a regression and a parameter vector can be constructed based on the estimated delay for the purpose of gain estimation. MMAC performs the same task by first assuming a prior knowledge of the time delay and gain interval of uncertainty. Then upon discretization of these intervals a set of candidate process models are constructed based on the elements of the discretized intervals. A model for the process is formed by calculating the weighted sum of the candidate process models where each weight reflects the probability of its corresponding candidate model to closely represent the true gain and time delay, and is calculated based on the magnitude of the prediction error computed for each model. Both FMVREIRLS and MMAC provide satisfactory dead-time compensation. However, in servo problems where the excitation in the system is rich FMVREIRLS is more versatile. Not only is it capable of estimating the process parameters accurately, but also it can estimate the parameters of a higher-order model, and therefore achieve good closed-loop performance. MMAC on the other hand always limits itself to assuming a first-order model of the true process. Furthermore, the true gain can not be estimated accurately since in the steady state conditions the estimates tend to converge to some median in the gain uncertainty interval. As result the controller’s performance can be sluggish 90 Chapter 3: MMAC as a dead-time compensator or introduce some oscillations to the output which can be a satisfactory perfonnance but not to the expectations. In regulating the chlorination tower around a fixed set-point both algorithms provide good dead-time compensation. However, due to low excitation conditions the parameter estimator in the FMVRE./RLS algorithm looses its activeness, and as a result, like the MMAC algorithm , it has to assume a first-order process model when in fact the true process has a higher order. This is achieved through lowpass filtering of the required signals for identification which will force the closedloop bandwidth to be equal or less than the filter’s cutoff frequency. Under such conditions both algorithms provide satisfactory control. However, under tight closed-loop conditions the MMAC algorithm shows greater robustness. The FMVREIRLS achieves parameter estimation by performing a small number of computations at every sample. The MMAC algorithm however, involves a large number of computations at every sample. This is a result of performing the same number of computations for every candidate process model. 91 Chapter 4: General Conclusions and future work Chapter 4 General Conclusions and future work 4.1 Thesis summary In the first half of this work a method based on the theory of Multi-Model Adaptive Control for the purpose of wood species identification and compensation in a kraft pulping process was presented. The ability of the algorithm to successfully identify the ratios of different wood species in the incoming wood chips feed stream was demonstrated. It was shown that an adaptive controller based on a multimodel approach can utilize the identified wood species ratios in order to take compensating actions for wood species variations. The control provided by the proposed algorithm showed better performance over a conventional adaptive controller from the output variance standpoint. In the remaining half of this thesis a different application of adaptive control based on a multi-model approach was investigated. The algorithm proposed by Gendron provides an effective way for time delay compensation while tracking a variable process gain. The performance of this algorithm was compared against the Fixed Model Variable Regression Estimator, FMVRE, proposed by Elnaggar, coupled with a recursive least-squares based estimator, RLS, to provide the same estimates of the process time delay and gain respectively. Both algorithms provide satisfactory adaptation of time delay under servo and regulation conditions. However, in servo problems where the change in the reference signal introduces strong excitations into the system, the RLS based estimator is capable of estimating the parameters of a higher order model which can lead to smaller model uncertainties and therefore a better control. MMAC on the other hand showed a greater robustness under regulation and tight closed-loop control conditions. 4.2 Future research One issue that has not been addressed throughout this work is the stability of the proposed multi model controller. The calculated control action is a weighted sum of the control signals calculated for each individual model. Though the overall control action resulting from a weighted sum might be 92 Chapter 4: General Conclusions and future work stable, each individual control signal could be growing infinitely, leading to numerical instabilities. Therefore a complete analysis of the stability of the multi-model controller is desirable, and this is left for future work. The theory of MMAC can be applied to many pulping processes that accommodate wood chips of different textures or pulps of different grades. Among these processes are mechanical pulping and bleaching where it is desirable to identify and compensate for the species variations. However, without having accurate mathematical models describing the kinetics or dynamics of individual wood species, the implementation of MMAC for the purpose of automatic wood species identification is not possible. Unfortunately, up until now, the only area where extensive research has been done on finding kinetic models of different wood species is kraft pulping. In the area of bleaching a lot of research has gone towards improving a kinetic model that describes bleaching of only one species rather than finding kinetic models for different species. In mechanical pulping it has been shown that the relation between the pulp freeness and the specific energy put into the system resembles an exponential decay. However, mathematical models showing this relation for different species have not yet been derived. 93 References [1] Ackert, J., L. Edwards and H. Norrstrom (1975). ‘How to use less chlorine in kraft pulp bleaching’. Pulp and Paper Canada 76(2), 96—98. [2] Aström, K. J. and B. Wittenmark (1989). Adaptive Control. Addison-Wesley Publishing Company. [3] Bryce, J. R. G. (1980). Alkaline pulping. In J. P. Casey (Ed.). ‘Pulp and Paper Chemistry and Chemical Technology’, third edn. Vol. 1. John Wiley and Sons. [4] Dahlin, E. B. (1968). ‘Designing and tuning digital controllers’. Instruments and Control Systems 41, 77—83. [5] Dumont, G. A. (1983). Analysis of the design and sensitivity of the Dahlin regulator. In ‘69th Anual Meeting, Technical Section, Canadian Pulp and Paper Association’. Montréal. pp. 37—50. [6] Dumont, G. A. (1989). System identification and adaptive control in papermaking. In C. F. Baker (Ed.). ‘Fundamentals of Papermaking’. Vol. 2. Mechanical Engineering Publication Limited. London. pp. 1151—1181. [7] Elnaggar, A. (1990). Variable Regression Estimation of Unknown System Delay. PhD thesis. Department of Electrical Engineering, University of British Columbia. Vancouver, B.C., Canada. [8] Elnaggar, A., G. A. Dumont and A. L. Elshafei (1992). Adaptive control with direct delay estimation. In ‘Control System’s 92, Dream vs. Reality:Modem Process Control in the Pulp and Papert Industry’. Whistler, B.C.. pp. 13—17. [9] Gendron, S. and A. P. Holko (1991). Simple adaptive digital dead-trtime compensators for low order siso processes. In ‘Preprints of the Ninth IFAC Symposium on System Identification and Parameter Estimation’. Budapest, Hungary. pp. 1135—1140. 94 [10] Gendron, S., 3. Barrette, M. Perrier and N. Legault (1992). Industrial experience with model weighting adaptive control. In ‘Control Systems 92, Dream vs. Reality:Modem Process Control in the Pulp and Papert Industry’. Whistler, B.C.. pp. 171—179. [11] Germgârd, U. (1982). ‘Stoichiometry of prebleaching of softwood kraft pulp with chlorine dioxide and small fractions of chlorine’. Cellulose Chemistiy and Technology 48(17), 35—48. [12] Germgkd, U., A. Teder and D. Tormund (1982). ‘The relative rates of consumption of chlorine and chlorine dioxide during (d+c) bleaching of softwood kraft pulp’. Tappi 65(5), 124—126. [13] Germgàrd, U. and Hans Lindberg (1982). ‘A mathematical model for (D+C)-prebleaching’. Svensk Papperstidning 85(18), R172. [14] Goodwin, G. C. and R. L. Payne (1977). Dynamic System Identfication:Experimental Design and data analysis. Academic Press. New York. [15] Gustafson, R. R., C. A. Sleicher, W. T. McKean and B. A. Finlayson (1983). ‘Theoretical model of the kraft pulping process’. Industrial and Engineering Chemistiy Process Design and Development 22(1), 87—96. [16] Hatton, 3. V. (1973a). ‘Application of empirical equations to kraft process control’. Tappi 56(8), 108—111. [17] Hatton, J. V. (1973b). ‘Development of yield prediction equations in kraft pulping’. Tappi 56(7), 97—100. [18] Hatton, J. V. (1976). ‘The potential of process control in kraft pulping of hardwoods relative to softwoods’. Tappi 59(8), 48—50. [19] Hatton, J. V. and J. L. Keays (1971). Effect of raw material and process variables on close control of the kraft pulping process. Technical report. Department of the Environment, Canadian Forestry Service, Western Forest Products Laboratory. [20] Kerr, A. J. (1970). ‘The kinetics of kraft pulping:progress in the development of a mathematical model’. Appita 24(3), 180—188. [21] Kocurek, M. J. (Ed.) (1989). Pulp and Paper Manufacture The Joint Textbook Committee of the Paper Industry. 95 Alkaline Pulping. Vol. 5. 3 edn. [22] Luyben, W. L. (1990). Process Modeling, Simulation, and Control for Chemical Engineers. second edn. McGraw Hill Publishing Company. [23] Salgado, M. E., G. C. Goodwin and R. H. Middleton (1988). ‘Modified least squares algorithm incorporating exponential resetting and forgetting’. International Journal of Control 47(2), 477— 491. [24] Shah, S. L. and W. R. Cluett (1991). ‘Recursive least squares based estimation schemes for self-tuning control’. The Canadian Journal of Chemical Engineering 69, 89—96. [25] Singh, R. P. (1979). The Bleaching of Pulp. 3 edn. TAPPI Press. [26] Smith, C. C. and T. J. Williams (1975). A kinetic model of kraft pulping reactions. In T. J. Williams (Ed.). ‘Modeling and Control of Kraft Production Systems’. Proc. ISA Symp.. Milwaukee. pp. 7—20. [27] Smith, C. L. (1972). Digital Computer Process Control. Intext Educational Publishers. [28] Smook, G. A. (1982). Handbook for pulp and paper technologists. The Joint Textbook Committee of the Paper Industry. [29] Vroom, K. E. (1957). ‘The H factor: The means of expressing cooking times and temperatures as a single variable’. Pulp and Paper Magazine of Canada 58:3, 228—23 1. 96
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of multi-model adaptive control to pulping...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of multi-model adaptive control to pulping processes Mohsenian, Navid 1994
pdf
Page Metadata
Item Metadata
Title | Application of multi-model adaptive control to pulping processes |
Creator |
Mohsenian, Navid |
Date Issued | 1994 |
Description | In this thesis, the theory of Multi-Model Adaptive Control is applied to pulping processes, particularly for wood species compensation and identification, and dead-time compensation. Wood species identification is performed on a kraft pulping process whose kinetics are described by the Hatton’ s equations where the ratios of the different wood species in the incoming feed stream are identified on-line by a recursive least-squares based estimation algorithm. Upon wood species identification, a multi-model adaptive controller takes action to compensate for the variability in the wood chips supply content and to control the kappa number. Dead-time compensation is applied to a chlorination tower in a bleaching process with chlorine and chlorine dioxide. The corresponding mass balance equations are developed. The process is simulated, sampled, and identified in open-loop. A multi-model adaptive controller is utilized in closed-loop to compensate for the variations in the pulp quality and transport time delay and to control the bleaching kappa number. The controller’s performance is compared against an adaptive controller which is coupled with a recursive least-squares based algorithm and the Fixed Model Variable Regression Estimation algorithm to compensate for the same variations. |
Extent | 2835829 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-02-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065058 |
URI | http://hdl.handle.net/2429/4935 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1994-0106.pdf [ 2.7MB ]
- Metadata
- JSON: 831-1.0065058.json
- JSON-LD: 831-1.0065058-ld.json
- RDF/XML (Pretty): 831-1.0065058-rdf.xml
- RDF/JSON: 831-1.0065058-rdf.json
- Turtle: 831-1.0065058-turtle.txt
- N-Triples: 831-1.0065058-rdf-ntriples.txt
- Original Record: 831-1.0065058-source.json
- Full Text
- 831-1.0065058-fulltext.txt
- Citation
- 831-1.0065058.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065058/manifest