UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

On dual control and adaptive Kalman filtering with applications in the pulp and paper industry Ismail, Ahmed Abdel-Rahman 2001

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2001-714772.pdf [ 5.74MB ]
Metadata
JSON: 831-1.0065524.json
JSON-LD: 831-1.0065524-ld.json
RDF/XML (Pretty): 831-1.0065524-rdf.xml
RDF/JSON: 831-1.0065524-rdf.json
Turtle: 831-1.0065524-turtle.txt
N-Triples: 831-1.0065524-rdf-ntriples.txt
Original Record: 831-1.0065524-source.json
Full Text
831-1.0065524-fulltext.txt
Citation
831-1.0065524.ris

Full Text

O N D U A L C O N T R O L A N D A D A P T I V E K A L M A N F I L T E R I N G W I T H A P P L I C A T I O N S I N T H E P U L P A N D P A P E R I N D U S T R Y B y Ahmed Abdel-Rahman Ismail B.Sc. Cairo University, 1993 M.Sc . Cairo University, 1996 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F E L E C T R I C A L A N D C O M P U T E R E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 2001 © Ahmed Abdel-Rahman Ismail, 2001 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of Br i t i sh 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. Department of Electrical and Computer Engineering The University of Br i t i sh Columbia 2356 M a i n M a l l Vancouver, Canada V 6 T 1Z4 Date: Abstract This thesis is about dual control of time-varying stochastic processes in the pulp and paper industry with parameter estimation carried out via adaptive K a l m a n filtering. It can be divided into three parts. The first part deals with the control of a chip refiner to produce wood pulp, where the process gain between the refiner motor load and the plate gap is both nonlinear and time-varying, with reversal in the sign of the gain indicating the onset of pulp pad collapse. The control objective is to regulate the motor load while avoiding pad collapse. The problem is principally stochastic in nature, since the gap at which gain reversal occurs can wander unpredictably. The proposed control strategy consists of an active subop-t imal dual controller coupled with an adaptive Ka lman filter for parameter estimation. The controller minimizes a nonlinear performance index designed especially to reflect the peculiarities of the process. Thus, no heuristic logic is needed. Simulations show the superior performance offered by this strategy. The second part is concerned with C D coat weight control on bent blade coaters. The coater is a coupled multivariable process whose gain drifts over time and often switches sign. Current industrial practice is to switch off automatic control when the loop becomes unstable due to gain sign reversal. Because of this, the standard industrial controller is rarely on for more than half of the blade life. The proposed control strategy relies on the general framework introduced for the chip refiner problem. The controller takes into consideration the loading and bending limitations imposed by the actuators and the blade. A n approximate analytical control law was derived to allow for easy implementation. The proposed strategy was successfully applied to an off-machine industrial coater. The results of a series of trials show the advantage of including probing in the control signal, and that the adaptive Ka lman filter was capable of tracking gain variations. The dual controller yielded substantial quality improvement and was able to control the process throughout the entire blade life. The developed strategy was well accepted by the company. The adaptive Ka lman filter used for both applications had the disadvantage that i i one has to wait for some time before starting to use the noise variance estimates for parameter estimation. In the third part, a novel adaptive K a l m a n filtering scheme is derived to resolve this problem. This scheme is based on the Expectat ion-Maximizat ion ( E M ) algorithm. It uses a number of fixed-point smoothers running in parallel to on-line estimate the variances of the process and measurement noise for a general linear, discrete, time-varying stochastic system. This approach implies that the estimates are used in the K a l m a n filter. The approach is evaluated through a number of simulation experiments. i i i Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement xi 1 Introduction 1 1.1 The Dual Control Concept 1 1.2 Problem Formulation 3 1.3 Optimal Dual Control 5 1.4 Suboptimal Dual Control 9 1.5 When to Use Dual Control? 10 1.6 The Thesis 10 1.6.1 Objectives 10 1.6.2 Contributions 11 1.6.3 Overview 11 2 Dual Control of Chip Refining 13 2.1 Introduction 13 2.2 Process Description 16 2.3 Process Modeling 17 2.4 Parameter Estimation 19 iv 2.5 Controller Design 23 2.6 Simulations 26 2.7 Conclusions 30 3 Dual Control of Paper Coating 35 3.1 Introduction 35 3.2 Blade Coating 37 3.2.1 Coat Weight Blade Pressure Relationship 38 3.2.2 Actuators and Measurements 40 3.3 Model ing of the C D Control Process 42 3.4 Parameter Estimation 45 3.5 Controller Design 49 3.5.1 Accommodating Constraints 52 3.6 M i l l Trials 54 3.6.1 Implementation Aspects 55 3.6.2 Experimental Results and Discussion 58 3.7 Conclusions 61 4 A New Approach to Adaptive Kalman Filtering 92 4.1 Introduction 92 4.2 Problem Statement 93 4.3 The E M Algor i thm 94 4.4 E M Estimation of Q and R 100 4.5 A Recursive E M Update 103 4.6 Simulations 107 4.7 Conclusions . . . 119 v 5 Summary and Future Work Bibliography A Approximation of the A S O D Controller B Proofs for Chapter 3 C Estimation of the Noise Filter Zero D Proofs for Chapter 4 v i List of Tables 3.1 Ga in estimation (data courtesy of Honeywell-Measurex) vn List of Figures 1.1 Self-tuning adaptive control system 2 1.2 Block diagram of the optimal dual controller 5 2.1 Schematic diagram of a chip refiner 14 2.2 Motor load versus plate gap 17 2.3 The performance index J(Au(t)) (load below setpoint) 27 2.4 The performance index J(Au(t)) (load above setpoint) 28 2.5 Safe operation, then unreachable setpoint 31 2.6 More of normal operation 32 2.7 Operation around zero gain 33 2.8 Behavior of the controller when the pad collapses 34 3.1 Blade coater [63] 37 3.2 Schematic diagram of a bent blade coater [37] 38 3.3 C D coat weight as a function of actuator position 39 3.4 Bent blade vs. stiff (bevelled) blade 40 3.5 C D coat weight actuators [58] 41 3.6 Diagonal scan path introduces M D variations into (what should be flat) measured profile 42 3.7 Flow of information in the C D control loop 42 3.8 Block diagram of the C D control process 43 3.9 Coating station where trials were performed 55 3.10 Identification of the actuator spatial response 57 v i i i 3.11 Tr ia l 1, the dual controller provides about 50% reduction in C D coat weight variations 63 3.12 3-D plots of tr ial 1 64 3.13 3-D plots of t r ia l 1 (continued) 65 3.14 Tr ia l 2, the performance of the dual controller is comparable to the Dahl in controller 66 3.15 3-D plots of tr ial 2 67 3.16 3-D plots of tr ial 2 (continued). 68 3.17 Noise variance estimation for tr ial 2 69 3.18 Tr ia l 3, the benefit of probing is illustrated 70 3.19 3-D plots of tr ial 3 71 3.20 3-D plots of tr ial 3 (continued) 72 3.21 Noise variance estimation for tr ial 3 73 3.22 Tr ia l 4, performance degrades as the coater drifts into stiff blade operation. 74 3.23 3-D plots of t r ia l 4 75 3.24 3-D plots of tr ial 4 (continued) 76 3.25 Noise variance estimation for tr ial 4 77 3.26 Tr ia l 5, the dual controller reduces coat weight variations in stiff blade mode. 77 3.27 3-D plots of tr ial 5 78 3.28 3-D plots of tr ial 5 (continued) 79 3.29 Tr ia l 6, the dual controller recovers from a major disruption 80 3.30 3-D plots of tr ial 6 81 3.31 3-D plots of tr ial 6 (continued) 82 3.32 Tr ia l 7, the dual controller outperforms the Dahl in controller 83 3.33 3-D plots of tr ial 7 84 ix 3.34 3-D plots of t r ia l 7 (continued) 85 3.35 Tr ia l 8, the dual controller is running for about 4 hours 86 3.36 3-D plots of tr ial 8, scan number 1-600 87 3.37 3-D plots of t r ia l 8, scan number 1-600 (continued) 88 3.38 Noise variance estimation for tr ial 8, scan number 1-600 89 3.39 Noise variance estimation for tr ial 8, scan number 600-1200 89 3.40 3-D plots of tr ial 8, scan number 600-1200 90 3.41 3-D plots of tr ial 8, scan number 600-1200 (continued) 91 4.1 Batch estimation of q and r for example 1 109 4.2 On-line estimation of q and r (solid) and Isaksson's method estimates (dashed line) for example 1 110 4.3 Batch estimation of Q and r for example 2 I l l 4.4 Batch estimation of Q and r for example 2 (continued) 112 4.5 On-line estimation of Q and r using the straight averaging scheme ((4.83) and (4.84)) for example 2 113 4.6 On-line estimation of Q and r using the straight averaging scheme ((4.83) and (4.84)) for example 2 (continued) 114 4.7 On-line estimation of Q and r (solid) and Isaksson's method estimates (dashed line) for example 2 115 4.8 On-line estimation of Q and r (solid) and Isaksson's method estimates (dashed line) for example 2 (continued) 116 4.9 On-line estimation of Q and r for a non-diagonal Q 117 4.10 On-line estimation of Q and r for a non-diagonal Q (continued) 118 A . l The loss function J together with its quadratic and nonlinear parts 129 x Acknowledgement I would like to express my sincere gratitude to my supervisor Dr . Guy Dumont. Without his encouragement, support and technical guidance this thesis would never have been completed. I would also like to thank Dr . Michael Davies and other past and present members of the Process Control group at the Pu lp and Paper Centre for creating such a stimulating research environment. I am also very grateful to Dr . A l f Isaksson whose enthusiasm for the subject of Adaptive Fi l ter ing was catching me. Alf , thank-you for providing the seed of Chapter 4, putting your trust in me and most importantly your friendship. The industrial work described in this thesis was carried out in collaboration with Honeywell-Measurex. Special thanks to Johan Backstrom for his technical help and valu-able discussions/arguments. Further I wish to thank my parents for their understanding and support, but most of al l for their unwavering love. Finally, I owe a great deal to my friends Ayman , Faten, Tamer, A m r , Hai tham and Eman who have been more than a family to me. M y only regret is that this means my days in Vancouver with them are over. x i Chapter 1 Introduction 1.1 The Dual Control Concept For a large number of control problems, process parameters are either unknown or partly known. One way to handle process uncertainty is to use an adaptive controller [7]. Fig-ure 1.1 shows the general structure of a self-tuning adaptive control system. The estimator uses the inputs and outputs of the process to identify the unknown process parameters. The design block then uses this information to implement the control algorithm. This type of adaptive controllers consists of an ordinary feedback loop and a controller updat-ing loop. Different algorithms are obtained depending on the design of the estimator and the controller. Classical adaptive control theory is almost entirely based on a heuristic approach. The unknown system parameters are estimated in real time, and the estimated parameters are then used as if they were the true ones, i.e. estimation uncertainties are not taken into account in the design. This is called the certainty equivalence principle. One case in which the certainty equivalence principle holds is the celebrated linear quadratic Gaussian case for known systems [7]. If, however, there are unknown parameters in the process to be controlled, then certainty equivalence does not hold. In this case, it is necessary to simultaneously control the system and estimate the unknown parameters. In order to properly estimate the unknown parameters, the control signal must disturb the system. This conflicts with the objective for the output signal to track the desired reference trajectory with as small deviations as possible. The compromise between these two goals -estimation and control- leads to a controller with dual features. Feldbaum presented the concept of dual control in a landmark series of papers in 1960-61 [28]. In those, he showed that to optimally control an unknown stochastic system with 1 Introduction 2 Specification Process parameters Design Estimator Reference signal Controller Control Process Output Figure 1.1: Self-tuning adaptive control system. time-varying parameters, one has to solve a functional equation, the Bellman equation, which can not be solved analytically. Further, it can be solved numerically only in very few simple cases, thereby making the optimal controller intractable for most cases. Yet, the controller exhibits a remarkable property. It attempts to drive the output to the desired value, but w i l l also introduce perturbations when the estimates are uncertain. This wi l l improve the accuracy of future estimates, consequently, leading to better control. The controller thus achieves the finest balance between good control and small estimation errors. This is why it is called optimal dual control. The computational complexity associated with the optimal dual controller has mo-tivated the development of suboptimal dual controllers that attempt to retain the dual characteristic of the optimal controller while being practically implementable. The dual property is of great importance since it introduces probing or active learning into the controller. This is in contrast to non-dual controllers where the learning is accidental or passive, i.e. the controller does not take any active measures to improve the estimates. The certainty equivalence controllers are typical of this class. A n analogy can be drawn between the operation of a dual controller and driving a car. The driver wants to drive a car at a steady speed to a specific destination. But in order to direct his/her own actions better, the driver must understand the car better. Therefore, sometimes he/she may slightly turn the wheel or nudge the pedals not in order to obtain direct advantage, but mainly with the aim of finding out how the car reacts. Introduction 3 Thus, driving a car and understanding how it reacts are closely linked to each other. 1.2 Problem Formulation The basic stochastic model used in stochastic adaptive control is a time-varying ARX-model y(t) + ai{t)y{t- 1) + • • • + an(t)y{t-n) = b0(t)u(t- 1) + • • • + 6 n _i (*)«(*-n) + e(t) (1.1) where y, u, and e are the process output, control signal and disturbance signal, respec-tively. The noise sequence e is assumed Gaussian white noise with zero mean and vari-ance R. The time-varying parameters 6(t) = are modeled by a Gauss-Markov process b0(t) ••• bn-i(t) ai(t) ••• an(t) 9(t + I) = $9(t) + w(t) (1.2) where $ is a known constant matrix and w is another sequence of Gaussian white noise with zero mean and variance Q. The special case of constant unknown parameters can be studied by setting $ = I and Q = 0. The initial state of the system in (1.1) is assumed to be normally distributed with mean value E {#(0)} = 90 and variance Var (0(0)) = P0. Further, it is assumed that w is independent of e and 0(0). The system in (1.1) can be modeled in a time-varying state-space form 9{t+l) = $9(t)+w(t) y(t) = ^(t)9(t) + e(t) (1.3) where <p(t) u(t — 1) • • • u(t — n) —y(t — 1) • • • —y(t — n) (1.4) In order to adaptively control the system (1.3) one has to estimate the time-varying parameters. Most adaptive controllers are based on the separation principle. This implies that the unknown parameters are estimated separately from the design part. The control Introduction 4 design part then uses the result from the estimator to determine the control signal, with or without any intention to introduce dual features. It is not hard to show that, under the Gaussian assumption about the disturbances, the conditional distribution of 9(t + 1) given all values of the output up to and including y(t) is Gaussian with mean 9(t + 1) and variance P(t+ 1), where 9(t + 1) and P(t + 1) satisfy the Kalman filter equations [6] 9(t + 1) = $6{t) + K(t) (y(t) - <pT(t)6{tj) (1.5) P(t + 1) = ($ - K(t)cpr(t)) P(t)$T + Q (1.6) K{t) = $P(t)<p(t) (R + ^WPtfipit))'1 (1.7) We then assume that the loss function to be minimized is JN = E jf; (y(t + i)- yr(t + z))2 j (1.8) where yr is the reference signal and E denotes mathematical expectation. This is called the TV-stage criterion. The loss function should be minimized with respect to the control signal u(t),u(t + l),... ,u(t + N — 1). At time t, old inputs and outputs are collected into Yt-i = {y(t — 1), u(t — 1),.. . , 2/(0), u(0)}. If the parameters were known, it is assumed that the optimal controller is u{t) = fknomt{y(t),Yt_u6(t)) (1.9) The adaptive controller based on the certainty equivalence principle is UCE(<) = /known t 1 ' 1 0 ) i.e. the strategy obtained in the case of known parameters is used and the true parameters are substituted by their estimates. Parameter uncertainties are not taken into account in the design. This approach is widely used in practice. When the loss function (1.8) is minimized one step ahead (i.e. N — I) based on the Introduction 5 Specification Reference signal Hyperstate Nonlinear control law Control Calculation of hyperstate Process Outni ut Figure 1.2: Block diagram of the optimal dual controller. current estimates, the resulting controller takes the parameter uncertainties into consider-ation. This type of controller is called a cautious controller [7], because its gain decreases as uncertainty increases. However, with small control signals less information will be gained about the process and parameter uncertainties will increase and even smaller con-trol signals will be generated. This vicious circle may lead to turn off of the control. The problem with turn-off makes the cautious controller unsuitable for practical applications. The cautious controller has the structure "cautious(*) = / c a u t i o u s Y t - U P(t)) ( I - 1 1 ) 1.3 Optimal Dual Control The difficulties of computing the optimal dual controller will now be illustrated. Remem-ber that we have the loss function (1.8) from the previous section. The controller obtained for iV = 1 is a cautious controller and is sometimes called a myopic controller since it is short-sighted and looks only one step ahead. The optimal controller will be very different when ./V is large. Define V(£(t + l),t + l ) = min E{ £ (y(k) - yr(k)f \Yt 1 (1.12) U (t) . . .u( iV- l) \ k = t + 1 I Introduction 6 V(£(t + 1), t + 1) is the minimum expected loss for the remaining part of the control horizon given data up to t and cf is the so-called hyperstate of the problem, which includes the parameter estimates, their uncertainties, old inputs and old outputs of the system. Using dynamic programming it can be shown that the optimal dual controller satisfies the Bel lman equation [7] V(£(t + ! ) , £ + !) = m i n E {(y(t + 1) - yr(t + l ) ) 2 + V(£.(t + 2), t + 2) \Yt} (1.13) The Bel lman equation contains nested minimization and mathematical expectation, and even in the simplest cases it has to be solved numerically. Since both V and u have to be discretized it follows that the storage requirements increase drastically with increasing grid size and the solution of the Bellman equation requires excessive computations. The optimal controller is a nonlinear function mapping the hyperstate into the control signal, see figure 1.2. Notice the similarity with the self-tuning regulator in figure 1.1. The structure simplicity of the controller is obtained thanks to the introduction of the hyperstate. To illustrate the different design methods, consider the following example. Consider an integrator in which the unknown gain is changing, i.e. we have the model where e is white noise N(0, v). Assuming zero a reference signal, the one-stage loss function is Example y(t)-y(t-l) = b(t)u(t-l) + e(t) (1.14) (1.15) Were b known, the solution is t r iv ia l u(t) 1 y(t) (1.16) b(t + l) Introduction 7 The minimum variance certainty equivalence controller is then and, b(t + 1) J l 0 E = i / 2 + ^ f t + 1 ) y 2 ( t ) (1.18) C E fe2(i + l) where P 6 is the variance of the estimate b. Noticing that the conditional distribution of y(t + 1) given Yt is Gaussian, J\ is eval-uated Jx = (y(t) + b{t + l)u(t)f + Pb(t + l)u2{t) + v2 (1.19) The cautious controller that actually minimizes J\ is McautiousW = b(t + 1) (1.20) W 62(t + l) + P6(t + l) W V ' Notice that the gain of the cautious controller approaches zero as the uncertainty increases. Notice also that the cautious controller reduces to the certainty equivalence one (1.17) when Pb = 0. The minimum loss function is given by •A t- =v2 + -~ n ( t + 1 ) y2(t) (1.21) 1 c a u t i o u s 6 2 ( t + 1 ) + p 6 ( t + 1 ) y W Let us now consider the two-stage optimization problem. The optimal dual controller that minimizes the two-stage loss function satisfies the Bellman equation V = mmE{y2(t + l)+ min E {y2(t + 2) \Yt+1} \Yt) (1.22) u(t) y u(t+l) ( ' J Assuming that an optimal value for u(t + 1) will be used, we get V = min E L2(t + 1) + v2 + , P ^ + 2 ) 2 ( t + 1 } ] y ] ( 1 2 3 ) Define Introduction 8 where the hyperstate £(£ + 1) is b(t + 1), Pb[t + Appl icat ion of the K a l m a n filter (1.5)-(1.7) to the process model of (1.14) is straightforward. It is then easy to find T(y(t + l ) ,u ( tU(< + l)) = [Kt + 1) + USfflW)) + !) " V® " + !))]2 + ^ iSfc) + ^  (1.25) where p2 is the variance of the white noise w in b(t + 1) = b(t) + w(t). Then, V = min E {y\t + 1) + v2 + T (y(t + 1), u ( t ) , f (* + 1)) y 2 ( t + 1) \Yt} (1.26) u(t) *• -1 Even though the distribution of y(t + 1) given Yt is Gaussian, the optimization problem of (1.26) requires extensive computations as the expectation has to be evaluated numeri-cally function of a discretized u(t). The optimal dual controller, although known in principle, is of lit t le practical use, as it is impossible to implement for realistic processes. Yet, it gives insight into the character and the advantages of the dual property. There are some simple examples where the optimal dual controller has been obtained numerically, see, for instance [4, 5]. Simulations show that the optimal controller, the cautious controller and the certainty equivalence controller are the same when the estimates are accurate. The optimal control law is also close to the cautious control for large control errors. For estimates with poor precision and moderate control errors the optimal control gives larger control actions than the other control laws. This is interpreted as the introduction of probing to achieve better estimates. The problem in finding the optimal solution has motivated research in developing approximations to the loss function or other ways to change the controller such that the dual characteristic is retained. This is the so-called suboptimal dual control. Introduction 9 1 .4 Suboptimal Dual Control Among the suboptimal approaches reported in the literature, two classes may be differ-entiated. There are those that try to approximate the optimal solution, and those that attempt to emulate the dual property by explicitly forcing the controller to do what the optimal controller does implicit ly. Approaches that approximate the Bel lman equation are, however, s t i l l quite complex and require large computational efforts in real-time. In the latter class, the loss function is usually reformulated in two ways, one, to make the problem tractable and two, to introduce an active learning feature. For surveys of subop-t imal approaches, see [29, 69]. A n interesting comparison between several non-dual and suboptimal dual controllers is found in [2]. A n excellent analysis and comparison of some suboptimal approaches is provided in [43]. A simple, yet effective, suboptimal approach, that falls into the second class mentioned above, is the so-called active suboptimal dual (ASOD) controller suggested by Wit ten-mark [68]. The idea behind the A S O D controller is to find the controller minimizing the loss function J = E{(y(t + l)-yT(t+ l ) ) 2 + Xf(P(t + 2)) \Yt} (1.27) where A is a design parameter and P(t + 2) is the first time at which the variance is influenced by u(t) (for a delay-free process). This is an extension, on an ad hoc basis, of the one step ahead loss function. Add ing a function of the future estimation variance matrix to the loss function wi l l reward good parameter estimates and prevent the controller from turning off the control. Since the crucial parameter is assumed to be b0 in (1.1), the function f(P(t + 2))=Pbo(t + 2) (1.28) is commonly used. In general it is not possible to minimize (1.27) analytically since P(t + 2) is nonlinear in u(t). However, in [70] some approximations were made to obtain an explicit controller expression (Appendix A ) . From simulation experiments it has been shown that the A S O D controller is close to the optimal dual controller for low order systems. However, its performance is worse when applied to more complex systems [42]. This was explained analytically by showing that Introduction 10 the A S O D loss function is a good approximation of the two-stage optimal loss function when the process is of low-order and under some other conditions [43]. 1.5 When to Use Dual Control? A s Wit tenmark has pointed out [69], there are at least two situations in which dual control is appropriate or even necessary: 1. When the time horizon is short and the ini t ial estimates are poor, it is then impor-tant to get good parameter estimates quickly. Examples of this situation include economic systems and post-surgical blood pressure regulation. In the latter applica-tion, the time horizon is typically less than 200 samples and the model parameters exhibit wide variations due to drug sensitivity of each patient [46]. 2. When the parameters are varying rapidly and possibly changing sign. Bioreactors and chip refiners are probably the most well known examples of this class. Some dual control approaches have been investigated for the chip refining problem [2, 20]. The application in [2] is believed to be the first application of dual control to process control. Ch ip refining and paper coating, the two processes considered in this work, share the fact that they can be modeled as a low order system with an unknown time-varying gain that often switches sign. Moreover, there typically exists a large uncertainty at startup. Neither fixed-parameter controllers nor classical adaptive controllers are able to provide consistent adequate performance in these situations. This motivates the application of dual control to these problems. 1.6 The Thesis 1.6.1 Objectives The objective of this work is to develop dual controllers for chip refiners and paper coaters in order to improve pulp and coated paper quality. Al though there has been a fair amount Introduction 11 of academic work for designing adaptive controllers with dual features, very few industrial applications have yet been reported. This thesis aims at bridging the gap between existing theoretical work and the industrial implementation of dual controllers. Another objective, denned later through the course of this work, was to develop a novel algorithm for adaptive Kalman filtering, used for parameter estimation in the dual control strategy proposed, to resolve some of the problems encountered in practice. 1.6.2 Contributions The major contributions of this work can be highlighted as follows: • The utilization of an adaptive Kalman filter for parameter estimation in a dual control strategy. • The industrial application of dual control to a coupled multivariable process (paper coating). • The development of a novel adaptive Kalman filtering scheme based on the Expectation-Maximization (EM) algorithm. Some other minor contributions are also summarized • The development of a chip refiner control strategy that does not require ad hoc rules. • The application of the ASOD controller to a process with time delay (chip refining). • Extending the ASOD controller and the adaptive Kalman filter to a multivariable process. 1.6.3 Overview This thesis is organized as follows. Chapter 2 describes the chip refining problem, dis-cusses the proposed dual control strategy and presents results of simulation experiments. Chapter 3 is devoted to the paper coating process. Process modeling, controller design, Introduction 12 implementation aspects and paper mi l l results are the subject of that chapter. In Chap-ter 4, a new approach to adaptive Ka lman filtering is derived and evaluated through a number of simulation examples. Finally, the overall summary and conclusions of the thesis are given in Chapter 5. Chapter 2 Dual Control of Chip Refining 2.1 Introduction Wood chip refiners are extensively used to convert wood chips into pulp in thermomechan-ical pulping ( T M P ) plants where up to 250 k W h / t of electrical energy can be consumed. Effective control in this part of the paper making process wi l l • improve pulp quality, • reduce energy consumption, and • minimize the faults occurring in the whole process. A chip refiner consists of either one fixed and one rotating or two counter-rotating grooved plates, with pressure exerted on one of them by a hydraulic cylinder (figure 2.1). Wood chips and dilut ion water are fed near the axis and forced to move outward between the plates by centrifugal and friction forces. Feeding of wood chips is realized via a transfer screw and the rotating plates are driven by electrical motors. Steam produced by evaporation and chips broken down into fibers by mechanical actions create a few hundred micrometers thick pad between the plates. This pad forms a major load to the driving motor. Steam and pulp are discharged at the periphery. To characterize the quality of the produced pulp, a number called freeness, F, is used. It describes the drainage characteristics of the pulp, i.e., the resistance of a fiber mat to the flow of water. Pu lp freeness is used as an indicator of how well the pulp wi l l drain on the paper machine [2]. The adverse effect of freeness variability on papermaking is well known. However, since the freeness measurement is only taken every one or two hours in most paper mills, such a long sampling period is not suitable for online control 13 Dual Control of Chip Refining 14 Chips Plate Gap Pulp Figure 2.1: Schematic diagram of a chip refiner. purposes. Instead, another quantity is used. Freeness is most strongly influenced by the specific energy imparted in the pulping process. The specific energy is defined as energy consumed per mass unit of wood fibers and is expressed as Specific energy = motor load dry pulp mass flow rate (2.1) This means that there is a deterministic function which relates freeness and specific en-ergy. Bu t specific energy alone does not determine pulp properties. Refining consistency, defined as the mass of dry fiber divided by the total mass of pulp, also has a major effect [2]. The overall system has three control elements: chip feed rate, dilution water feed rate and plate gap, and generates three outputs: pulp throughput, specific energy and refining consistency. The system is also subject to disturbances that include varying feed characteristics and plate wear. Al though alternatives were reported [66], a typical refiner control strategy is as follows: the chip feed rate is held constant to meet production requirements, di lut ion water feed rate controls the refining consistency and the plate gap, v ia its effect on the motor load, controls the specific energy. The plate gap is adjusted by manipulating the hydraulic pressure applied to the shaft. Al though plate gap sensors are available, the common practice is to simply measure the shaft displacement to indicate plate movement. However, because of plate wear Dual Control of Chip Refining 15 and thermal expansion, this does not provide an absolute measurement of the gap. The hydraulic pressure is generally manipulated by sending a pulse to a microdial commanding a servo valve directing oil to the cylinder [56]. Al though there has been a fair amount of activity in the development of control systems for chip refiners, many refiners are not under closed-loop control. One difficulty is that the motor load control process is nonlinear and time-varying. When a linearized model is used, the model gain is subject to a slow drift as well as to sudden sign reversal indicating the onset of pulp pad collapse. This is discussed in more detail in the next section. In the past, there has been a number of different approaches to that problem. Horner and Korhonen [36] proposed to use gain scheduling to compensate for plate wear, assuming exponential decay of the gain and relying on operator judgment for the severity of plate clashes. Plate gap sensors and plate clash detectors using vibration monitors have been used in conjunction with heuristic logic to try to detect pad collapses. Despite being very useful as last resort safety devices, they are not used much in practice. A first attempt at applying adaptive control to the chip refiner was made by D u -mont [18], using a self-tuning regulator consisting of a recursive least-squares parameter estimator wi th a variable forgetting factor [30] and a Dahl in controller. A set of rules was used to back out the plates and restore the pulp pad in the case of pad collapse. Al though some success on an industrial refiner was reported, the strategy was deemed unreliable for continuous, unsupervised operation. The main problem was in adjusting the variable forgetting factor scheme to track both slow drifts and abrupt changes [19]. Dumont and Astrom [20] investigated several alternatives including active suboptimal dual control with a nonlinear performance index chosen to reflect the nonlinear behavior of the process. Their technique proved reliable and performed well in simulations. Their work, however, was l imited to simulation studies of simple, delay-free systems with white noise disturbances. Al l i son [2] developed what is probably the first application of dual control to process control. The controller is an active adaptive controller, which consists of a constrained certainty equivalence approach coupled with an extended output horizon, a cost function to get probing and some heuristic logic to deal with the nonlinearity. Dual Control of Chip Refining 16 More recently, Owen et al . [49] - t ry ing to avoid the complexity of an adaptive controller-implemented a fixed-parameter PI controller which is disabled when a simple linear re-gression detects a probable gain sign reversal with sufficient confidence. The plates are then retracted for a preset amount (determined by experience). The dual feature of the controller was retained by injecting a probing signal at the process input. 2.2 Process Description Owing to the complicated nature of the process, the characteristic curve between the motor load and the plate gap is both nonstationary and nonlinear. The nonstationarity is due to the variability in the feed characteristics, such as wood species, chip size and density as well as to plate wear. As the plates wear, the noise level increases and the process gain drops dramatically [51]. Wear generally occurs gradually over a period of several hundred hours but can also occur quickly in case of plate clash, i.e., metal to metal contact between the plates. In manual operation, plate clashing is a very fast phenomenon primarily related to feeding problems. It is highly undesirable as it disrupts production and damages the plates. A n idealized steady-state operating curve for a chip refiner is shown in figure 2.2. A s the gap decreases, a maximum load is reached. If the gap is reduced further, the load sharply drops as the pad -no longer able to sustain the pressure- collapses. This is believed to correspond to the point where the refiner starts shifting from fiber surface development to fiber cutting. Compared with the friction, breakage needs much less energy. As a result, a further decrease in the gap wi l l lead to more fibers being broken and less energy being consumed. The incremental gain between motor load and plate gap has indeed reversed its sign. Then, for the pad to rebuild and the gain to become negative again, the gap has to be opened past the point where the collapse occurred. This corresponds to the hysteresis pattern in figure 2.2. During refiner operation, the control objective is to regulate the motor load while operating solely in normal operating zone (to the right of the dashed line in figure 2.2) where closing the gap increases the load. To avoid plate clashes and poor quality pulp, a controller should never attempt to operate in the pad collapse zone. Product ivi ty and Dual Control of Chip Refining 17 6.6 1.5 i 2 Plate gap, mm 2.5 3 Figure 2.2: Motor load versus plate gap. pulp quality considerations may dictate operating the refiner close to the maximum load where pad collapse is always a possibility. Further complications corne from the facts that a significant amount of noise is present and that the curve in figure 2.2 depends on feed characteristics, plate wear, etc. and thus, the maximum load and the collapse point vary and are unpredictable. 2.3 Process Model ing The chip refiner dynamics are essentially due to the hydraulic system. Thus, the chip refiner can be modeled by a discrete linear system with an output nonlinearity as described by the following equations. where A = 1 — q"1 is the increment operator, A n is the duration of the pulse sent to the hydraulic cylinder, x is the plate gap, a is the process pole, and d is the dead time. The load y is given by Ax(t) k(l + a) 1 + aq-1 Au{t - d) (2.2) y(t) = g(x(t)) (2-3) Dual Control of Chip Refining 18 where g(x) represents the load-gap nonlinearity. Then, Ay{t) = ^ W) (2.4) This linear model is a reasonable approximation of the actual system for small variations around the operating point. Using (2.2), we can write Ay(t) = ~aAy(t - 1) + b(x)Au(t - d) (2.5) w here b{x) = k(l + a ) ^ - (2.6) Typica l values of the dead time and the time constant were found to be respectively 4 sec and 7 sec on a double disc refiner [51]. For a sampling interval Ts = 2 sec, this gives Ay(t) = 0.75Ay(t - 1) + b(x)Au(t - 3) (2.7) W i t h exact knowledge of k, x and g(x), the control of the system (2.7) is t r iv ia l . Unfor-tunately, not only is the exact nonlinearity unknown, but it is also time-varying. A simplistic approach to this problem is to ignore the nonlinearity and use a fixed-parameter linear controller. As long as a reasonable gain is used for the particular op-erating point, the performance is satisfactory. However, if, for some reason, the refiner enters the zone where the gain is positive, the closed-loop system becomes unstable. This situation is inevitable if the setpoint were made greater than the maximum load, in which case the controller would actually accelerate pad collapse and induce plate clashing as a result of the gain sign change associated with traversing the crit ical gap. A solution to prevent plate clashing in case of pad collapse could consist in opening the gap when it is below a given value. However, the collapse point is riot constant and is rather unpre-dictable. Moreover, this can induce a cyclic behaviour, creating oscillations between the open position and the lower gap l imit . What is really needed is a control scheme that can recover on its own from a pad collapse without plate clash and that also tracks the slow drift in gain due to plate wear. Dual Control of Chip Refining 19 2.4 Parameter Estimation The linearized stochastic model of a typical chip refiner is (1 + aq-l)ky{f) = b(t)Au(t - 3) + (1 + cg _ 1)e(t) (2.8) where e(t) is a white noise sequence N(0, v) and c is the noise filter zero which is assumed known. B y defining the filtered output v(t) as v(t) = l ± ^ A y ( t ) (2.9) 1 + cq 1 and the state variable s(f) as s{t) = l-^Tb(t)Au(t - 3) (2.10) 1 I Then, (2.8) can be rewritten as v(t) = s(t) + e(t) (2.11) Assuming that the gain b is described by b(t + 1) = b{t) + w(t) (2.12) where w(t) is additional white noise N(0, p), s(t + 1) can then be expressed as s(t + 1) = -cs{t) + Au(t - 2)b(t) + Au(t - 2)w(t) (2.13) Using (2.11)—(2.13), the chip refiner is represented by the time-varying state-space model x(t + l) = F(t)x{t)+j(t) (2.14) v(t) = Hx(t) + e(t) (2.15) Dual Control of Chip Refining 20 where x(t) is the state vector x(t) = b(t) s(t) (2.16) F(t) is the state transition matrix F(t) 1 0 Au(t - 2) - c (2.17) 7(r) is the state noise vector 7(*) w{t) Au(t - 2)w(t) (2.18) and H is the output matrix H = 0 1 (2.19) The process model, defined by (2.14) and (2.15), looks at the gain b as being time-varying and ignores the fact that it is nonlinearly dependent on the gap x. The gain b can then be estimated using a standard Ka lman filter. e(t) = v(t) - Hx(t\t - 1) K{t) = P(t\t - 1)HT (y2 + HP(t\t - 1 ) H T ) - 1 x(t\t) = x(t\t- 1) + K{t)e(t) x(t + l\t) = F{t)x(t\t) P(t\t) = (I - K(t)H) P(t\t - 1) P(t + l\t) = F(t)P(t\t)FT{t) + Q{t)p2 where Q(t)p2 = E 7 T W } 1 Au(t - 2) Au(t-2) Au2{t-2) (2.20) (2.21) (2.22) (2.23) (2.24) (2.25) (2.26) v2 and p2, as defined earlier, represent the variances of the noise sequences e and w Dual Control of Chip Refining 21 respectively. The performance of this estimator is quite influenced by the choice of p [20]. A small value for p allows the tracking of a slow drift but wi l l impair the speed of identification in the case of a sudden sign reversal. A large value for p, on the other hand, may lead to a "blow-up" of the estimates. It was suggested that p might be time-varying due to the fact that the gain undergoes two types of changes: slow drifts due to plate wear and sudden jump changes due to pad collapse. Since no apriori knowledge is available, direct estimation of p2 w i l l be carried out. In Ka lman filtering, the problem of estimating the noise covariances is known as adaptive filtering. Adaptive filtering was a very active field in the late sixties and early seventies. For a survey of different approaches see [48]. A n enhanced version of Isaksson's adaptive Ka lman filter [38] wi l l be used. Isaksson's method is based on expressing e2(t) as a linear function of p2 and u2 and using a least-squares (LS) approach to estimate the noise variances. Here, this method is augmented with two more linear regressions derived for e(t)e(t — 1) and e(t)e(t — 2) to boost the accuracy of the adaptive filter. We now derive the adaptive filter scheme. If we define the estimation error as We obtain, using (2.14), (2.22) and (2.23), x(t) = F(t - 1) (I - K(t - 1)H) x{t - 1) + 7(t - 1) - F{t - l)K(t - l)e(t - 1) (2.28) Now, define the conditional variance of x(t) as where Yt-\ denotes al l data observed up to and including time t — 1. Note that, since the true values of p and v are unknown, the Ka lman filter is not optimal, and so, P(t\t — 1) x(t) = x(t) - x(t\t - 1) (2.27) (2.29) Dual Control of Chip Refining 22 does not represent the actual error variance M(f). Using (2.28), M(t) is evaluated M(t) = F ( * - 1 ) ( I - K[t-l)H)M[t-l){l- K{t-1)H)T FT(t-l) + Q(t - l)p2 + F(t - l)K(t - l)KT(t - l)FT(t - l)u2 (2.30) Starting at time t0, we may rewrite (2.30) as M(t) = f(M(t0),t)+T(t-l)p2 + Z ( t - i y (2.31) where T(t) and Z(i) can be computed recursively as T(t) = Q(t) + F(t) (I - K{t)H) T(t - 1) (I - K(t)H)r FT(t) (2.32) Z{t) = F(t)K(t)KT(t)FT(t) + F(i) (I - K{t)H) Z(t - 1) (I - K{t)H)T FT(t) (2.33) If we have a stable Ka lman filter 1 , the first term in (2.31) wi l l tend to zero as t increases. Hence, we may neglect the influence of the ini t ial values and M{t) is then a known linear function in p2 and u2 M(t) = T{t - l)p2 + Z(t - l)v2 (2.34) We can use this fact in the following way. The prediction error, defined in (2.20), can be expressed as e(t) = Hx(t) + e{t) (2.35) We get the conditional variance of e(t) E [e2{t) | y t _ i } = HM(t)HT + v2 (2.36) Hence, it is possible to approximate e2(t) as e2(t) = [HT(t - 1)HT] p2 + [HZ(t - 1)HT + l ] v2 (2.37) 1 The Kalman filter is always stable provided that the process dynamic model is correct. Dual Control of Chip Refining 23 In the same way, two additional linear regressions are derived for the correlation functions E{e{t)e(t - 1) \Yt_2) and E{e(t)e(t - 2) \Yt-3}- Then, e(t)e(t - 1) = [HF(t - 1) (I - K(t - l)H) T(t - 2)HT] p2 + HF{t - 1) [(I - K(t - 1)H) Z(t - 2)HT - K(t - 1)] v2 (2.38) and, e(t)e(t-2) = [HF(t-l)(I- K(t - 1)H) F(t - 2) (I - K(t - 2)H)T{t - 3)HT] p2 + HF(t - 1) (I - K(t - 1)H) F{t - 2) • [(I - K(t - 2)H) Z(t - 3)HT - K(t - 2)] v2 (2.39) The implemented R L S estimator utilizes the three linear regressions (2.37), (2.38) and (2.39). A s suggested by Isaksson, the estimates p and v are not used to improve the K a l m a n filtering of the gain bit) until they are likely to have converged. This is carried out in the simulations by waiting 500 sec before using p and 0 in the Ka lman filter. 2.5 Controller Design Most control laws suggested so far for the refiner motor load control process are al l linear in Au and, because of that, require the addition of heuristic logic to handle pad collapses. Al though the exact nonlinearity is not known, its general shape is known apriori. It is thus natural to use a performance index that yields a nonlinear control law [20]. A n important feature of the control criterion is that control in the positive gain zone should be severely penalized if not prohibited. When the gain estimate is negative and is relatively of large amplitude, closing the gap should obviously be allowed. However, as the estimate approaches zero, negative control actions should be increasingly penalized. Furthermore, when the estimate is close to zero and the uncertainty is large, probing should not tend Dual Control of Chip Refining 24 to further close the plates. Given al l these considerations, the loss function is chosen as J = E{{y{t + 3)-yr(t + 3))2 + toxh{b{t + 3 ) ) / ( A M ( * ) ) + u2Pb{t + 4\t) \Yt } (2.40) where yr is the setpoint, ui and o>2 are positive weighting constants and the variance Pb Pbs matrix P is represented by P = Pbs Ps The minimization is carried out with respect to Au(t). The function h(b) should be of small amplitude when b is negative and should increase as b approaches zero and then becomes positive. As for f(Au), it should penalize closing the gap, i.e., it should be high when Ait is negative. Various choices are possible, but the functions below satisfy the preceding requirement while being simple [20]. h(b) = eb (2.41) f(Au) = e" - A M (2.42) The future variance Pb(t + 4\t) is included in the loss function -following Wittenmark 's A S O D controller [68]- as an artificial way to introduce probing. Using (2.25) and the inversion form of (2.24), we may write P{t + 2\t) = F(f+l) - l HTH + P~\t + l | i ) FT{t + 1) + Q(t + l ) p 2 (2.43) Then, P(t + 3\t) = F(t + 2) ;HTH + p-l(t + 2\t) FT(t + 2) + Q(t + 2)p2 (2.44) Knowing Au(t- 1), P(t + 2\t) is determined. Al though P(t + 3\t) is dependent on Au(t), Pb{t + 3|t) is not. Pb{t + 4\t) is given by Pb(t + 4\t) 1 0 ;HrH + p-l(t + 3\t) 1 0 (2.45) Pb(t + 4|t) is thus a deterministic function of Au(t). Pb(t + 4\t) penalizes inaction in the presence of large uncertainty, so, preventing the turn-off phenomenon. In fact, the effect Dual Control of Chip Refining 25 of Pb(t + A\t) is to increase the magnitude of the optimal control that minimizes the loss function J. This excess in magnitude can be viewed as a probing component to enhance the quality of future estimates. Using (2.8), y(t + 3) can be written as y(t + 3) = ( l - a + a 2 - a 3 ) y(f) + (a - a 2 + a 3 ) y{t - 1) + b(t + 3)Au(t) + b(t + 2)(1 - a ) A u ( i - 1) + b(t + 1) ( l - a + a 2 ) Au(t - 2) + e(t + 3) + (1 - a + c)e(t + 2) + ( l - a + a 2 + c ( l - a)) e(i + 1) + c ( l - a + a 2 )e( t ) (2.46) Provided that a and c are deterministic, the probability distribution of y(t + 3) given Yt is Gaussian. (2.11) and (2.12) directly lead to E{e(t)\Yt} = v(t)-s(t\t) (2.47) and, b(t + 3\t) = b{t + 2\t) = b(t +l\t) = b(t\t) (2.48) Using (2.46)-(2.48) and the facts: E {z2} = p2 + a2 and E{e*} N ( / i , a ) , the expectation in (2.40) is evaluated where z is bit + l\t)Au(t) + b(t + l\t)(l - a)Au(t - 1) + b{t + l\t) (1 - a + a2) Au(t - 2) J = + (1 - a + a 2 - a 3 )y( t ) + ( a - a 2 + a 3 ) y ( £ - 1) + c (1 - a + a 2) (>(£) - s(«|t)) - y r ( t + 3)] 2 + Pb(t + 3\t)Au2(t) + Pb(t + 2 | t )( l - a)2Au2(t - 1) + Pb(t + l | t ) ( l - a + a 2 ) 2 A u 2 ( t - 2) + 2Pb(t + 2 | t ) ( l - o )Au( t )Aw(t - 1) + 2Pb{t + l | t ) ( l - a) ( l - a + a 2 ) Au(t - l ) A u ( t - 2) + 2Pb(t + l\t) ( l - a + a 2 ) A u ( t ) A u ( t - 2) - 2c ( l - a + a 2 ) 2 P 6 , ( t | t )Au( t - 2) - 2c( l - a) ( l - a 4- a 2 ) Pbs(t\t)Au{t - 1) - 2c ( l - a + a 2 ) P 6 s ( t |* )Au( t ) 1 + (1 - a + c) 2 + ( l - a + a 2 + c ( l - a ) ) 2 + [c ( l - a + a 2 ) ] 2 P s ( t | t ) 6 ( l + l | t)+ift ( l + 3 | t ) e - A « ( t ) + W 2 p ^ t + 4 | t ) ( 2 4 9 ) W i e Dual Control of Chip Refining 26 Typica l curves of this performance index for several values of b, P and for both a positive and negative control error are plotted in figures 2.3 and 2.4. Some comments on the loss function and the resulting control law are now in order: 1. The influence of the nonlinear function eMH-3) e-Au(t) is negligible when the estimator is confident that the refiner is operating well inside the desirable zone. A s the confidence decreases and the gain estimate increases, it is seen that the optimal control is pushed towards positive values. So, The resulting control law tends to open the gap when closing it is likely to induce a collapse. 2. The resulting control law wi l l never attempt to control the refiner in the pad collapse region. 3. The use of probing is justified by the fact that the gain estimate is the key to identifying a pad collapse and that probing targets a portion of the input energy at continuously identifying that parameter. 4. No heuristic logic is needed, which makes the resulting controller more appealing for implementation. 5. The optimal control can not be computed analytically. For simulation purposes, it is obtained through a Newton-Raphson numerical optimization scheme. 2.6 Simulations The behavior of the strategy was thoroughly tested by way of computer simulations. The process was simulated by (2.2) and (2.3) with the output nonlinearity g(x) represented In normal case a = 1. In collapsed state, i.e., when x < 0.9 mm, a = 0.77. Once in collapsed state, the refiner stays in that state until the plates are opened past the point where the two curves intersect (see figure 2.2). Max imum load is achieved at x = 1 mm. For both curves gi = 7600 and g2 = 800. by (2.50) Dual Control of Chip Refining b=-o.l 4000 3800 3600 q$ 3400 C § 3200 £ o t C L 3000 2800 2600 2400 -0.5 0 0.5 Control signal b=0.1 \p=o. D6 ^ ^ 0 03 P^OTCr i i i -1.5 -0.5 0 0.5 Control signal 1.5 Figure 2.3: The performance index J(Au(t)) (load below setpoint). Dual Control of Chip Refining b=-o.i I 1 i i \=0.06 - V \P=0 .03 P^ CVOT : i i i i --1.5 -0.5 0 0.5 Control signal b=0.1 1.5 I i i f?=b .06 V=0.03 -1.5 -0.5 0 0.5 Control signal 1.5 Figure 2.4: The performance index J(Au(t)) (load above setpoint). Dual Control of Chip Refining 29 Other process settings are: process pole a = —0.75, noise filter zero c = —0.99, Ts = 2 sec, v = 0.5, and ini t ial gap x0 = 1.6 mm. The control signal A u is l imited to the region [—2, 2] mm. The controller parameters are: u>i = 12 and u>2 = 100. Adaptive K a l m a n filter in i t ia l estimates p0 = 0.05 and u0 = 0.4. The controller is not too affected when there exists an error in the modeling of the noise filter zero c. The following simulations were all carried out with the controller having a value of —0.95 chosen for c. Figures 2.5 and 2.6 demonstrate normal operation, i.e., when the refiner is operating well into the negative gain zone. The control problem is relatively easy in this case as the gain varies li t t le. This is reflected by p being very small. For 1000 < t < 2000 sec, the setpoint is raised by 250 k W to an unreachable value. If the controller tries to achieve the new setpoint, it wi l l induce a pad collapse. In both trials, the controller performs very well as it does not provoke any collapse. It keeps the load as close as possible to the unreachable setpoint while not closing the gap below the collapse point. The adaptive K a l m a n filter successfully tracks the rapid variations of the gain as the gap alternates around the maximum load point x = 1. In figure 2.7, the refiner is required to operate near maximum load, a situation pre-ferred in practice in order to maximize production but that carries a high risk of collapse. The controller successfully maintains the motor load close to the setpoint while avoid-ing any collapse. Again , the controller behaves very well as the setpoint is raised to an unreachable value for 1000 < t < 2000 sec. Final ly , figure 2.8 depicts the occurrence of a pad collapse. The adaptive K a l m a n filter is very fast to react as seen from the p and the gain estimate graphs. The setpoint is reduced shortly after the gain sign has become positive. This situation (load above setpoint and positive gain) usually poses a challenge to any chip refiner control scheme as the controller wi l l tend to further close the gap i f it directly tries to achieve the new setpoint. However, as evident from the gap plot, the controller opens the gap, recovers from the pad collapse and brings the load to the new safe setpoint. Al though the controller is designed so as not to induce a pad collapse, it did not succeed in avoiding the pad collapse in figure 2.8. In fact, the proposed strategy does not guarantee absolute prevention of pad collapse. The occurrence of a pad collapse is Dual Control of Chip Refining 30 always a possibility when the refiner is operating at the extremum of the nonlinearity (due to unreachable setpoint, for instance). However, the controller has the advantage that, in the case of pad collapse, it can rapidly recover on its own. In order to guarantee prevention of pad collapse, constraints have to be added for the setpoint high l imit and the gap low l imit . Yet, in light of the characteristic curve nonstationary nature, these constraints would have to be set conservatively and may thus jeopardize performance. 2.7 Conclusions This chapter has demonstrated that it is possible to obtain a reliable controller for the chip renner motor load (in the presence of a time delay). The use of an adaptive K a l m a n filter for parameter estimation proved to be worthwhile as it was capable of tracking both slow drifts and abrupt sign changes in the gain. The major advantages of the proposed approach are that it does not require additional heuristic logic or expensive dedicated sensors and that, in the case of pad collapse, the controller rapidly recovers on its own. The disadvantage, however, is that the controller is obtained through numerical optimization. Economic concerns at the time of this work made it very difficult to find a mi l l to test the developed controller. Dual Control of Chip ReGning 31 Figure 2.5: Safe operation, then unreachable setpoint. Dual Control of Chip Refining 32 Figure 2.6: More of normal operation. Dual Control of Chip Refining 33 Figure 2.7: Operation around zero gain. Dual Control of Chip Refining Figure 2.8: Behavior of the controller when the pad collapses. Chapter 3 Dual Control of Paper Coating 3.1 Introduction Coating has been applied to paper for over 100 years. L I F E magazine, born in 1936, was the first to use al l coated paper in a magazine. Today, almost one out of every four tons of printing paper produced in North America is coated [33]. Coating is clearly a form of surface enhancement. Printabi l i ty of paper is improved by leveling the face of the sheet and impart ing a high degree of smoothness and gloss to it. Addit ionally, brightness and opacity are increased where needed. Blade coaters [63] are firmly established as the most popular means of producing high quality coated paper. They are becoming more and more sophisticated and there are numerous types of construction, which differ in the mode of application of the coating material and in the leveling and coat weight adjustment mechanism. The uniformity of the coating material on the paper is extremely important from both a product quality viewpoint and material savings viewpoint. Increased demand by cus-tomers and end-users for better quality coated paper has forced producers to continually improve their coating processes. Many have taken a serious interest in automatic coat weight control and some are already adopting the technology. Coat weight uniformity can be improved indirectly, by manipulating a base sheet property [64], or directly by manipulating some aspect of the coating process where coating is applied to a continuously moving web. Because coat weight varies both across the sheet and along the sheet, direct coat weight control is a two-dimensional control problem and is carried out separately in the machine-direction (MD) and the cross-direction (CD) . M D is the direction parallel to the motion of the sheet, while C D refers to the direction perpendicular to the sheet travel. 35 Dual Control of Paper Coating 36 M D control is concerned with controlling the average coating of each scan, where the average is taken across the width of the sheet. The role of the M D controller is constrained (by the nature of the actuator assigned to it) to affect the entire width of the sheet simultaneously, though perhaps not uniformly. Achieving an even distribution of coating across the width of the sheet is the task of the C D controller and is handled by an array of actuators arranged across the width of the coating station. Automat ic M D coat weight control is well known in the industry 1 [1, 25, 52, 53, 54], but C D coat weight control is s t i l l relatively new. Over the last several years, coated paper producers have begun to implement C D coat weight control with promising results [24, 34, 35, 55, 59]. Most ofthe literature available, however, concentrates on issues associated with the implementation and operation of C D coat weight control systems rather than the design of control algorithms. This chapter addresses the C D control problem on bent blade coaters. O n bent blade coaters, the C D coat weight profile is regulated by manipulating the blade pressure in zones across the blade. The process gain varies over time due to blade wear and changing operating point and often switches sign. This poses a problem to any fixed-parameter linear controller which wi l l get caught in a positive feedback loop as the gain reverses its sign. Nevertheless, closed-loop control of C D coat weight is operational in numerous mills. A common industrial practice is to control C D coat weight with a tradit ional fixed-parameter Dahl in controller with a decoupling matrix with the aid of some heuristic logic that simply switches off automatic control when it detects that the loop has become unstable. The coater is then running in open-loop. Whi le this technique has proven useful, there is clearly a wide room for improvement. The time-varying nonlinear relationship between C D coat weight and blade pressure requires the use of some form of advanced control. A couple of research reports [32, 65] claim to have successfully implemented adaptive controllers for C D coat weight control, but with the absence of any technical details no assessment of their work can be done. 1 Commonly performed on blade coaters by adjusting the blade angle. Dual Control of Paper Coating 37 Figure 3.1: Blade coater [63]. 3.2 Blade Coating A typical bent blade coater is illustrated in figure 3.1. A n applicator roll running in a pan of coating applies excess coating to the moving paper web as it passes through a flooded nip formed by the applicator roll and the backing roll which supports the sheet. Excess coating is removed with a thin steel bent blade pressed against the paper and backing roll (figure 3.2). The coating material is essentially troweled off the base sheet with the side of the blade. The bent blade is flexible, square edged, and is engaged at low angles, usually below 10°. The blade can be loaded either by a solid backing bar or by a pneumatic pressure inside an elastic air tube. Coat weight applied by a blade coater is influenced by a number of process parameters, the most important ones of which are base paper quality (roughness, moisture), coating material viscosity, blade geometry, blade angle and blade loading. A l l these factors are also present when local coat weight changes are discussed. The effects of these factors have been studied by researchers over a period of about thirty years. During that time, Dual Control of Paper Coating 38 Figure 3.2: Schematic diagram of a bent blade coater [37]. many models have been presented, combining these factors in different ways [10, 21, 31, 61, 62, 67]. 3.2.1 Coat Weight Blade Pressure Relationship O n a bent blade coater, the amount of coating applied is dependent on an equilibrium between the blade pressure and the hydrodynamic forces from the coating color acting against the blade. Because the blade is flexible, increased blade pressure bends the blade further, increasing the area of the blade tip pressing against the paper sheet. Since more blade surface area is exposed to the coating color, the blade is forced away from the web by the coating hydrodynamic pressure and coat weight increases [31]. However, the increase of the blade pressure eventually surpasses that of the hydrodynamic pressure, and coat weight starts decreasing after reaching a maximum. A t this point, the blade pressure is so great that it starts forcing the blade against the web. Simultaneously, the structure of the coating changes drastically from being a thick smooth layer to a coating Dual Control of Paper Coating 39 5 ro Actuator position Figure 3.3: CD coat weight as a function of actuator position. layer with a pattern similar to that of a roll coater. This phenomenon is often called film splitting [21]. Moreover, the rate of change of coat weight is not linear but has been found to be much higher at higher blade pressures [31]. The static characteristic between CD coat weight and actuator position has the general profile of figure 3.3 [21, 23]. When blade loading is changed, it leads to a change in coat weight in a way described by this nonlinear curve. Further complication comes from the fact that this relationship is nonstationary since it depends on coater speed, coating material, coated paper grade, blade angle, blade wear as well as CD position. Bent blade coating normally operates in the negative gain zone2, where an increase in blade loading leads to a higher coat weight (figure 3.3). Positive gain zones designate stiff blade operation [21, 31, 40]. The major difference is that in bent blade operation the blade trowels -rather than scrapes- the coating with the side of the blade (figure 3.4). The bent blade is lifted by the coating hydrodynamic pressure and is not necessarily in contact with the sheet. The mode of operation is related to the blade angle (a in figure 3.2). Bent blade oper-ation is characterized by a much smaller blade angle than stiff blade (figure 3.4) and hence 2The industry standard is that a positive actuator movement reduces the local pressure. Dual Control of Paper Coating 40 Bent blade Stiff blade Figure 3.4: Bent blade vs. stiff (bevelled) blade. it is often referred to as low-angle coating. In addition to the geometrical differences, the blade configurations exhibit various coating behaviors. Stiff blade mode is more conve-nient for achieving light coat weights, usually from 5 to 15 grams per square meter (gsm), and bent blade mode for thicker ones (15-35 gsm). Moreover, Bent blade operation offers improved optical properties as compared to stiff" blade [26] and allows coating of smooth substrates, something which is impossible with the stiff blade coater [22]. As Ek lund [21] and other researchers have observed, there is a continuous progression from stiff to bent blade behavior in coating. 3.2.2 Actuators and Measurements C D coat weight actuators are generally electric or pneumatic stepper motors [17, 37]. They are closely spaced (75-150 mm) across the blade as demonstrated by figure 3.5. A profile actuator pushes or pulls the loading element so that the local pressure is altered. If, for example, a solid backing bar is pulled away from the blade at some position along the sheet, the coat weight is locally decreased (bent blade operation). The response of one actuator typically spreads out wider than its corresponding area in the sheet. Since the responses overlap, the actuators can not be controlled independently of each other. Therefore, the C D control process forms a coupled multivariable problem. Dual Control of Paper Coating 11 Figure 3.5: C D coat weight actuators [58]. A number of techniques for on-line measurement of coat weight have been imple-mented [41]. Determination of the applied coat weight usually involves measuring the properties of the sheet before and after the coating station by gauges that scan across the sheet. Because the web is moving rapidly in the machine-direction while the gauges are moving in the cross-direction, the web is actually sampled along a diagonal path, i.e., the "raw" measurements obtained from a scanning gauge contain both C D and M D compo-nents (figure 3.6). A number of approaches have been applied to process the signal from the scanning gauge so as to separate these two components [66]. A widely used approach is to apply exponential smoothing to successive scans and to use this as an estimate of the C D profile. In practice, the number of measurement points is usually greater than the number of actuators. For C D control purposes, the measurement data are mapped to the cor-responding actuator zones. There are several mapping methods available [8, 13], the simplest of which is to take the average of the measurements in a zone. After mapping, the number of points in the low resolution coat weight profile is the same as the number of actuators. A block diagram showing the functional components and the flow of information in Dual Control of Paper Coating 42 Pure MD variation CD •4 • Figure 3.6: Diagonal scan path introduces M D variations into (what should be flat) measured profile. the C D control loop is given in figure 3.7. 3.3 Modeling of the C D Control Process The model of the C D control process has to encompass the dynamics of both the coating station and the exponential smoothing filter. Since the response of a coating station to a change in an actuator setpoint is known to be almost instantaneous, the coating station Specification CD controller Actuator setpoints Actuators Coating Scanning • station w sensor Low resolution CD profile Alignment & mapping High resolution CD profile CD/MD decoupling Raw measurements Figure 3.7: Flow of information in the C D control loop. Dual Control of Paper Coating 43 u actuator positions ,-1 white noise disturbances A l + cq 1 1 _ „ - l <±> 1 + 1 + a q -1 C D coat weight profile Figure 3.8: Block diagram of the CD control process. can be described as a static nonlinear function. Thus, the process can be modeled by a multivariable discrete-time linear system with an input nonlinearity as illustrated in figure 3.8 and described by the following equation. y(t) = n ^ r / ( « ( * - i ) ) + f cq-lAe(t) aq-1 A (3.1) where A is the increment operator, u are the actuator positions, y is the CD coat weight profile, a is the exponential smoothing filter pole, c is the noise filter zero, e are.indepen-dent Gaussian white noise sequences E |e(t)eT(t)| = z / 2 I , A is the matrix describing the correlation between process noise observed at adjacent measurement points within a scan and /(•) represents the unknown coat weight-blade pressure nonlinearity (figure 3.3). Then, the linearized model [1 + ag_1)Ay(*) = G(t)Au(t - 1) + (1 + cg_ 1)Ae(i) (3.2) is a reasonable approximation of the actual system for small variations around the operat-ing point. G is a square matrix that describes the spatial response (i.e., the steady-state CD response) of the actuators. Each column of J T ^ G contains the effect on the coat Dual Control of Paper Coating 44 weight due to a change in the position of one of the actuators, once the transient response has decayed away. G is band-diagonal due to the fact that, in coating, the response to an actuator movement has a concentrated and localized area of influence [14, 58]. The spatial response of the actuator is usually determined via "bump" tests, where step changes are applied to the setpoints of several actuators across the sheet while the process is running in open-loop (i.e., the CD controller is disabled). By observing the effect of these changes, the spatial response can be estimated. The following assumption, which is supported by the results of bump tests performed at different times throughout the blade life, will be made: the shape of the spatial response of an actuator does not change over time (in essence, this means that process conditions simply scale the spatial response without changing its shape). This assumption allows the off-diagonal elements of G to be expressed as a fixed fraction of the process gains. As an example, the following equation demonstrates the special structure of G when an actuator response is limited to its two nearest neighboring zones. G(t) = (3.3) 6i(<) ab2{f) 0 0 abi(t) b2(t) ab3(t) 0 ••• 0 0 ab2{t) b3{t) ! :. 0 ab3{t) '•• 0 ! 0 ••• abn(t) 0 0 ••• bn(t) where, n is the number of actuators, a is a coupling coefficient and the gain 6, is propor-tional to the slope of the nonlinear relationship between coat weight and blade pressure (figure 3.3) at the local operating point (actuator position) for each CD position. Taking advantage of this assumption, G is expressed as G(t) = CB{f) (3.4) Dual Control of Paper Coating 45 where C is the n x n constant actuator coupling matrix, for the previous example C = 1 a 0 • • 0 a 1 a 0 0 a 1 • a 0 0 1 and B is the n x n diagonal process gain matrix B(t) h(t) o 0 b2(t) 0 0 0 0 bn(t) (3.5) (3.6) A unique challenge presented by CD coat weight control is that the gains bi vary due to blade wear and operating conditions and may differ in sign from one CD position to the next and over the life of the blade. Another challenge is due to the subtle loading and bending limitations imposed by the actuators and the blade. Finally, relative to the paper machine, coating is a discontinuous process as it is often interrupted for blade changes. 3.4 Parameter Estimation Before discussing the dual control approach to this problem, an adaptive Kalman filter is designed for parameter estimation. Assume that the gain bi can be modeled by a Gauss-Markov process bi(t + 1) = bi(t) + Wi(t) (3.7) where Wi(t) is a Gaussian white noise sequence Dual Control of Paper Coating 46 By denning the filtered output Vi as and the state variable S j as si(t) = ——^bi(t)Aui(t-l) 1 + CO 1 the coater of (3.2) can be represented by a time-varying state-space model x{t + 1) = F(t)x(t) + j{t) v(t) = ifx(t) + Ae(t) where x(t) is the 2n-dimensional state vector x(t) hit) ••• bn(t) s i ( i ) • • • sn{t) F(t) is the 2n x 2n state transition matrix F(t) 1 0 0 1 0 0 0 Aui(t) 0 0 Au 2(t) 0 0 0 0 0 1  ••• 0 0 - c 0 ••• 0 0 0 - c 0 Aun(t) 0 - c 7(t) is the 2n-dimensional state noise vector (3.8) (3.9) (3.10) (3.11) (3.12) (3.13) 7(*) wx{t) ••• wn(t) Aui(£)iui(t) ••• Aunit)wn(t) (3.14) Dual Control of Paper Coating 47 and Pf is the n x 2n output matrix H 0 C The gains bi can then be estimated using a standard Kalman filter. e(t) = v(t) - Hx(t\t - 1) K(t) = P(t\t - 1)HT (HP(t\t - 1)HT + ^ 2 A A T ) _ 1 x(t\t) = x(t\t- 1) + K(t)e(t) x(t + l\t) = F(t)x(t\t) P(t\t) = (I - K(t)H) P(t\t - 1) Pit + \\t) = F(t)P(t\t)FT(t) + Q(t)p2 where Q(t)p2 = E{ 7 ( t ) 7 T ( t )} 1 0 0 Aui(t) 0 0 Au2(t) 0 Awi(t) 0 0 0 Au2(*) 1 0 0 Au2(t) 0 0 0 Aul{t) Aun(t) 0 0 0 Aun{t) 0 0 (3.15) (3.16) (3.17) (3.18) (3.19) (3.20) (3.21) (3.22) Following the same procedure as in Chapter 2. A multivariable version of Isaksson's adaptive Kalman filter [38] is derived. The method is based on expressing e(t)eT(t) as a linear function of p2 and u2 and using a least-squares approach to estimate the noise Dual Control of Paper Coating 48 variances. Using (3.10), (3.11) and (3.16)—(3.19), an approximate expression for the conditional variance of the prediction error E{e{t)eT(t) \Yt-i}, where Yt _ i denotes al l data observed up to and including time t — 1, is obtained E {e(t)sr{t) \Yt-! } = [HT(t - l)HT] p2 + \HZ{t - 1)HT + A A T ] v2 (3.23) where T{i) and Z(t) are computed recursively as T(t) = Q{t) + F{t) (I - K(t)H) T(t - 1) (I - K(t)H)T FT(t) (3.24) Z{t) = F(t)K(t)AATKT(t)FT(t) + F(t) (I - K(t)H)Z{t-l) (I - K(t)H)T FT(t) (3.25) Hence, it is possible to approximate e(t)eT(t) as e(t)eT{t) = [HT(t - 1)HT] p2 + [HZ(t - l)Hr + A A T ] v2 (3.26) Then, col (e{t)eT{tj) = col ( # T ( t - 1)HT) p2 + col (HZ{t - 1)HT + A A T ) v2 (3.27) where col of an m x n matrix A is the mn x 1 vector col(A) arranging the columns of A one after another in a longer column co\(A) O n a 21 " ' " flml °12 a 22 ' " ' O m 2 " ' ' a l n 0-2n " ' " ar. (3.28) Adaptive K a l m a n filtering is implemented by running two estimators in a boot-strap approach, a R L S estimator based on (3.27) to estimate the noise variances and the Ka lman filter (3.16)—(3.21) to estimate the gains. Again , the estimates p and 0 are not used to improve the filtering of the gains b{ until they are likely to have converged. This is carried out by waiting 150 scans before using p and v in the Ka lman filter. Dual Control of Paper Coating 49 3.5 Controller Design The aim of the C D controller is to improve the uniformity of the sheet by reducing coat weight variations. The choice of the loss function to minimize is crucial as it must reflect several important aspects of the control problem. A n important feature of the control criterion is that control in the positive gain zones should be penalized. To avoid actuator picketing 3 , the second difference of the actuator setpoints should also be penalized. F i -nally, the controller should often perturb the process to improve the identification of the gains at different C D positions. Given al l these considerations, the loss function is chosen as J = E{yT(t + l)y{t+l) + X1AuT{t)Au(t) + X2uT(t)MTMu(t) + A 3 E ^ i t + 1 ) ( ~ ^ - + U t ) ) 2 + A 4 E Pbi (t + 2\t) \Yt } (3.29) i m a x i ) where A i , A 2 , A 3 and A 4 are positive constants, M , sometimes called the moment gener-ating matrix, is given by - 1 1 0 0 1—» - 2 1 0 . . . o 0 1 - 2 1 0 1 - 2 1 0 0 1 -such that the elements of Mu(t) are the second differences of the actuator setpoints. A « m a x is the maximum feasible actuator position change in one scan and Q = ± 1 de-pending on whether the local operating point U{ is in the vicinity of a maximum or a minimum. Minimiza t ion of J is carried out with respect to Au(t). The first term in the loss function represents the variance of the measurements across the sheet. XiAuT(t)Au(t) weights the control effort. A 2 u T ( i ) M T M u ( t ) penalizes the 3 A C D control problem observed in practice, defined as large differences between adjacent actuator setpoints. Dual Control of Paper Coating 50 bending moment of the blade, thus, smoothing the actuator setpoint profile. ebi ( A^"' a + Ci) 2 penalizes operation in stiff blade (positive gain) mode. The flag Q is set to 1 if the local operating point Ui is closer to the nonlinearity minimum (to pe-nalize further reduction of local pressure), and to —1 if it is closer to the maximum (see figure 3.3). To determine that, the knowledge of a fixed actuator position well into the normal operating (negative gain) zone or a halfway point ^halfway is assumed. Then, - 1 Ui(t) < ^halfway ^ ^ 1 Ui(t) > ^halfway The future variance Pbi(t + 2\t) is included in the loss function in a multivariable ver-sion of Wit tenmark 's A S O D controller [68]. Neglecting the interaction between adjacent actuators, an approximate expression for Pbi{t + 2\t) can be obtained „ (f , = (KM*) + P 2 ) ( P s M t y + v2)-Plti(t\ty , 2 tow foA ' ^ v2 + ( P b M t ) + P2) A n 2 ( t ) + PSi(t\t)c2 - 2PbiSi{t\t)cAul[t)+P 1 • ' Pbi(t + 2\f) is thus a deterministic nonlinear function of Aiij(£). Using (3.2) and (3.4), J is evaluated J = AuT(t)[L(t) + X1I + X2MTM)Au(t) + 2AuT{t)GT(t+l\t)[(l-a)y(t) + ay(t - 1) + cv[f)\ + 2cAur{t)p{t) + 2X2AuT(t)MT Mu(t - 1) + A 3 e ^ 1 ^ ^ 1 ^ ^ ^ + Ut))2 + A 4 E P^t + 2\t) + J0 (3.33) j ""max j where G( t + l | t ) = E {C7(* + l\t)\Yt} = CB{t + l\t) (3.34) E { e 6 i ( t + 1 ) \Yt} = ek(t+i\t)+^ph(t+i\t) ^3>35^ L(t) = E { c 7 T ( t + l)C?(t + l ) | F 4 } (3.36) /?(£) = - E { G T ( i + l ) i fx^) |r t } (3.37) and J0 includes the terms independent of Au(t). L(t) and (5{f) are given by (see A p -Dual Control of Paper Coating 51 pendix B) and, col (L(t)) = E {B{t + 1) ® B{t + 1) \Yt} col ( C T C ) (3.38) (3.39) /?(<) = - E {s T ( t ) ® B{t + 1) | F t } col (CTC) with <g> denoting the kronecker matrix product [9]. The optimal control, which minimizes J , can not be found analytically. However, i f Pbi(t + 2\t) is expanded around a nominal input A u n o m ; , up to the quadratic term, an approximate analytical control law can be obtained. The simple way to find a suitable A u n o m i as proposed in [70] (Appendix A) wi l l be adopted here. Then, Pbi(t + 2\t) = hi(Aui(t)) = hi(Aunomiit)) + ^ ( A u n o m i ( t ) ) ( A ^ ( i ) - Aunom{t)) + h'^Aunomi(t))(Aut(t) - Aunomi(t))2 (3.40) Using (3.40), J is expressed as 4 J = AuT{t) (L(t) + A1I + A 2 M T M + A3#(£) + A4$(t)) Au(t) 1 GT(t + l\t) [(1 - a)y{t) + ay(t - 1) + cv(t)] + c(3(t) v + A 2 M T M u ( t - 1) + X3ipit) + \A<j>{t) + 2AuT{t) where e6i(t+l|t)+ip 6 l(t+l|t) 0 0 eb2(t+l\t)+^Pb2(t+l\t) Q 0 0 (3.41) 0 ebn(t+l\t)+±Pbn(t + l\t) (3. 'Omitting the terms independent of Au(t). Dual Control of Paper Coating 52 h'l(Aunomi{t)) 0 0 h'l(Aunom2(t)) 0 0 0 0 0 h'f(Aunomn(t)) Ci(t)eSl('+1l')+i^i(*+1l<) C2(t)e*2(t+1l*)+lA2(t+i|0 Cn(t)ei^t+l^+^p>>n(t+i\t) and 0(t) = 1 / i - ( A u n o m i ( t ) ) - ^ ' ( A u n o m i ( i ) ) A u n o m i ( t ) / i - ( A u n o m 2 ( i ) ) - h'l{Aunom2(t))Aunom2(t) ^ ( A n n o m n ( t ) ) - ^ ' ( A « n o m n ( t ) ) A u n o m r i ( i ) Then, the optimal control law is Au(t) - l - (L(t) + + A 2 M T M + A 3 * ( * ) + A 4 $( t ) ) " ' G T ( i + l | t ) [(1 - a)y(t) + ay{t - 1) + cv(*)] + c/5(t) v + X2MTMu{t - 1) + A 3 ^ ( t ) + A4</>(£) (3.43) (3.44) (3.45) (3.46) 3.5.1 Accommoda t ing Constraints The above control law inherently assumes that there are no constraints on the system. In practice, this is unrealistic as there are both loading and bending limitations imposed by the actuators and the blade. Specifically, 1. the setpoints that can be applied to the actuators must lie in the range [—Umax,"max], so that \ui(t)\<um&x (3.47) Dual Control of Paper Coating 53 2. setpoints changes are limited, that is, \Aui(t)\ < A * w (3.48) 3. to avoid sharp deflection of the blade, the difference between adjacent actuator positions is l imited \ui+1(t) - Ui(t)\ < 6lmax (3.49) 4. to prevent excess bending being applied to the blade, the second difference of the actuator setpoints is bounded K _ i ( i ) - 2ui{t) + ui+1(t)\ < 52max (3.50) A simple, yet effective, constraint handling method wi l l be utilized. The unconstrained setpoint changes A u u n c o n s t given by (3.46) are computed. If no constraint is violated, ^unconst is feasible and is applied to the process. When one or more of the constraints are violated, the following loss function is minimized V = (u(t) - U u n c o n s t ( £ ) ) T (u(t) - X W o n s t W ) (3-51) subject to the above constraints. This is a quadratic programming (QP) [45] problem that can be solved for u(t). The standard form of the Q P problem is mm zT Qz + 2 fTz (3.52) z subject to Az < c (3.53) where z, f and c are column vectors, Q and A are matrices of appropriate dimensions. Dual Control of Paper Coating 54 The loss function V can directly be formulated as (3.52). The foregoing constraints are rewritten as where I " m a x - I " m a x I Au max 4" u(t - 1 ) - I u(t) < Au max tl(t - 1 ) D ^ l r a a x -D ^ l m a x M ^2max -M -^2max 1 - 1 0 0 0 1 - 1 0 ••• 0 D = 0 1 - 1 0 ••• 0 1 (3.54) (3.55) and M is given by (3.30). Efficient algorithms exist for solving Q P problems [15]. The solution procedure in-volves an iteration. To ensure that the procedure converges to a solution within a short time, the algorithm is interrupted if it exceeds a specified maximum number of iterations and the intermediate solution is used as an approximation to the optimal one. This constraint handling method performs better than the scaling method recom-mended by Braatz et al . [12] which fails when the interactions between adjacent actuators are not negligible. The proposed method has been found to work well in practice. 3.6 Mill Trials M i l l trials were performed in a Scottish m i l l on an off-machine bent blade coater producing an average of 15 tons of paper per hour. The coater had 4 coating stations, two of which did not have a C D control system. The other two had 35 actuators each, spaced at 100 m m intervals. The trials were held at the first coating station (figure 3.9) associated Dual Control of Paper Coating 55 Figure 3.9: Coating station where trials were performed. with the wire side of the sheet. The coating station was a Beloit S-Matic. Some technical information about this machine can be found in [44]. The sheet normally extended from the 4th to the 32nd actuator zone. Measurement points were mapped to the corresponding actuator zones using a spatial antialiasing filter followed by downsampling. The actuator setpoint range was from —700 to 700 microns. The blade was loaded using a solid backing bar. The average blade life was 9 hours and the scan time was 12 seconds. The coater had a traditional Dahl in controller with a so-called expert mode that turns off automatic control when it detects that the loop has become unstable. The coater was then allowed to run under manual control in stiff blade mode, and wide swings in process conditions and product quality were permitted to occur. In fact, at the time of the trials, automatic control was not operational for more than half of the blade life. 3.6 .1 Implementat ion Aspects The algorithm was coded as M A T L A B function M-files, translated into C + + source code using the M A T L A B Compiler, and compiled generating a stand-alone dynamic link library ( D L L ) for use on P C s running Microsoft Windows. A L a b V I E W application was developed to: call the algorithm D L L , provide the graphical user interface needed to interact with the algorithm in the mi l l , and communicate with the control processor Dual Control of Paper Coating 56 Experiment 1 zone 9 zone 18 zone 27 Dua l controller on-line gain estimates -0.001 -0 .002 -0 .003 Ga in estimates obtained v ia bump tests -0.001 -0 .002 -0 .002 Experiment 2 Dua l controller on-line gain estimates -0.0006 -0.002 -0.001 Ga in estimates obtained via bump tests -0.0004 -0.0012 -0.0012 Table 3.1: Ga in estimation (data courtesy of Honeywell-Measurex). retrieving various data and sending the actuator setpoints. The graphical user interface displayed al l relevant profiles and allowed the user to change the values of the algorithm parameters. Two important functions were implemented in the graphical user interface: a "switch" to enable/disable the controller and another switch to enable/disable sending the actuator setpoints. When the latter switch was off, the user could sti l l monitor the estimated gain profile and the generated actuator setpoints despite not sending them to the control processor. The implementation process took nearly four months. In a short pre-trial study, it was confirmed that the implemented algorithm can nor-mally communicate with the control processor via the L A N connection. It was also verified that the parameter estimator works very well. The verification was done by com-paring the estimated gains at some zones just before bump tests with those obtained from the bump tests using the Honeywell-Measurex system. The result was a significant matching (table 3.1). The process pole a, which could directly be obtained from the settings of the expo-nential smoothing filter, was set at —0.8. The noise filter zero c was identified to be —0.6 (see Appendix C ) . As described in [66], A was set such that A A T has a banded structure of the form 1 Q 1 this application, a suitable Q 1 or Q Q 1 . . . Qn-2 Q n - l . . . gn~3 Qn-2 . . . gn-4 gn-3 (3.56) [ g n ~ 2 £ n _ 3 • • • Q where Q is a scalar parameter in the range 0 < g < 1. For Dual Control of Paper Coating 5 7 . 2 I 1 _ J 1 1 1 1 -3 - 2 - 1 0 1 2 3 CD position Figure 3 . 1 0 : Identification of the actuator spatial response. value of Q lies in the region of 0 .2 . The average actuator spatial response (figure 3 .10 ) was identified using a series of bump tests. The response shape was five times as wide as the actuator zone size for this coating station. The spatial response was normalized so that its centre gain is unity and the actuator coupling matrix C was then constructed. A u m a x was set to 5 0 , 5 i m a x to 4 0 0 and ^ m a x to 6 0 0 . The controller design parameters were set as follows: A i = 0 . 0 0 0 1 , A 2 = 0 . 0 0 0 0 0 1 , A 3 = 0 and A 4 = 5 0 0 0 . These values were determined by scaling the parameters so that their corresponding tuning range is from zero to one, selecting reasonable values and finally fine-tuning. A3 = 0 meant that operation in stiff blade mode was allowed. Al though that was not favorable, this choice was imposed due to the inability to determine ^halfway ( m (3 -31) ) and in order to minimize production upsets. In essence, the concept that the operat-ing point can always be brought back into normal operation (by selecting appropriate actuator position) could not be verified for this coater. A t several times, al l the zones were persistently operating in stiff blade mode regardless of their actuator setpoints. The idealized nonlinearity of figure 3 .3 , though well studied in the literature, did not seem to hold (at the existing operating conditions) for this coater. Dual Control of Paper Coating 58 3.6.2 Experimental Results and Discussion A team consisting of Johan Backstrom from Honeywell-Measurex and the author of this thesis spent two weeks in the m i l l conducting experiments to assess the performance of the developed controller. The C D spread (2-o~), defined as double the standard deviation of the measurements across the sheet, is an important quality factor in the paper industry. For each tr ial , a plot of 2-o against time outlines the performance of the coating station under study. In a l l 3-D figures, the coat weight and actuator setpoints graphs display 35-zones blade-width profiles (off-sheet zones are included), while the C D coat weight and gain estimates graphs display 29-zones sheet-width profiles (from zone 4 to 32). Actuator setpoints refer to the setpoints generated by the dual controller. In figures 3.11-3.13, the Dahl in controller had already been turned off (by the expert mode) before the tr ial started. The coater was running under manual control producing a single-sided 80 gsm grade when the dual controller was switched on. The machine speed was 1000 m / m i n . The average 2-o was 0.3 gsm. The dual controller quickly achieved about 50% reduction in 2-o. The controller reacted to the streak at zone 22 by gradually moving the actuators while avoiding actuator picketing. The process gains were a l l positive confirming that the coater was operating in stiff blade mode. In the second tr ial (figures 3.14-3.17), the Dahl in controller was operating normally when the dual controller took over at scan 52. The grade was a double-sided 216 gsm coated paper and the machine speed was 570 m / m i n . No significant improvement could be obtained in this situation. The coater was operating in bent blade mode. A l l gains were negative except at the sheet near edge (zone 4) where the gain was almost zero. Consequently, the controller was ineffective at that location and the zone error could not be eliminated. The reason why the dual controller could not reduce the 2-o in this case is believed to be that the Dahl in controller was well-tuned and was yielding close-to optimal performance. To validate this statement, assume that a minimum variance controller were feasible, y(t) would then be a white noise sequence and 2-o would be about 0.19 gsm (because 0 was about 0.095 in figure 3.17). That is, the theoretical minimum value of 2-cr Dual Control of Paper Coating 59 was 0.19. The Dahl in controller was in fact providing an average 2-cr of about 0.26 gsm during this t r ia l (excluding the near edge error 5). When the Dahl in controller is in operation, the edge zones are often excluded from automatic control. That was the case in tr ial 3 (figures 3.18-3.21). The grade was a double-sided 170 gsm coated paper and the machine speed was 570 m / m i n . The dual controller was turned on at scan 28. Initially, no major improvement was observed but when A 4 was increased from 5000 to 10000 at scan 110, a significant drop in 2-tr was achieved. Al though probing increased, the output variance decreased due to improved estimates. This clearly illustrates the benefit of including probing in the control signal. More probing enabled the controller to accurately identify the gains at several C D posi-tions particularly at the far edge (zone 32). The error at that zone could then be regulated easily. Note that the gain at the other edge (zone 4) had reversed its sign. Tr ia l 4 (figures 3.22-3.25) lasted for an hour and a half. The results showed the coater drifting from bent blade to stiff blade operation. The grade was a single-sided 88 gsm coated paper and the machine speed was 1000 m / m i n . The dual controller was turned on at scan 57 after a manual start-up of the machine. The gain profile was mostly negative and the blade was operating in bent blade mode. A t that time, the performance of the dual controller was comparable to that of the Dahl in controller. A 4 was increased to 20000 at scan 170 but no improvement was noticed. It was reset back to 10000 when the dual controller was switched back on at scan 258. As the gain profile moved towards zero, the process became harder to control and the performance of the dual controller degraded considerably. Not before long, the gains were all positive and the coater was operating in stiff blade mode. To verify that, the Dahl in controller was switched on at scan 481. As it was expected, the loop became unstable and the Dahl in controller had to be switched off. The coater was then left running under manual control for an hour before tr ial 5 (figures 3.26-3.28) was authorized. The operator had manually reduced 2-cr to about 0.6 gsm before the dual controller was turned on at scan 800. The dual controller was able to further reduce it by more than 60%. The gains were al l positive and larger in magnitude than what had been observed in trial 4. 5 Whi le figure 3.14 indicates a value of 0.5-0.7 for the Dahlin controller 2-tr, one should note that the edge zones are included in all 2-cr plots. Dual Control of Paper Coating 60 The results of tr ial 6 (figures 3.29-3.31) demonstrate how the dual controller can recover from a major disruption. The grade was a double-sided 220 gsm coated paper and the machine speed was 550 m / m i n . The dual controller was turned on at scan 8 after another start-up. A sheet break occurred at scan 27. The coater had some of the zones operating in bent blade mode and the others operating in stiff blade mode. A t scan 75, a major disruption was about to enter the process. The machine operator was manually adjusting the front and back end stop actuators to eliminate the M D coat weight error. That had the effect of moving and t i l t ing the whole blade beam. Not aware of what the operator was doing at that time, we were facing an unexplained performance degradation to which we reacted by switching to the Dahl in controller at scan 84. As 2-a grew even larger, the dual controller was switched back on at scan 101. The coater was then operating in bent blade mode. As shown by the 2-a plot, the dual controller could smoothly recover from the disturbance. Figures 3.32-3.34 show the results of 180 scans of operation toggling between the Dahl in controller and the dual controller. Most of the zones were operating in stiff blade mode. 2-cr clearly deteriorated when the Dahl in controller was turned on at scan 47. The dual controller was switched back on at scan 85. A significant reduction in 2-cr could then be obtained. A t al l times, the dual controller offered a significant reduction in coat weight variations over the manual control policy currently adopted when operating in stiff blade mode. In bent blade mode, the dual controller did not always provide a performance improvement over the Dahl in controller. In fact, the dual controller was originally tuned to meet the normal operating performance of the Dahl in controller which was quite acceptable by the mi l l . The major benefit of the dual controller is that this satisfactory performance is attainable for most of the blade life, which clearly is not the case with the Dahl in controller. If improvement over the Dahl in controller is sought, more work on parameter tuning wi l l be required. Finally, in figures 3.35-3.41, the dual controller was allowed to run for a much longer period. The tr ial started with an old blade. The dual controller was turned on at scan 5. The coater was operating in stiff blade mode. A t scan 34, sudden sharp changes in Dual Control of Paper Coating 61 2-cr occurred. That was caused by the machine operator manually adjusting the M D coat weight. The controller was shut down at scan 425. Actuator positions were then reset at scan 447 and a new blade was installed at scan 475. The process was back on automatic control at scan 503. The gain profile became much smaller (in magnitude) and had different signs across the sheet. This contradicts the general belief that a new blade always starts operation in bent blade mode. Al though st i l l operating at the same grade, the same low 2-cr could no longer be maintained due to fading gains. A n important conclusion drawn from this tr ial is that the process is rather poorly set up. If operation in bent blade mode is truly desired, process settings (in particular, blade angle) should be adjusted properly. One more interesting thing to notice is that the dual controller started wi th the wrong sign for the ini t ia l gain estimate. The controller could, however, detect the correct sign of the gains and a valid gain estimate profile was obtained within a few scans. 3.7 C o n c l u s i o n s Prel iminary results of applying the proposed dual controller to a bent blade coater were generally successful. Applications of dual control to industrial processes is s t i l l in its infancy. The application discussed in this chapter is certainly a novelty. The developed controller was well received by the company. Encouraged by the results, the mi l l crew asked to continue using it - on an experimental basis- after the tr ial period. The controller was operational in the m i l l for about two months before it was suspended due to lack of technical support and tuning guidelines. Nevertheless, moving this controller into a commercial product is currently underway. Several contributions of this work are now outlined: • this is believed to be the first application of dual control to a coupled multivariable industrial process, • the potential usefulness of including probing in the control signal has been demon-strated, Dual Control of Paper Coating 62 • the use of an adaptive Ka lman filter to track gain variations proved to be worthwhile, and • the controller yields substantial quality improvement as it is able to regulate the process throughout the blade life. Dual Control of Paper Coating Figure 3.11: Tr ia l 1, the dual controller provides about 50% reduction in C D coat wei variations. Dual Control of Paper Coating 64 20 > Figure 3.12: 3-D plots of trial 1. Dual Control of Paper Coating 65 Figure 3 . 1 3 : 3 - D plots of trial 1 (continued). Dual Control of Paper Coating 66 1.5 1 1 0.5 -Hah Lin Dual c ontroller 1 nqhlin controller controller A . .A «A j \r I 50 100 150 200 250 300 350 400 450 500 Number of scans x 10 CD -10 150 .2 100 o 5 50 c o t o I -50 < -100 -150 scan 350 J sc an 250 //\ ^ . % sc£n 150 10 15 20 CD position 25 30 10 15 20 CD position 25 30 35 A s can 350 \ \ scan 25 / • \ V 0 \ . // . . \ N \ \ \ Y ^z^y—~\ /• J i e f , a n - T ^ f l - Y N o c / ' v ' . . . / . . . \ -35 1 \ ,\ : i scan 350 / ; 1 v \ • ^ :\ [:::::' ' " :scan 150 i scan 250 i 1 0) B 1 o o Q -1 O 10 15 20 CD position 25 30 Figure 3.14: Trial 2, the performance of the dual controller is comparable to the Dahlin controller. Dual Control of Paper Coating 67 Figure 3.15: 3-D plots of trial 2. Dual Control of Paper Coating 68 150 > Figure 3.16: 3-D plots of trial 2 (continued). Dual Control of Paper Coating 6 9 100 150 200 250 300 Number of scans I L L i 350 400 100 150 200 250 300 Number of scans 350 400 Figure 3.17: Noise variance estimation for trial 2. Dual Control of Paper Coating 1.5 1 1 0.5 I Duai controller i i Dahlin controller probing increased i i 50 100 150 200 Number of scans 250 300 350 2 • i 1 o Q. s  0 ra E CD 1 c 'CO CD_2 -3 0 400 2 300 2 Z7 200 _c £ 100 to 1 o < -100 -200 x 10 ' : 1 ' : scan 100 K / \ ^ / \ / / V Scan 200 \ / scan 300 10 15 20 C D position 25 30 35 I I I I i scan 300 scan 100 : / ' \ • • / v. 4 * scan 200 I I I I i 10 15 20 C D position 25 30 35 | o Q -1 O I i scan 200 // V / — ' "7 scan 100 i i i / scan 300 10 15 20 C D position 25 30 Figure 3.18: Trial 3, the benefit of probing is illustrated. Dual Control of Paper Coating 71 20> Figure 3.19: 3-D plots of trial 3. Dual Control of Paper Coating 72 Figure 3.20: 3-D plots of tr ial 3 (continued). Dual Control of Paper Coating 73 Figure 3.21: Noise variance estimation for trial 3. Dual Control of Paper Coating 74 manual 1 1 . .rii ial .pnntrnNpr - • cbntrol • manual, control Dahlin controller controller i 1.5 1 1 0.5 50 100 150 200 250 300 Number of scans 350 400 450 500 x 10 ' =5 2 E 0 CO o scan 495 I - V x V / y , < — \ \ i / scan 200 • : scan 400 i i I 0 400 2 200 Q. a-200 < -400 10 15 20 CD position 25 30 10 15 20 CD position 25 30 "S 1 '(B 0 Q -1 O 35 SC£ I in 200 \ : s / / / / / / \ \ \ \ \ \ M / / \ — y \ _ 1 \ • _ ' scan 400 i 35 - I 1 — - . _ scan 495 i / ^ / / ~~ ' • -~ ' " : V " C . ::-;:::.rv.vr ' > - ^ \ — scan 200 i : scan 400 i 10 15 20 CD position 25 30 Figure 3.22: Trial 4 , performance degrades as the coater drifts into stiff blade operation. Dual Control of Paper Coating 75 Figure 3.23: 3-D plots of trial 4. Dual Control of Paper Coating 76 Figure 3.24: 3-D plots of trial 4 (continued). Dual Control of Paper Coating 77 Number of scans Number of scans Figure 3.25: Noise variance estimation for trial 4. 0.81 1 1 1 1 r 800 850 900 950 1000 1050 1100 1150 Number of scans x 10"3 o Q. £ 3 a E CO 0 1 scan 1100 I i \ \ . . Y • • • \ • • -x- • • • — \ / • — • f. \ .x./.„ < -\ '• ""7 N - — • / v ' Y -scan 900 i scan 1000 i I 0 5 10 15 20 25 30 35 CD position Figure 3.26: Trial 5, the dual controller reduces coat weight variations in stiff blade mode. Dual Control of Paper Coating 78 Dual Control of Paper Coating 79 Figure 3.28: 3-D plots of trial 5 (continued). Dual Control of Paper Coating Dahlin 50 100 150 200 Number of scans 250 300 x 10 ' 15 20 CD position i i i / \ s \ scan 300 \ / \ _ ' ...Jl 1 ^ ^. scan 70 A . . . / . j scan 120 i i • i 2 •5 Q- 1 '» 0 o _1 D ° - 2 -3 10 15 20 CD position 25 30 Figure 3.29: Trial 6, the dual controller recovers from a major disruption. Dual Control of Paper Coating 81 Figure 3.30: 3-D plots of trial 6. Dual Control of Paper Coating 82 Figure 3.31: 3-D plots of trial 6 (continued). Dual Control of Paper Coating 1.5 manual Dual Dahlin manual 0.5 control controller controller control 1 ' . ' ' Dual controller sheet break 20 40 60 80 100 120 Number of scans 140 160 180 4 JB 3 1 2 CD 1 1 x 10 0 -1 -2 0 300 I : scan 160 i -v\ ; U V \ Y \ J — scan 120 \\ '/ V/ V/ i scan 45 10 15 20 CD position 25 30 35 2 200 3 100 h < -200 -300 scan 120 ^ i / ' A scan 160 / "A • ^ > scan 45 vV / / V V V / / \ \ \ fi \ A V/ > \ v r ^ - J \ V . . . \ . . ' A ' / V 10 15 20 CD position 25 30 35 1.5 I 1 o CL Z 0.5 15 0 o CJ D O -0.5 i i scan 45: / ^ — j - . _ / C -/ / • x / / / scan 120 i \ scan 1 60 i 10 15 20 CD position 25 30 Figure 3.32: Trial 7, the dual controller outperforms the Dahlin controller. Dual Control of Paper Coating 84 Figure 3.33: 3-D plots of trial 7. Dual Control of Paper Coating 85 Figure 3.34: 3-D plots of trial 7 (continued). Dual Control of Paper Coating 2 1.5 0.5 0 \ : i b l a d e sheet change : break • ; . VW*N^ .... V 1 1 1 1 1 0 200 400 600 800 1000 1200 Number of scans x10~3 -2 I sc i an 300 1 scan 1000 - - ^ . . . . . T - . r - X „ . _ . ' < ' ' ' • ^ K — ' ' ~-i 1 ^/ \ : " • scan 800 i 0 5 10 15 20 25 30 35 CD position 400 a> 300 I 200 § 100 a. aS cn o o | -100 o -200 -300 i i i i - - - cr*-ar* .1.nfifV / -N ) scan n uuu: / \ S p a n soo : ' s ^ / /. scan 300 : / , : \ \ • . - , \ . / ^ • v v _/_/ ^___x / \ y " i i i i i 0 5 10 15 20 25 30 35 CD position 5 10 15 20 25 30 CD position Figure 3.35: Trial 8, the dual controller is running for about 4 hours. Dual Control of Paper Coating 87 Figure 3.36: 3-D plots of trial 8, scan number 1-600. Dual Control of Paper Coating 88 Figure 3.37: 3-D plots of trial 8, scan number 1-600 (continued). Dual Control of Paper Coating 89 Figure 3.38: Noise variance estimation for trial 8, scan number 1-600. 600 700 800 Numl iber of 1000 scans 0.1 oP-08 0.04 1100 120C 600 700 800 900 . 1000 Number of scans 1100 120C Figure 3.39: Noise variance estimation for trial 8, scan number 600-1200. Dual Control of Paper Coating 90 Figure 3.40: 3-D plots of trial 8, scan number 600-1200. Dual Control of Paper Coating Figure 3.41: 3-D plots of trial 8, scan number 600-1200 (continued). Chapter 4 A New Approach to Adaptive Kalman Filtering Classical adaptive controllers are often criticized for their start-up performance. This is due to the fact that certainty equivalence adaptive controllers -not paying attention to the commonly high uncertainty at start-up- boldly attempt to control the process straight off. As discussed in Chapter 1, dual control provides a subtle solution to this problem. The dual control strategy implemented in Chapters 2 and 3 utilizes an adaptive K a l m a n filter to track parameter variations. The adaptive Ka lman filter scheme adopted, however, has the disadvantage that one has to wait for some time before starting to use the noise variance estimates in the Ka lman filter. This means that at start-up, the K a l m a n filter is relying on the fixed ini t ia l estimates of the noise variance. This clearly contradicts the foregoing advantage of dual control. The need for a better adaptive K a l m a n filtering scheme was perceived during the m i l l trials of Chapter 3. The development of this scheme is the subject of this chapter. 4.1 Introduction The Kalman-Bucy formulation of the filtering problem assumes complete a priori knowl-edge of the process and measurement noise variances. Clearly, this is not a realistic assumption. Using wrong statistics in the design of a K a l m a n filter can lead to large estimation errors or even to divergence of the state estimates [39]. There exists a number of different approaches to adaptive K a l m a n filtering [48]. Mehra [47] classifies the various approaches into four groups: Bayesian estimation, max-imum likelihood, correlation methods and variance matching methods. Our aim is to use adaptive filtering as a tool for identification of time-varying systems. Most of the methods reported in the literature are for various reasons not applicable to our problem. 92 A New Approach to Adaptive Kalman Filtering 93 We need a method that recursively can obtain estimates of the noise variances when the state-space system is time-varying. The Expectat ion-Maximizat ion ( E M ) algorithm formulated by Dempster et al . [16] is a modification of the maximum likelihood ( M L ) method. It is a two step iterative batch scheme for maximization of the likelihood function avoiding the numerical optimization required in many M L problems. A batch E M scheme is first developed to determine the estimates of the noise variances. This scheme is related to the work of Tanaka and Katayama [57], and is based on a Ka lman smoothing reconstruction performed using fixed-interval smoothing formulas. The estimates are obtained via an iterative procedure. This , together with the fact that the entire set of observations is used, makes it a challenging problem to calculate the estimates on-line. This problem is addressed and a recursive scheme based on the batch one is derived. This recursive scheme is the E M adaptive K a l m a n filter. 4.2 Problem Statement Consider a linear, discrete, time-varying dynamic process described by the state-space model x(t+l)=F{t)x(t) + w(t) (4.1) y(t) = H{t)x(t)+e(t) (4.2) where x is an n-dimensional state vector, y is an r-dimensional measurement vector, with the matrices having appropriate dimensions. We assume w(t) and e(t) to be mutually independent sequences of zero mean Gaussian white noise wi th variance matrices E{w(t)wr{t)} = Q (4.3) E{e{t)eT{t)} = R (4.4) in which Q and R are nonnegative definite matrices that are assumed unknown. The adaptive Ka lman filtering problem consists of obtaining on-line estimates of x(t) given the observation set Yt = {y(0),..., y(t)}. Let x(i\j) be an estimate of x(i) based A New Approach to Adaptive Kalman Filtering 94 on the observation set Yj and P(i\j) = E{(x(i) - x(i\j))(x(i) - x(i\j))T} (4.5) The solution to this problem when Q and R are exactly known is given by a K a l m a n filter of the form e(t\t-l) = y(t)-H(t)x(t\t-l) (4.6) K(i) = P(t\t - l)HT{t) (H(t)P(t\t - l)HT{t) + R) _ 1 (4.7) x{t\t)=x(t\t-l)+K(t)e{t\t-l) (4.8) P(t\t) = (l-K(t)H(t))P(t\t-l) (4.9) x{t+ l\t) = F(t)x(t\t) (4.10) P{t + l\t)= F(t)P(t\t)FT{t) + Q (4.11) wi th the in i t ia l conditions x0 = x(0\ — 1) and P0 — P(0\ — 1). In the following sections we derive an E M method to estimate Q and R. 4.3 The E M Algorithm The E M algorithm is related to the Max imum likelihood ( M L ) method. M L is a principal method to solve an important problem in system analysis, that is to estimate values for the unknown parameters of a system given measurements which are some functions of these parameters. Let Y denote a measurement vector with the associated probability density function /y(y|#), conditioned on the (unknown) parameter vector 9. /y(y|#), commonly called the likelihood function, is a function of y and 9. If 9 were known, the most likely value of Y can be computed by maximizing /y(y|#) wi th respect to y . Given an observation y* of Y, the M L estimate, #ML> is the value of 9 that maximizes the log-likelihood function, that is, ^ML = arg max log fY (y* \9) (4.12) A New Approach to Adaptive Kalman Filtering 95 The problem is that the maximization in (4.12) is often complicated. The procedure and difficulties of the M L method are best illustrated by the following examples [60]. Example 1 Assume that we have the data set = {yi, • • • ,VN} where y are independent Gaus-sian random variables N ( / U , < J ) , then Y^ has a multivariate normal distribution with the probability density function (pdf) fvMo) = n (7^ e _ (" ) 2 / 2 < Tj <4-13) where 9 = {p, a}. We need to find a value for 9 given a set of observations y * of Y^. The log-likelihood function is given by log/yw(y|0) = -TVlogx/27- iVloga - f: ( y * J 0 ' (4.14) fc=i Differentiation yields ^ l o g / , „ ( y | 0 ) = g ^ and, d l o g / y w W ) . = - ^ + E ^ # (4.16) da ' a ^ aA Equating this to zero gives and, £ = i $ > (4-17) k=l a = 1 N M ^ E ^ - A ) 2 (4-18) which is a well-known result. Notice that, in this case, we get a closed form solution. Example 2 Assume that we have the data set YN = {y\,..., yN} where y^ wi th the probability p is N (^ i , tT i ) and with the probability (1 — p) is N(^2,er 2), i.e., the pdf of yu is fyk(y\0) =p / i (y ,Mi , cr i ) + (1 -p)h(y,p*,o-2) (4.19) A New Approach to Adaptive Kalman Filtering 96 where 9 = {pi,&i, p2,o~2,p}, fi and f2 are the corresponding normal pdf's. The log-likelihood function is log/yw(y|0) = £ log ( p ^ L - e - ^ W + ( i _ p ) L - e - ^ - ^ 2 ^ ) (4.20) £Ti V v2iro1 V2TTO2 J This time logfYN{y\9) can not be maximized analytically and one has to resort to nu-merical optimization. The E M algorithm is a two step iterative scheme for solving M L problems, avoiding the numerical optimization required in many situations. It was first introduced in [16] for the case of missing data. The basis for the E M algorithm is the notion of incomplete and complete data. Incomplete data is the data actually being observed. One then has to find some convenient completion of the data, leading to a simpler likelihood function given that the complete data was observed. Each iteration of the algorithm involves two steps called the Expectation step (E-step) and the Maximizat ion step (M-step). A brief review of the mathematical background to the E M algorithm is given below. For an exhaustive presentation of the E M algorithm and its application to several signal processing problems, see [27]. Suppose that the data vector Y can be viewed as being incomplete, and that there is a way to extend Y wi th some non-observed "fictitious" data into an extended data set Z related to Y by T(Z) = Y (4.21) where T(-) is a noninvertible (many to one) transformation such that i f Z is observed, an observation of Y is available too, but not vice versa. Z w i l l be referred to as the complete data set. The probability density function of the complete data, denoted fz{z\9), is also conditioned by the parameter vector 9. The idea is that the complete data Z is chosen such that the solution of <W = argmaxlog/ z(z*| cr) (4.22) is straightforward, typically in closed form. The E M algorithm, presented below, wi l l use the simple procedure for M L estimation in the complete data model (4.22), as a part of an iterative algorithm to solve the M L problem in the observation model (4.12). A New Approach to Adaptive Kalman Filtering 97 The E M algorithm • Start, n = 0 : guess 6>(0) • Iterate (until convergence) - The E-step: calculate T(P\eW) = E { l o g / z ( z | 0 ) \Y = y*, 0™ } (4.23) - The M-step: solve 0(^+1) = a r g m a x ^ l ^ ) (4.24) The iterative algorithm consists of two steps. In the E-step the sufficient statistics of the complete data are estimated, i.e., the expected value of l og / z (z | 0 ) conditioned on Y is computed. This estimate is then used to obtain a new estimate of 6 by maximizing the log-likelihood of the complete data in the so-called M-step. This parameter estimate is then fed back to the E-step, and so the method iterates until convergence. The iteration is justified intuitively as follows. We would like to choose 6 that maxi-mizes log/z(z*|(9), the log-likelihood of the complete data. However, since log fz(z*\0) is not available to us (because the complete data is not available), we maximize instead its expectation, given the observed data y*, and the current value of the parameters 0(n\ The choice of the complete data may critically affect the complexity and the conver-gence properties of the algorithm. A n unfortunate choice of the complete data w i l l yield a completely useless algorithm. The algorithm, on the other hand, converges rapidly if the complete data is chosen such that it can be predicted well given the observations [11, 71]. Thus, it takes creativity to apply the E M algorithm to a given problem. Application of the E M algorithm to example 2 Define WN = {wi,..., w^} where r wk = S (4.25) 0 if yk e h A New Approach to Adaptive Kalman Filtering 98 if yk and wk were observed at al l times, the M L problem would have been straightforward. The complete data set is therefore chosen as ZN = {YN,WN} (4.26) The likelihood function is given by fZN(z\8) = /Z w(y,w|0) = /yN(y|w,c9)/^(w|c9) (4.27) But , /V„(y|w,0)= n ( - ^ e - l y - r t ' M ) n ( - ^ e - ^ - " ^ / ^ l (4.28) and Wk is a Bernoull i random variable which is 1 with probability p and 0 wi th probability 1 — p, then fwN(^)=UPWk(1-P)1'Wk > «>* = 0,1 (4-29) k=l In the E-step of the E M algorithm we compute the expected value of the log-likelihood function, that is, Ho\o [n)) = E ( - A n o g v ^ - ( f > W i - E [ v k ~ ^ ? Vrr2 \k=l J k,(wk=l) 1 - f i V - ^ % j logf f j " E V k=l I k,(uik = (Vk - P2)' (Wk=0) 2 C J2 + (j2 ™*) !ogP + ( E 1 - ™ * ) l o § ( l - P)\YN = Y*, 0{n) } (4-30) This gives T{o\eW) = - T v i o g v ^ - ( E ^ ) l Q g a i - E ( ? / f c o 2iy™k \k=i J k=i 2<7i - NV - E ™* logcr2 - 2^  — — 2 — ( i - wk) \ k=l / k=l za2 + ( E l o g p + ^ E 1 - ^ ) l o g ( ! - P) ( 4 - 3 1 ) A New Approach to Adaptive Kalman Filtering 99 where wk = E{wk\y*k,eW} P ( B ) / i ( y U i B V r ) (n) (n)s k = l,...,N (4.32) ^•((9|^ ( n )) has then to be maximized with respect to 9 in the M-step. Differentiating (4.31) and equating the derivatives to zero gives E N k=lWk Ol Efc=i wk(yk - Ai)2 \l is* Ek=i(l-™k)yk E f = 1 ( i - ^ ) 0-2 Ef = 1 ( l-^)(y f c -A 2 )2 TV *;=i The iterative scheme can then be encapsulated as follows • Start, n = 0 : guess 6>(0) • Iterate (until convergence): - Compute wk (4.32) - Compute 0 ( n + 1 ) ((4.33)-(4.37)) - Increment n (4.33) (4.34) (4.35) (4.36) (4.37) The convergence properties of the E M algorithm are discussed in [11, 71]. If there exists a unique global maximum of log/z(z* |#) , and no other stationary points, then the algorithm converges to this global maximum. If not, however, there is no guarantee that the estimate converges to the global maximum. The rate of convergence of the E M A New Approach to Adaptive Kalman Filtering 100 algorithm in general depends on the amount of incomplete data in relation to complete data, see [27]. The more non-observed data the algorithm has the more iterations that have to be performed to reach convergence. 4.4 E M Estimation of Q and R Recall from section 4.2 that we have the model x(t + 1) = F(t)x(t) +w(t) (4.38) y(t) = H(t)x(t) + e(t) (4.39) If x(t) and y(t) were observed at all time instants, the problem of determining Q and R would have been straightforward. Therefore, the choice of the complete data set ZN for this problem seems obvious. Let ZN = {YN,XN} where YN = {y(0), • • • ,y(N)} and XN = {x(0),... ,x(N)}. The probability density function of ZN is given by fzA*\°) = /zw(y,x|0) = /y„(y|x,0)/* w (x|0) (4.40) where 9 is {Q,R}. Using the fact that given x(t), y(t) is Gaussian with mean H(t)x(t) and variance R, then fy ( y | x = TT | 1 c-h_(y(k)-H(k)x(k))TR-Uy(k)-H(k)x(k)) j ^ ^ Similarly, N ( i fv C x i m = FT _ _ | ^ _ p - | ( x ( A ! ) - F ( f e - l ) x ( f c - l ) ) T Q - l ( i ( f c ) - F ( f c - l ) x ( A : - l ) ) I e - i (s;(0)- i o ) T p- 1 (a;(0)- i o ) (4.42) A New Approach to Adaptive Kalman Filtering 101 App ly ing the E M algorithm given an observation y* of YN : the E-step H0\0{n)) = E {log fyN (y|x, 0) + log fXN (x|0) \YN = y*, } (4.43) Using (4.41) and (4.42) H0\0(n)) = E J - ^ ^ l o g 2 7 r - ^ ± ^ l o g | J R | 1 N N N - ^ !>(*) - H^xWfR-'iyik) - H(k)x(k)) - - log27r - - log\Q\ k=0 N Y,Wk) - F i k ~ l X k ~ l)fQ~l{<k) - F(k - l)x{k - 1)) YN = y\e^} (4.44) fc=i l- log 27T - l- log\P 0\ - \ (x(0) - Xof P-1 (x(0) - x0) Simplifying (4.44) using the matrix property: zTCz = tr (zTCz} = tr (CzzT^J where z is a column vector, C is a square matrix of appropriate dimensions, and tr denotes the trace of a matrix, we get F{6\e^) = - - £ > [it!"1 (e(A; |N)e T ^|N) + H(k)P(k\N)HT(k)) 2 fc=0 1 ^ - - E t r [Q'1 {Mk\N) + F{k - l)A(k - l\N)Fr{k - 1) 2 fc=i - B(k\N)FT(k - 1) - F(k - 1)B T(A;|A^))] + O (4.45) where 0 includes the terms independent of Q and R, and e(fc|jV) = y(k) - H(k)x(k\N) (4.46) A(Jfc|iV) = E {x{k)xT(k) \YN } = x(k\N)xT(k\N) + P(k\N) (4.47) B(k\N) = E {x(A;)x T (^ - 1) \YN } = £(A;|A^)xT(A: - 1|7V) + P(k, k - 1\N) (4.48) A New Approach to Adaptive Kalman Filtering 102 wi th P(k, k — 1\N) = E {{x(k) - x(k\N))(x{k - 1) - x(k - 1\N))T} (4.49) In (4.46)-(4.49), x(k\N), P(k\N) and P(k,k - 1\N) are the state estimate and the error covariance matrices given al l the measurements YN. This implies the use of some kind of fixed-interval smoothing. W i t h al l observations stored, the fixed-interval smoother works - i n a batch procedure- by running a Ka lman filter forward in time followed by the backward algorithm [3]: for k = N — 1, . . . , 0 S(k) = P{k\k)Fr(k)p-1(k + l\k) (4.50) x(k\N) = x(k\k) + S(k) (x{k + 1\N) - x{k + l|fc)) (4.51) P(k\N) = P{k\k) + S{k) (P{k + 1\N) - P{k + l\k)) ST(k) (4.52) and for k = N — 1 , . . . , 1 P(k, k - 1\N) = P{k\k)ST(k - 1 ) + S{k) [P{k + 1, k\N) - F(k)P{k\k)) ST(k - 1) (4.53) with the boundary value P(N, N - 1\N) = (I - K(N)H(N)) F(N - 1)P(N — l | iV — 1) (4.54) In the M-step of the E M algorithm we have to maximize T{9\9^) wi th respect to Q and R. Differentiating (4.45) [50] and setting the derivatives to zero gives 1 ^ ( n + 1 ) = ^ ^ E ^ | i V ) (4-55) and where Q { n + 1 ) = ijEL(k\N) (4-56) V(k\N) = e(k\N)eT(k\N) + H(k)P{k\N)HT(k) (4.57) 1 Making use of the fact that Q and R are both symmetric. A New Approach to Adaptive Kalman Filtering 103 and L{k\N) A(k\N) + F(k - l)A(k - l\N)FT(k - 1) - B(k\N)FT(k - 1) - F(k - l)BT(k\N) (4.58) The batch E M method can then be summarized as follows • Start, n = 0 : guess Q ( 0 ) and i ? ( 0 ) • Iterate (until convergence): - R u n the fixed-interval smoother to obtain x(k\N), P(k\N) and e(k\N) for k = 0 , . . . , N and P(k, k - 1\N) for k = 1 , . . . , N - Compute V(k\N) for k = 0 , . . . , TV (4.57) - Compute A(k\N), B(k\N), L(k\N) for k = 1 , . . . , N ((4.47), (4.48) and (4.58)) - Determine R^n+1^ and Q^n+l^ ((4.55) and (4.56)) ' - Increment n 4.5 A Recursive E M Update To obtain a recursive formula for the update of the EM-solut ion, one should try to recursify the iteration steps ((4.55) and (4.56)). Assume that, at time t — 1, the estimates are given by R(t-l) = 7 E V(k\t - l ) (4.59) and 1 t-i Q[t~1] ( t - l ) (4.60) Then, at time t (4.61) and (4.62) A New Approach to Adaptive Kalman Filtering 104 A simple way to recursify (4.61) and (4.62) would be R[t) « Tj^Yy + V(t\t)] (4.63) and gW *j[(t- l)Q{t~l) + L(t\t)] (4.64) B u t since the expectations in (4.59) and (4.60) are only conditioned on measurements unti l time t — 1, not t, this approximation often causes the simple recursive scheme of (4.63) and (4.64) not to converge. We need to express V(k\t) and L(k\t) in terms of V(k\t — 1) and L(k\t — 1). To do so, the following fixed-point smoother wi l l be utilized [3]. For some point in time k, and t = k + 1, k + 2,... Pa(k, t) = Pa(k, t - 1) (F(t - 1) - F(t - l)K(t - l)H(t - 1)) T . (4.65) Ka{k,t) = Pa(k,t)HT{t) (R® + H(t)P(t\t- l ) H T ( t ) ) - 1 (4.66) x(k\t) = x(k\t - 1) + Ka(k, t)e(t\t - 1) (4.67) P{k\t) = P(k\t - 1) - Pa(k, t)HT{t)Kj(k, t) (4.68) P(k,k-l\t) = P(k,k-l\t-l)-Pa(k,t)HT(t)Kj(k-l,t) (4.69) with the init ial ization conditions: at t = k Pa(k,k) = P(k\k-l) (4.70) and P{k, k - l\k) = (I - K(k)H(k)) F(k - l)P(k - l\k - 1) (4.71) A proof of (4.69), which does not belong to the standard fixed-point smoothing formulas, is provided in Appendix D . Substituting (4.67)-(4.69) into (4.46)-(4.48) yields e{k\t) = e(k\t - 1) - H(k)Ka{k: t)e{t\t - 1) (4.72) A New Approach to Adaptive Kalman Filtering 105 A(k\t) = A(k\t-1)+ W{k\t) (4.73) B{k\t) = B[k\t-l)+0{k\t) (4.74) where W{k\t) = Ka(k, t)e(t\t - l)eT(t\t - l)Kj(k, t) + Ka(k, t)e(t\t - l)xT(k\t - 1) + x(k\t - l)eT(t\t - l)Kj(k, t) - Pa(k, t)HT(t)Kj(k, t) (4.75) and 0{k\t) = Ka(k,t)e(t\t - l)eT{t\t - l)Kj{k - l,t) + Ka(k,t)e(t\t - l)xT(k - l\t - 1) + x{k\t - l)eT(t\t - l)Kj(k - 1, t) - Pa(k, t)HT(t)Kj(k - 1, t) (4.76) Then, using (4.57) and (4.58), we get V(k\t) = V(k\t-l) + N{k\t) (4.77) and L(k\t) = L(k\t - 1) + M(k\t) (4.78) where N(k\t) = H(k)Ka(k, t)e{t\t - l)eT{t\t - l)Kj{k, t)HT{k) - e(k\t - l)eT(t\t - l)Kj(k, t)Hr(k) - H{k)Ka{k, t)e(t\t - l)eT(k\t - 1) - H(k)Pa(k, t)Hr(t)Kj(k, t)Hr(k) (4.79) and M(k\t) = W{k\t) + F(k-l)W(k-l\t)FT{k-l) -0{k\t)FT(k - 1) - F(k - l)0T(k\t) (4.80) A New Approach to Adaptive Kalman Filtering 106 Finally, the estimates are given by and Q<*> + t - i t J R ( * - 1 ) + y ( t | t ) + E ' A r ( A : | t ) fc=0 t - i (t-l)Q^ + L(t\t) + ^M(k\t) k=l (4.81) (4.82) To have a truly recursive scheme, the summations in (4.81) and (4.82) can not go over the entire set of observations. Fortunately, for some k, N(k\t) and M(k\t) vanish as t increases. This is due to the dependence of N(k\t) and M(k\t) on Ka(k,t) and the fact that all the smoothing that is possible to achieve can be achieved in two or three times the dominant time constant of the Kalman filter [3], i.e., Ka(k,t) —>• 0 as t increases. Consequently, (4.81) and (4.82) can be rewritten as (t + l) t - i tBSt-V + V(t\t)+ £ N(k\t) and t k=t-T t-1 (4.83) (4.84) (t-l)Q^+L(t\t)+ £ M(k\t)\ k=t-T J where T is a design parameter. T can safely be set to three times the maximum expected value of the Kalman filter time constant. Simulations show that (4.83) and (4.84) often converge to biased estimates. This can be explained as follows. First, notice that the dependence of V and L on R and Q is omitted in the derivation of the scheme, for instance, V(k\t — 1) and V(k\t) in (4.59) and (4.61), should read V(k\t - l,R^) and V(k\t, R ^ ) . However, in the recursive update formulas (4.77) and (4.78), it is implicitly assumed that V(k\t), V(k\t— 1), L(k\t) and L(k\t — 1) are computed using the same values of R and Q. Then, in (4.83) and (4.84), this assumption is exploited and El=o V(k\t—l) and £fc=i L(k\t—1) are replaced by tR^~1^ and (t — 1)Q^~^ with the result that (4.83) and (4.84) will converge to biased estimates especially if the initial estimates are far from the true values. To fix this problem, a recursive least squares (RLS) estimator with a so-called forgetting factor A [6] will be A New Approach to Adaptive Kalman Filtering 107 used. Equation (4.84) then becomes KQ(t) = PQ(t) (AI + PQW)- 1 col ( Q « ) =col (QW) + KQ(t) t-i col i V(t\t) + £ N(k\t) - col PQ(t + l) = (l-KQ(t))PQ(t)/\ Similar update formulas for are easy to obtain. (4.85) (4.86) (4.87) To reduce the computational effort, one may take advantage of the symmetric nature of Q and R and only estimate their unique elements. If Q is n x n , it has (n2 + n)/2 unique elements. We therefore replace the n2-vector col by applying the col operator only to the lower triangular half of the matrix. 4.6 Simulations Assume that we have the following model for y(t) y(t) = b(t)u(t - 1) + e(t) (4.88) where b(t) can be modeled as b(t + 1) = b(t) + w(t) (4.89) This time-varying system is then simulated with E {w2(t)} = q = 0.5xlO - 5 and E {e2(t)} — r = 0.01 and the input u(t) as a square wave. In Figure 4.1, we see the results of the application of the batch EM variance estimation method. On-line estimation of q and r is shown in figure 4.2 where it is compared to Isaksson's method [38]. Isaksson's method requires that one waits for at least 75 samples before updating the estimates and for another 75 samples before using the estimates in the Kalman filter. Such a requirement is not necessary for the EM scheme. The initial-estimates were r^ = 0.16 and g(°) = 0.0001. A New Approach to Adaptive Kalman Filtering 108 Consider the following time-varying process as another example y(t) + a(t)y(t - 1) = b{t)u(t - 1) + e{t) where a(t) and b(t) can be modeled as (4.90) " Kt) ' " Kt -1) " + = a(t - 1) _ w2(t) _ 0.01 0 (4.91) The results of the with the true parameter values r = 0.01 and Q 0 0.0001 batch E M scheme are shown in figures 4.3 and 4.4. Figures 4.5 and 4.6 show that the straight averaging scheme ((4.83) and (4.84)) converges to biased estimates. The results of the recursive EM scheme are shown in figures 4.7 and 4.8. For both examples, the parameters of the recursive estimation were A = 0.96 and T = 50 samples. Based on the simulations presented, we conclude that the EM scheme provides a smoother transient than Isaksson's method. Although both methods appear to converge in about the same number of samples, the EM scheme has the advantage that one does not have to wait for some time before starting to use the noise variance estimates in the Kalman filter. Such an advantage would prove useful in situations where the parameter estimator is often interrupted, as is the case in paper coating (Chapter 3) where the controller is shut down for sheet breaks, blade changes and grade changes. Checking back the results of trial 8 in Chapter 3, one realizes that gain estimation can be boosted by employing the EM scheme developed here as it eliminates the "wait" period each time the controller is restarted. Finally, in figures 4.9 and 4.10, the recursive EM scheme is reapplied to the process of T 0.01 0.001 example 2 (4.90) but in this case Q is non-diagonal and is given by Q = [ 0.001 0.0001 This is a more interesting case as the noise sequences W\ and w2 are now correlated. As seen in the figure, the application is successful and the EM scheme could get the directions right. A New Approach to Adaptive Kalman Filtering Figure 4.1: Batch estimation of q and r for example 1. A New Approach to Adaptive Kalman Filtering 110 0.1 I 0.05 0 -0.05 I /. ^ ^ -s. ^ j i 0 100 200 300 400 500 600 700 800 900 1000 Number of samples I i I I i i I i 1 1 1 0 100 200 300 400 500 600 700 800 900 1000 Numberof samples Figure 4.2: On-line estimation of q and r (solid) and Isaksson's method estimates (dashed line) for example 1. A New Approach to Adaptive Kalman Filtering 111 Figure 4.3: Batch estimation of Q and r for example 2. A New Approach to Adaptive Kalman Filtering 112 Figure 4.4: Batch estimation of Q and r for example 2 (continued). A New Approach to Adaptive Kalman Filtering 113 Figure 4.5: On-line estimation of Q and r using the straight averaging scheme ((4.83) and (4.84)) for example 2. A New Approach to Adaptive Kalman Filtering 114 Figure 4.6: On-line estimation of Q and r using the straight averaging scheme ((4.83) and (4.84)) for example 2 (continued). A New Approach to Adaptive Kalman Filtering 115 0.4 0.35 0.3 0.25 0.2 •I 0.15 0.1 0.05 Oh - 0 . 05 -0.1 ! I " I I I I I I I I I l ! \ ! - \ - 1 . - ^ i i I-•i _ i-, •i i i 0.12 0.1 0.08 a> 0.06 to E o- 0.04 0.02 0 h -0 .02 0 100 200 300 400 500 600 700 800 900 1000 Numberof samples III III lil III. Ill III I "' I lil V " \ " . \ ii \ i i V i i \ n II in in in 1 -i r " * _ . _ • -0 100 200 300 400 500 600 700 800 900 1000 Number of samples Figure 4.7: On-line estimation of Q and r (solid) and Isaksson's method estimates (dashed line) for example 2. A New Approach to Adaptive Kalman Filtering 116 10 8 6 CD 4 cr 2 0 -2 -4 I- I ,' ; I • I : •I ; \\ 1; •• x ; - . . . I'-V I I J V i i 0 100 200 300 400 500 600 700 800 900 1000 Number of samples I o_i i i i i i i i i I 0 100 200 300 400 500 600 700 800 900 1000 Number of samples Figure 4.8: On-line estimation of Q and r (solid) and Isaksson's method estimates (dashed line) for example 2 (continued). A New Approach to Adaptive Kalman Filtering 0.35 0.3 h 0.251-0 100 200 300 400 500 600 700 800 900 1000 N u m b e r of samples Figure 4 . 9 : On-line estimation of Q and r for a non-diagonal Q. A New Approach to Adaptive Kalman Filtering 0.01 0.009 r-0.008 h 0.007 h 0.006 -B 15 E w 0.005 -CD CM 0.004 -0.003 h 0 100 200 300 400 500 600 700 800 900 1000 Number of samples ,l i i i I i i i 1 1 1 0 100 200 300 400 500 600 700 800 900 1000 Number of samples Figure 4.10: On-line estimation of Q and r for a non-diagonal Q (continued). A New Approach to Adaptive Kalman Filtering 119 4.7 Conclusions In this chapter, we presented an E M approach for adaptive K a l m a n filtering of a time-varying system. A recursive scheme based on the off-line method was derived. A number of fixed-point smoothers running in parallel is utilized to estimate on-line the noise variance matrices. The developed scheme implies that the estimates are used in the K a l m a n filter. Simulations show that the recursive E M update converges quickly. The advantage of this scheme is that one does not have to wait for some time before starting to use the noise variance estimates in the Ka lman filter. Chapter 5 Summary and Future Work A common approach in dual control is to describe the time-varying parameters by a random walk. Further, the parameters are treated as state variables and therefore can change as rapidly as other system states. This thesis started with the idea that an adaptive K a l m a n filter would then naturally fit into the dual control framework. Two major applications were analyzed: chip refining and paper coating. O n a chip refiner, the incremental gain between motor load and the plate gap is subject to a slow drift due to plate wear and to sudden changes in sign due to pad collapse. A pad collapse can be caused by a change in operating point, or may occur suddenly due to a feed rate or consistency disturbance. The control objective is to regulate the motor load while avoiding pad collapse. The proposed dual controller minimizes a nonlinear loss function, designed especially to deal with the process nonlinearity. Act ive learning is achieved - i n a suboptimal fashion- by including the future estimation variance in the loss function. The main difficulty in this part was to extend the A S O D approach to handle a process with time delay and correlated disturbances. The controller is coupled with an adaptive K a l m a n filter for gain estimation. The controller, tested by computer simulations, proved successful at regulating the motor load and dealing wi th pad collapse. The advantages of the proposed strategy is that it does not require additional heuristic logic or expensive dedicated sensors and that, in the case of pad collapse, the controller rapidly recovers on its own. The disadvantage, however, is that the controller is obtained through numerical optimization. Future work should focus on implementing and testing this strategy on an industrial refiner. The dual control strategy presented for the chip refiner problem was extended to cope with a coupled multivariable process such as paper coating. The challenge presented by C D coat weight control on bent blade coaters is that process gains vary due to blade wear 120 Summary and Future Work 121 and operating conditions and may differ in sign from one CD position to the next and over the life of the blade. Analytical control law, actuator constraints, blade constraints and implementation issues were all discussed. Preliminary results of applying the proposed strategy to an industrial coater were generally successful. The results show the advantages of including probing in the control signal, and that the adaptive Kalman filter was capable of tracking gain variations. The dual controller yielded substantial quality improvement and was able to control the process throughout the entire blade life. Future work will concentrate on moving this controller into a commercial product. A lot of work has still to be done in order to reach that stage. Reducing the number of tuning parameters, generating tuning guidelines, and the development of a graphical user interface easily accessible to non-experts are all required before gaining mill acceptance. Constraint handling may also be the focus of further development aimed at directly optimizing the original loss function subject to different constraints. This work constitutes the first known application of dual control to a multivariable industrial process. The dual control strategy developed here may have applications other than chip refining and paper coating. For example, it may be applied to bioreactors which exhibit similar nonlinear characteristics. The practical theme of this thesis was balanced by the development of a new adaptive Kalman filter scheme. The motivation behind this work was to eliminate the wait period at start-up where the filter is stuck with the fixed initial estimates of the noise variances. A batch E M scheme was first derived, and a novel way to recursify it was developed. The recursive scheme uses a number of fixed-point smoothers running in parallel. It implies that the noise variance estimates are used in the Kalman filter. Simulations show that the recursive E M update converges quickly. Future work should include incorporating the new filter into the dual control strategy, proving the convergence of the recursive scheme and developing a method to force the estimates to be positive semidefinite. Some theoretical issues in dual control are yet to be addressed. The ASOD controller contains a design parameter A that must be carefully chosen to obtain good performance. Although the ASOD loss function is believed to be convex with respect to A, no work has been done to verify this concept or to determine an optimal value for A. On the other Summary and Future Work 122 hand, the development of new suboptimal dual controllers that are easy to implement, can accommodate process constraints, have analytical expressions and are applicable to higher order as well as multivariable processes is anxiously anticipated. Bibliography [1] Edi tor ia l staff report, "Off-machine coater control system reduces coat weight varia-tion," Pu lp & Paper, pp. 150-151, May 1982. [2] B . Al l i son , Dua l adaptive control of chip refiner motor load, P h . D . Thesis, Depart-ment of Chemical Engineering, University of Br i t i sh Columbia, 1994. [3] B . Anderson and J . Moore, Opt imal filtering, Prentice-Hall, 1979. [4] K . Astrom and A . Helmersson, "Dual control of a low order system," Lund Institute of Technology, Department of Automatic Control , Report TFRT-7244 , 1982. [5] K . Astrom and B . Wittenmark, "Problems of identification and control," J . M a t h . A n a l . A p p l . , vol. 34, pp. 90-113, 1971. [6] K . Astrom and B . Wittenmark, Computer controlled systems - theory and design, Prentice-Hall , New Jersey, 1990. [7] K . Astrom and B . Wittenmark, Adaptive control, Addison-Wesley, Mass., 1995. [8] S. Bale, C . Fu , S. Nuyan, J . Mkinen and D . Johnson, "Towards a better coat weight profile control," X I V I M E K O Wor ld Congress, Tampere, F in land, 1997. [9] R. Bel lman, Introduction to matrix analysis, M c G r a w - H i l l , 1960. [10] W . Bliesner, "Basic mechanisms in blade coating," Tappi Journal, vol. 54, pp. 1673-1679, October 1971. [11] R. Boyles, "On the convergence of the E M algorithm," J . Roy. Stat. Soc. B , vol. 45, pp. 47-50, 1983. [12] R . Braatz, M . Tyler and M . Morar i , "Identification and cross-directional control of coating processes," A I C h E Journal, vol. 38, pp. 1329-1339, 1992. 123 BIBLIOGRAPHY 124 [13] Y . Chen, Adaptive paper coating weight control, M.Sc. Thesis, Department of Elec-trical Engineering, University of British Columbia, 1993. [14] K . Christianson and K . Depoy, "Cross direction control: a case study," Tappi Coating Conference, pp. 175-179, 1994. [15] T. Coleman, M . Branch and A. Grace, M A T L A B optimization toolbox user's guide, The MathWorks Inc., 1999. [16] A . Dempster, N . Laird and D. Rubin, "Maximum likelihood from incomplete data via the E M algorithm," J. Roy. Stat. Soc. B, vol. 39, pp. 1.-38, 1977. [17] D. Doerschuk, "An advanced actuator for improved control of cross direction coating weight variations," Tappi Journal, vol. 76, pp. 83-88, October 1993. [18] G. Dumont, "Self-tuning of a chip refiner motor load," Automatica, vol. 18, pp. 307-314, 1982. [19] G. Dumont, "Experience from adaptive control applications," Symp. Contr. Syst. in Pulp and Paper industry, Stockholm, 1986. [20] G. Dumont and K . Astrom, "Wood chip refiner control," IEEE Control Systems Magazine, vol. 8, pp. 38-43, 1988. [21] D. Eklund, "Influence of blade geometry and blade pressure on the appearance of a coated surface," Tappi Journal, pp. 66-70, May 1984. [22] D. Eklund, "Review of surface application," in Fundamentals of Papermaking, vol. 2, Trans, of the 9th Fundamental Research Symposium, pp. 833-870, 1989. [23] M . Ellila, "Automatic coatweight profiling," Paper and Timber, vol. 76, pp. 238-242, 1994. [24] M . Ellila, "Practical studies of coat weight profile control," Tappi Coating conference, pp. 165-173, 1994. [25] A . Enomoto, "Computer control system for coater," Japan Pulp and Paper, vol. 18, pp. 51-60, May 1980. BIBLIOGRAPHY 125 [26] J. Faure, P. Pluvinage and J. Pouyet, "Optimization of mill coater for a low-angle bent blade coating process," Tappi Coating conference, Dallas, pp. 19-26, 1995. [27] M. Feder, Statistical signal processing using a class of iterative estimation algorithms, Ph.D. Thesis, MIT, 1987. [28] A. Feldbaum, "Dual control theory I-IV," Automation and Remote Control, vol. 21, pp. 874-880, vol. 21, pp. 1033-1039, vol. 22, pp. 1-12, vol. 22, pp. 109-121, 1960-61. [29] N . Filatov and H. Unbehauen, "Survey of adaptive dual control methods," IEE Proc. Control Theory Appl., vol. 147, pp. 118-128, 2000. [30] T. Fortescue, L. Kershenbaum and B. Ydstie, "Implementation of self-tuning regu-lators with variable forgetting factors," Automatica, vol. 17, pp. 831-835, 1981. [31] P. Gartaganis, A. Cleland and T. Wairegi, "Blade mechanics of extended blade coaters - theory and practice," Tappi Journal, vol. 61, pp. 77-81, April 1978. [32] J. Ghofraniha, M. Heaven, L. Lee and G. Sidhu, "Optimizing quality and produc-tivity through coordinated CD/MD coater control," XIV IMEKO World Congress, Tampere, Finland, 1997. [33] T. Greiner, "Emerging trends in coating," Paper Trade Journal, pp. 44-46, May 1985. [34] R. Heeg, "On line measurement and CD control of coat weight," Tappsa Journal, pp. 25-30, July 2000. [35] U. Hoeke and H. Sommer, "Setting and controlling CD profiles on blade coaters - a research report," Tappi Coating conference, pp. 353-359, 1990. [36] M. Horner and S. Korhonen, "Computer control of thermomechanical pulping," CPPA 66th Annual Meeting, Montreal, 1980. [37] M. Hurst, "CD coat weight control," Tappi Journal, vol. 78, pp. 248-251, May 1995. BIBLIOGRAPHY 126 [38] A . Isaksson, On system identification in one and two dimensions with signal pro-cessing applications, Ph.D. Thesis, Department of Electrical Engineering, Linkoping University, 1988. [39] A . Jazwinski, "Adaptive filtering," Automatica, vol. 5, pp. 475-485, 1969. [40] J. Kuzmak, "Bevelled-blade coating," Tappi Journal, pp. 72-75, February 1986. [41] D. Lang, M . Elli la, L . McNelles and R. MacHattie, "Requirements for cross direction coat weight measurement and control," Tappi Coating conference, pp. 521-527, 1991. [42] B. Lindoff and J. Hoist, "Suboptimal dual control of stochastic systems with time-varying parameters," Lund Institute of Technology, Department of Mathematical Statistics, Report TFMS-3152, 1997. [43] B . Lindoff, J . Hoist and B. Wittenmark, "Analysis of approximations of dual control," Int. J . Adapt. Control Signal Process., vol. 13, pp. 593-620, 1999. [44] G. Loehrl, "Flooded nip coater," Tappi Blade Coating Technology, Chapter 3, pp. 7-9, 1978. [45] N . Loomba and E. Turban, Applied programming for management, Holt, Reinhart & Winston, New York, 1974. [46] A . Maitelli and T. Yoneyama, "Suboptimal dual adaptive control for blood pressure management," IEEE Transactions on Biomedical Engineering, vol. 44, pp. 486-492, 1997. [47] R. Mehra, "Approaches to adaptive filtering," IEEE Transactions on Automatic Control, vol. AC-17, pp. 693-698, 1972. [48] A . Moghaddamjoo, "Approaches to adaptive Kalman filtering," Control Theory Adv. Technol., vol. 5, pp. 1-18, 1989. [49], J. Owen, A . Roche and K . Miles, "A practical approach to operator acceptance of advanced control with dual functionality," Proceedings of the 1998 Conference on Control Systems, pp. 184-194, 1998. BIBLIOGRAPHY 127 [50] G. Rogers, Matrix derivatives, Lecture Notes in Statistics, vol. 2, Marcel Dekker, 1980. [51] J. Rogers et al., "Automatic control of chip refining," Pulp and Paper Canada, vol. 81, pp. 89-96, 1980. [52] H. Ruckert, "Automatic coat weight control at the blade," PPI, pp. 54-57, Apri l 1980. [53] S. Shapiro, "Computer control of blade coaters," Tappi Blade Coating Seminar Notes, pp. 151-155, 1988. [54] R. Snyder, D. Streifthau and J. Moore, "Computer process control of coating weight," Tappi Journal, vol. 63, pp. 79-81, December 1980. [55] H. Sollinger, "New developments in coating CD profile control," Tappi Coating con-ference, pp. 347-351, 1990. [56] J. Stebel and D. Aeby, "Optimization mode control for double-disc refiner systems: design and experience," Int. Symp. Fundamental Concepts of Refining, Appleton, 1980. [57] M . Tanaka and T. Katayama, "Robust identification and smoothing for linear sys-tem with outliers and missing data," Proceedings of the 11th IFAC Triennial World Congress, Tallinn, USSR, 1991. [58] B. Taylor and D. Doerschuk, "Coat-weight profiling system improves uniformity, increases coater uptime," Pulp and Paper, pp. 45-50, May 1993. [59] B. Taylor and D. Doerschuk, "Mills move to improve coating control," PPI, pp. 50-51, May 1993. [60] D. Titterington, A . Smith and U . Makov, Statistical analysis of finite mixture distri-butions, John Wiley and sons, 1985. [61] H. Tomimasu, K . Suzuki, T. Ogura and P. Luner, "The effect of basestock structure on coating weight distribution," Tappi Journal, pp. 179-187, May 1990. BIBLIOGRAPHY 128 [62] L. Turai, "Analysis of the blade coating process," Tappi Journal, vol. 54, pp. 1315-1318, August 1971. [63] T. Vanya, "Blade coaters," from The Coating Processes, Chapter 3, Tappi Press, pp. 71-81, 1993. [64] R. Vyse and D. Carrington, "Moving towards closed loop control of coat weight," Paper Technology, pp. 24-27, August 1992. [65] R. Vyse, J. King, M . Heaven and H. Mononen, "Advances in cross direction coat weight control," Pulp and Paper Canada, vol. 97, pp. 44-48, 1996. [66] H. Wang, A . Wang and S. Duncan, Advanced process control in paper and board making, Pira International, U K , 1997. [67] W. Windle and K . Beazley, "The mechanics of blade coating," Tappi Journal, vol. 50, pp. 1-7, January 1967. [68] B . Wittenmark, "An active suboptimal dual controller for systems with stochastic parameters," Auto. Contr. Theory Appl ic , vol. 3, pp. 13-19, 1975. [69] B. Wittenmark, "Adaptive dual control methods: an overview," 5th IFAC Symp. Adaptive Systems in Control and Signal Processing, pp. 67-73, 1995. [70] B. Wittenmark and C. Elevitch, "An adaptive control algorithm with dual features," 7th IFAC Symp. Identif. and Syst. Param. Estim., pp. 587-592, 1985. [71] C. Wu, "On the convergence properties of the E M algorithm," The Annals of Statis-tics, vol. 11, pp. 95-103, 1983. Appendix A Approximation of the ASOD Controller h(u) Figure A . l : The loss function J together with its quadratic and nonlinear parts. As stated in [70], the ASOD loss function J (1.27) is nonlinear and therefore can not be minimized analytically. However, if J is expanded around a nominal control w n 0 m up to the quadratic term, an approximate analytical solution can then be found. The question is how to determine unom such that it is as close as possible to the global minimum of J. As illustrated in figure A . l , J has a quadratic part and a nonlinear one. These parts will be called g and h, respectively. J typically has two local minimas and one local maximum. Two important attributes of J are known. First, the global minimum -umin lies on the same side of u0 (the input maximizing h) as the minimum of g. This latter minimum is given by the cautious input uCautious- The point u c a u t ; o u s , however, does not give a satisfactory approximation, since it can fall at a point where J is concave. We may then utilize a second important characteristic of J , that J is convex for all |it| > \UA\ 129 Uo UA U U . ~u "B cautious m i n nom APPENDIX A. APPROXIMATION OF THE ASOD CONTROLLER 130 where uA is the infliction point of h, i.e., it satisfies h"(uA) = 0. The point uA is easily determined analytically. The point UB can then be obtained. It is reasonable to assume that g(uB) = J{UA) since at UB the function h can be neglected. In typical cases umm lies somewhere inbetween uA and uB- Since this section of the curve approximates a quadratic, the nominal point is chosen as •nom 2 (A.l) It has been demonstrated through simulations that the approximate analytical ASOD control law using (A.l) gives results close to those of numerical minimization [70]. Appendix B Proofs for Chapter 3 L(t) was denned in (3.36) as Using (3.4), we get L(t) = E{Gr(t + l)G(t + l)\Yt} L(t) = E {B(t + l)CTCB{t + 1) \Yt} (B.l) (B.2) One useful application of the kronecker product is linear matrix equations. We can write the matrix equation V = FXG (B.3) in the form Then, col (V) = (GT eg) F) col (X) col (L(t)) = E {B(t + 1) <8> B(t + 1) \Yt} col (C T C) which is the expression shown in (3.38), where (B-4) (B.5) E {B(t + 1) eg) B(t + 1) \Yt} bx(t + l)B{t + l) 0 0 b2(t + l)B(t + l) E < 0 0 bn(t + l)B(t + l) Yt (B.6) 131 APPENDIX B. PROOFS FOR CHAPTER 3 132 with E{bi(t+l)B(t + l)\Yt} = ' ( k(t + l\t)b\(t + l\t) \ q _ Q ^ +pblbi(t + l\t) ) 0 ••• I ( h ^ t + l \ t ) \ o { +pbi(t + l\t)) I b^t + i\t)hn(t + i\t) \ { +pbibn(t + i\t) ) (B.7) Proof of (3.39) is obtained in the same manner. Appendix C Estimation of the Noise Filter Zero Recall from Section 3.3 that we have the following linearized model (1 + aq-l)Ay(t) = G{t)Au(t - 1) + (1 + cq-l)Ae(t) (C.l) When the coater is running under manual control, actuator setpoints are fixed (for most of the time). For such periods, one can write {l + aq-1)Ay(t) = (l + cq-1)Ae(t) (C.2) Define Then, and z(t) = (l + aq-1)Ay(t) (C.3) A , = Var Wt)) = ^ 2 ( A A T ) ( l + c2) ' (C4) A1 = Cov(z(t),z(t-l)) = v2 (AAT)c (C.5) where Var denotes the variance and Cov the covariance. This leads to col (A0) = col (Ai) (C6) A simple batch LS procedure can then be applied to estimate c based on (C.6). 133 Appendix D Proofs for Chapter 4 Consider the augmented state vector for t > j x(t + l) x W ( H l ) (D- l ) satisfying x^(t + 1) = x^{t) = x(j) and x^{t + 1) = x^{t) = x(j - 1). Then, w(t) (D.2) x{t + l) F(t) 0 0 x(t) I xW(t + l) = 0 I 0 xW(t) + 0 xW(t+l) 0 0 I 0 and x(t + 1) y(t)= H(t) 0 0 xM(t+l) +v(t) (D.3) xW{t + l) In a similar manner to the derivation of the fixed-interval smoother in [3], the error covariance matrix of the augmented state vector is given by 134 APPENDIX D. PROOFS FOR CHAPTER 4 135 P(t + l\t) p W ( t + P<2>(t + l\t) [ p ( D ( i + i | f ) ] T pWW(t + P ( 1 ) ( 2 ) ( * + ll*) [ p ( 2 ) ( i + 1 | t ) ] T [p^)(t + l\t)]T PW(t + l\t) P(t\t-l) P^(t\t-l) P^(t\t-l) p ( D ( t | t _ i ) ] T PMi)(t\i-i) pw(t\t-i) [pW(t\t-i)]T {pW(t\t-i)]T pw>{t\t-i) ' F{t) 0 0 0 I 0 0 0 I / 0 0 " HT(t)" X 0 I 0 — 0 V 0 0 I 0 0 0 \ 0 I 0 0 0 I / I 0 0 R I 0 0 (D.4) After simplification, pW{t + l\t) = pW(t\t - 1) - [pW{t\t - l)]THT(t) [K^(t) But [pW(t\t - 1)] T - Pa{j,t) and K<$(t) = Ka(j - l,t). This yields P(J,J - l|t) = P(j,j - l\t - 1) - Pa(j,t)HT(t)Kj(j - l,t) (D.5) (D.6) which is (4.69). 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0065524/manifest

Comment

Related Items