AUTOMATED TUNING OF PAPER MACHINE CROSS-DIRECTIONAL FEEDBACK CONTROLLERS by MING ZHANG B.Eng. Xian University of Technology, 1982 M.Eng. Xian University of Technology, 1987 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December, 1997 © Ming Zhang, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of £Uctf}caJ'.CLwL CoAn^^T ^ ^ c ^ T ^ J The University of British Columbia Vancouver, Canada Date DE-6 (2/88) Abstract Automatic tuning of the feedback controllers in a commercial paper machine CD control system is considered in this thesis. A feedback controller tuning tool has been designed and developed for tuning the Dahlin controller and control filter in the Honeywell-Measurex CD control system. This tool includes the process identification module, the disturbance identification module and the controller and filter tuning module. The tool is coded in MATLAB and some algorithms have been implemented as "C" modules embedded into a product currently on the market. The tuning strategy used in this tool is based on modern system identification and controller parameter optimization techniques. A global search method combined with the Levenberg-Marquardt refinement is used for identification of the paper machine time-response model (first order plus delay). To obtain a longer disturbance realization sequence and increase data available for identification of the process disturbance, CD residuals, extracted from bump test data, are used for identifying an Integrated Moving Average disturbance model. The identification method is based on the Recursive Extended Least Squares. Based on the identified process and disturbance models, the Dahlin controller and control filter are tuned to minimize a LQG based performance index. Since disturbance dynamics for a paper machine vary with time and with the operating conditions, the possible changes of the disturbance characteristics should be taken into account in the developed tuning algorithms. In the modified tuning scheme, the controller and filter are tuned in such a way that a performance lower bound is guaranteed for all disturbances in a certain set. The feedback controller tuning algorithms were tested extensively using simulation models, Honeywell-Measurex Devron hardware-in-the-loop paper machine simulator and many sets of mill data. A mill trial was carried out for validation of the disturbance identification and control loop ii variation prediction. Simulation and mill data tests showed that the developed identification algorithms were applicable to real-life mill data and the process variations were predicted with satisfactory accuracy. For better tractability of problem, in the developed algorithms it was assumed that the disturbances in adjacent zones were independent and the actuator response was narrow (about 1 or 2 zones wide). In most cases, the developed algorithms gave reasonable tuning parameters as long as these assumptions held. Validation using mill data also showed that the modified tuning algorithm sometimes overestimated the process and actuator variances and resulted in too conservative a tuning. iii Table of Contents ABSTRACT H LIST O F T A B L E S VI LIST O F FIGURES V H A C K N O W L E D G M E N T S X 1 INTRODUCTION 1 1.1 CONTROL OF PAPER MACHINE CD PROCESS 1 1.2 MOTIVATION 4 1.3 CONTRIBUTION 5 1.4 THESIS OUTLINE 7 2 M O D E L S AND ASSUMPTIONS 8 2.1 PAPER MACHINE CROSS-DIRECTIONAL PROCESS MODEL 8 2.2 FEEDBACKLOOP INMEASUREX CD CONTROL SYSTEM 10 2.3 TUNING PROBLEM 14 3 IDENTIFICATION OF T H E PROCESS TIME-RESPONSE M O D E L 17 3.1 FORMULATION OF THE PROCESS IDENTIFICATION 17 3.2 GLOBAL OPTIMIZATION OF THE LOSS INDEX 20 3.3 LOCAL OPTIMIZATION USING LEVENBERG-MARQARDT METHOD 24 3.3.1 Levenberg-Marquardt Optimization Scheme 24 3.3.2 Computation of the Jacobian 25 3.3.3 Scaling of the Parameters 26 3.3.4 Identification of Integrator Process 29 3.4 IDENTIFICATION FOR DYNAMICALLY FILTERED DATA 31 3.5 COMPARISON OF THE DEVELOPED ALGORITHM AND LINEAR LEAST SQUARES METHOD 34 3.6 CONCLUSION 37 iv 4 IDENTIFICATION OF T H E DISTURBANCE M O D E L 38 4.1 DISTURBANCE MODEL STRUCTURE 38 4.2 DISTURBANCE IDENTIFICATION USING RESIDUAL OF PROCESS RESPONSE IDENTIFICATION 39 4.3 IDENTIFICATION ALGORITHM 40 5 TUNING OF DAHLIN C O N T R O L L E R AND C O N T R O L FILTER. . . 43 5.1 EVALUATION OF THE PERFORMANCE INDEX 43 5.2 OPTIMIZATION OF THE PERFORMANCE INDEX 46 5.3 PERFORMANCE INDEX MODIFICATION FOR ROBUST TUNING 48 5.3.1 Notation and Assumption 48 5.3.2 Definition of a Disturbance Set 49 5.3.3 A Robust Tuning Performance Index 51 5.4 FEEDBACK CONTROLLER TUNING TOOL 52 6 VALIDATION OF T H E TUNING ALGORITHMS , 55 6.1 VALIDATION APPROACH 55 6.2 SIMULATION VALIDATION 58 6.3 VALIDATION USING PAPER MILL D A T A 63 6.3.1 Case 1: Steambox Actuator, Moisture Process 63 6.3.2 Case 2: Slice-Lip Actuator, Weight Process 71 6.3.3 Case 3: Dilution Actuator, Weight Process 78 CONCLUSIONS 86 REFERENCES 89 v List of Tables 3.1 Simulation model parameters and their estimates. The responses are shown in Figure 3.12. The Display filter factor is 0.3 and the noise amplitude is 0.3 33 3.2 Estimates obtained by LLS and by the identification algorithm; Noise to signal ratio is 0.2 36 3.3 Comparison of performance and features between the identification algorithm and LLS 37 6.1 Comparison of the simulation model and estimated model 59 6.2 Parameter values of each component transfer function used in the example 2 61 6.3 Implemented controller parameters 65 6.4 Comparison between the measured process/actuator variation and the predicted process/actuator variation 66 6.5 Parameter values of each component transfer function used in spectral check 68 6.6 Parameters of two tested feedback controllers 73 6.7 Comparison of the predicted 2-sigma and measured 2-sigma for the process variation 74 6.8 Comparison of the predicted 2-sigma and measured 2-sigma for the incremental actuator moves 74 6.9 Parameter values of each component transfer function used in the spectrum check 76 6.10 Implemented feedback controller parameters 80 6.11 Comparison between the measured process/actuator variation and predicted process/actuator variation for Case 3 81 6.12 Parameter values of each component transfer function used in the spectrum check 84 vi List of Figures 1.1 Overview of a paper machine process 2 2.1 Simplified schematics of paper machine CD control system 10 2.2 Simplified feedback loop of Measurex CD control system 11 3.1 Dataflow of me process time-response identification 17 3.2 Loss surface is plotted versus delay and logarithm of time constant. There are too many search points for the time delay and time constant 22 3.3 For some ideal cases the search points can be reduced without missing the global minimum. This decreases computational load in the global search 23 3.4 In this noisy data example loss surface has many local minima. It is possible to miss the global minimum in global search if the search points are not enough 23 dY / 3.5 Before scaling of time constant , the derivative with respect to time constant has smaller amplitude (< 0.14) than other plots 27 3.6 Before scaling of time constant loss surface converges to the minimum unevenly 28 3Y / 3.7 After scaling of time constant the amplitude of increases to 0.38 compared with the previous value 0.14 shown in Figure 3.5 28 3.8 After scaling of time constant loss surface converges to the minimum more evenly than in the case shown in Figure 3.6 29 3.9 Identification of an integrator process: thin line - Process response; thick line - Identified model response; the data is collected from a caliper process 30 3.10 Loss function of identification of an integrator process 3 0 vii 3.11 Time-response identification with a MD filter 31 3.12 Identification with filtered data. Simulation model and estimated model shown in Table 3.1 34 3.13 Comparison of the developed algorithm and LLS. Thick solid - The developed algorithm model response; Dashed - LLS model response 35 3.14 Comparison of the identification algorithm and LLS using digital simulation data 36 5.1 A search domain with arbitrary grid size for finding optimal tuning parameters a and /? 47 5.2 A block diagram of the feedback loop in the Measurex CD control system 49 5.3 Gamma specifies a disturbance set 50 5.4 Main screen of the feedback controller tuning tool 54 6.1 Identification of disturbance model for a simulated process 59 6.2 The process model identification using simulated model 60 6.3 The disturbance model identification using simulated model 60 6.4 Optimality of the Dahlin controller 62 6.5 Identification of the process time-response model 64 6.6 Identification of the disturbance model 65 6.7 Spectral check for the disturbance model: Open-loop spectrum for the model - solid; spectrum computed directly from the data - dashed 67 6.8 Closed-loop process spectral check: Process spectrum predicted using the closed-loop model - solid; spectrum computed directly from the closed-loop data - dashed 69 6.9 Closed-loop actuator spectral check: Actuator spectrum predicted using the viii closed-loop model - solid; spectrum computed directly from the closed-loop data - dashed 70 6.10 Identification of the process temporal model using the controller tuning algorithm 71 6.11 Identification of the process disturbance model 72 6.12 Comparison between the predicted spectrum and measurement spectrum of the open-loop process variations 75 6.13 Comparison between the predicted process spectrum and measurement spectrum on closed-loop 77 6.14 Comparison between the predicted spectrum and measurement spectrum for the actuator variations 77 6.15 Process time-response model identification 79 6.16 Disturbance model identification 80 6.17 Spectral check for the disturbance model 82 6.18 Closed-loop process spectral check 83 6.19 Closed-loop actuator spectral check 85 ix Acknowledgment I would like to thank my supervisor Dr. Guy A. Dumont for his teaching and guidance during my study at the University of British Columbia. I would also like to thank my thesis co-supervisor Dr. Dimitry Gorinevsky who, together with Dr. Dumont, gave me valuable advice and careful correction to my thesis. I wish to express my thanks to Honeywell-Measurex Devron Inc. for providing supervision, computer facilities, data and assistance in development and testing of the feedback controller tuning tool. Special thanks are given to Peter Tremblay who helped to carry out a mill trial and to Johan Backstrom, Scott Morgan, Michael Ward, Cameron Marr and Stephen Chu who provided mill data for validation of the developed feedback controller tuning algorithms. I am very appreciative of the advice and discussion of Chris Sung, Jahan Ghofraniha and my colleagues in the Pulp and Paper Centre. I would like to thank Dr. Dumont and Networks of Centres of Excellence for providing financial support. I gratefully acknowledge the IRAP Support through Technology Enhancement Grand - 42785W to Honeywell-Measurex Devron. And to my family, Jichun and Sai. Chapter 1: Introduction Chapter 1 Introduction 1.1 Control of Paper Machine Cross-Directional Processes Paper machines produce a continuous web of the paper from pulp stock. As shown in Figure 1.1 the headbox delivers the stock onto a mesh conveyor, where the paper moisture is reduced through gravitation and with the help of vacuum drainage elements underneath the wire, thus forming the sheet. Further dewatering takes place as the paper sheet passes through the press and dryer sections. Finally, paper properties such as basis weight, moisture and caliper are measured with a scanning gauge and the sheet is wound up onto a reel. The paper quality is assessed in two dimensions, the machine direction (MD), i.e. the direction in which the paper sheet moves, and the cross-direction (CD), i.e., across the paper web. The goal of paper machine control systems is to maintain the paper properties on target in both MD and CD. CD profiles of paper sheet properties are controlled by various CD actuators. Each type of CD actuator includes a set of identical actuators usually located at evenly spaced points along the cross-direction. Depending on the application and the actuator type, there can be 20 to 300 individual actuator units in one CD actuator set. An example of important CD actuator is the weight actuator, which adjusts the stock distribution across the machine by changing the opening of different sections of the slice lip in the headbox. Sensor measurements are located at a distance down the machine-direction from the actuation. Due to the high cost of sensors, a limited number of sensors (1-2) measures only a zigzag portion of the paper sheet. From this limited number of measurements, the entire sheet profile must be estimated at each sampling time for feedback 1 Chapter 1: Introduction control. This estimation can be performed in a straightforward manner using Kalman filtering and averaging techniques [24]. The control problem is to calculate the actuator moves based on the estimated profile at each sampling instant. Figure 1.1: Overview of a paper machine process For better control of papermaking processes, various advanced control strategies have been proposed and have achieved a certain degree of success [1][28][30][31]. However, most CD control systems still use some kinds of well established simple controllers, such as PI, PID, Dahlin and so on. Tuning of those controllers plays an important role in reducing impact of the process variability on the product uniformity and in ensuring that the process is operated at the chosen target. Industrial experience has shown that controller auto-tuning is a highly desirable and useful feature. 2 Chapter 1: Introduction Controller auto-tuning has been an important research topic for a long time and various tuning strategies have been proposed in the literature. The Ziegler-Nichols method for tuning PID regulators [2] is a popular and widely accepted scheme. It is based on detection of the critical gain and critical period and a quarter decay criterion for controller parameter design. The Dahlin controller is a well known dead-time compensator and is widely used in industries. Its tuning requires the process model being identified and the closed-loop time constant being chosen. E.B. Dahlin [3] gave some guideline for tuning the closed-loop time constant. G.A. Dumont presented a thorough sensitivity analysis of Dahlin controller when subject to modelling errors [4]. He also proved optimality of Dahlin controller for a first order plus delay process perturbed by a first-order integrated-moving-average disturbance. Cluett and Wang [5] proposed a PID controller tuning method based on a specification of the desired control signal trajectory. In this scheme the users are provided with a tuning parameter that specifies the closed-loop response speed. It does give the users one more degree of freedom to chose the closed-loop performance, but it also increases the tuning burden. Commercial single-loop controllers with an auto-tuning feature have been available since the early 80's. They are used in most of typical industrial processes, including temperature, flow, pressure, level and pH control loops. One example of such auto-tuning controller is SattControl ECA400, announced by SattControl in 1988. It is based on Astrom and Hagglund research work [6]. The auto-tuning of the ECA400 is performed using the relay method in the following way. The process is brought to a desired operation point, either by the operator in manual mode or by a previously tuned controller in automatic mode. When the loop is stationary, the operator presses a tuning button. After a short period, when the noise level is measured automatically, a relay with hysteresis is introduced in the loop, and the PID controller is temporarily disconnected. The hysteresis of the relay is determined automatically from the noise level. During the oscillation, the 3 Chapter 1: Introduction relay amplitude is adjusted so that a desired level of the oscillation amplitude is obtained. When an oscillation with constant amplitude and period is obtained, the relay experiment is interrupted and Gpficoo), i.e., the value of the transfer function at oscillation frequency a>o, is calculated using the describing function analysis. The ECA400 controller also has other adaptation features, such as gain scheduling, adaptive feedback and adaptive feedforward. It has been successfully applied to a wide range of industrial processes. But this tuning approach might not be a good solution for paper machine CD control, where high quality control and low variation requirements must be satisfied. Since this tuning method requires to generate an oscillation in the control loop, such oscillation propagates on 2 dimensional paper sheet and its amplitude is hard to confine within certain limits, as even a small amplitude oscillation could cause a serious degradation in the paper product quality or a sheet break. 1.2 Motivation This work was partially supported by Honeywell-Measurex Devron Inc. and was focuses on real-life issues of tuning the feedback controller in Honeywell-Measurex CD control systems. Parameter settings of Honeywell-Measurex CD control systems may crucially influence the paper product quality. In particular, tuning of the Dahlin feedback controller and the control filter in Honeywell-Measurex CD control systems plays an important role in compensating the process variations and mamtaining high quality product. The current tuning method is based on manual adjustment of the feedback loop parameters by trial and error. This is time-consuming and can result in sub-optimal operation, therefore it is very important to equip Honeywell-Measurex CD control systems with a reliable automated controller tuning tool. This would reduce dependence of the field personnel on specialized process control expertise and improve product quality. 4 Chapter 1: Introduction The objective of this thesis was thus to develop a feedback controller tuning tool for the Honeywell-Measurex CD control systems. This tool targets papermaking processes and ensures the continuous production process without being interrupted by a special identification test. In particular, the following design specifications were accepted. (1) The tool should provide robust and automated joint tuning for the Dahlin controller time constant a and the control filter factor /3for the Honeywell-Measurex CD control systems. (2) A reliable process time-response identification algorithm should be developed and included in this tool. It should be robust to the process variation, which might differ significantly from white noise. (3) The tool should be able to identify a dynamic model of the process disturbance in order to tune the controller optimally. (4) The tuning computed by the tool should provide a satisfactory balance between minimizing process variations and avoiding excessive actuator moves. 1.3 Contribution In this thesis, A feedback controller tuning tool was designed for tuning the Dahlin controller and control filter in the Honeywell-Measurex (HMX) CD control system. This tool includes a process identification module, a disturbance identification module and a controller and filter tuning module. The tool is coded in MATLAB. The process identification algorithm has been implemented as 'C modules embedded into the company's CD alignment and tuning tool (a product currently on the market). 5 Chapter 1: Introduction The developed process identification algorithm has enhanced functionality, robustness and performance as compared with other algorithms used in the industry. It allows reliable and fast identification of first order as well as integrator processes in the presence of strong noise. This algorithm works well with nonuniform sampling time and is able to identify the process dynamics by using either the MD-filtered scanner data or the raw data. In the course of this thesis work, a Recursive Extended Least Squares algorithm for identifying the disturbance model was studied. An entire 2D array bump test data are used to identify an integrated moving average model of the disturbance. This gives more consistent identification results than using zone disturbance sequences. The algorithm validation results show that the disturbance identification algorithm is capable of identifying disturbance dynamics model with satisfactory accuracy. In the developed algorithms, the Dahlin controller and the control filter are tuned to minimize a LQG performance index. The engineering tuning "knob" parameter p enables the users to balance different control objectives. By using a modified tuning scheme, the controller and filter can be tuned in such a way that a performance lower bound is guaranteed for a wide class of disturbance characteristics. The feedback controller tuning algorithms have been tested extensively using simulation models, Honeywell-Measurex Devron hardware-in-the-loop paper machine simulator and a few sets of mill data. A mill trial was carried out for validation of the disturbance identification and control loop variation prediction. The validation results show that the identified disturbance model captures basic characteristics of the process variations. Overall the process variations were predicted with a satisfactory degree of accuracy. 6 Chapter 1: Introduction 1.4 Thesis Outline The thesis is organized as follows. Chapter 2 introduces models and assumptions used in design of the feedback controller tuning algorithms. Chapter 3 describes the method and algorithm used for identification of the process temporal model. Chapter 4 discusses the identification of the disturbance model and presents a Recursive Extended Least Squares based approach. Chapter 5 gives the tuning strategy and formulae for the performance index evaluation. In order to cope with the disturbance model uncertainty, a robust tuning approach is introduced. The procedures and results for validation of the developed algorithms are presented in Chapter 6. Chapter 7 summarizes this research results and proposes future studies in this area. 7 Chapter 2: Models and Assumptions Chapter 2 Models and Assumptions The foundation of this research is presented in this chapter. First, the paper machine CD process model is introduced and the assumptions used in design of conventional CD control systems and the developed tuning algorithms are reviewed. Next, the feedback loop of the Honeywell-Measurex CD control system is described and the model for each of components in the feedback loop is presented. Finally, the controller tuning problems to be solved are formulated, thus providing a starting point for the remaining chapters. 2.1 Paper Machine CD Process Model Paper machine CD control is a complicated and challenging task. It involves measuring a few hundreds of data boxes across the sheet and adjusting dozens of control elements (actuators). Adjustment of each actuator affects not only its corresponding measurement zone on the paper sheet but also that of its neighboring actuators. Therefore a paper machine is inherently a multi-input-multi-output (MIMO) system with different numbers of inputs and outputs. For most of paper machine CD control systems currently used in paper mills, including those of Honeywell-Measurex Inc., a relatively simple control scheme based on profile transformation and assumption of actuator independency is adopted. Figure 2.1 shows the block diagram of a traditional CD profile control system. In Figure 2.1 a scanner system measures a particular paper property across the sheet and produces a high resolution measurement profile. This high resolution profile is then pre-processed and filtered to remove invalid measurements and high frequency noise. The filtered profile is compared with the target and the error profile is determined. To make the 8 Chapter 2: Models and Assumptions number of process outputs equal to the number of inputs, the high resolution error profile is ( transformed into a low resolution profile through a mapping algorithm, producing a profile with the same number of error elements in the error profile as the number of actuators being controlled. Based on the low resolution error profile the controller determines control moves that correct the profile deviation from the target. Such a mapped CD process can be described as a square MIMO system, where the numbers of inputs and outputs are the same: y(t) = g(z-l)Gu(t) + i(t), (2.1) where z_1 is a backward shift operator; t is discrete time; the output vector y(t) has elements (jy.(r), i-l,...,n} such that each element is an averaged paper property measurement associated with the corresponding actuator zone; u(t) is n by 1 actuator input vector; g(t) is n by 1 process disturbance vector, each of its elements representing the lumped amount of process disturbance at a corresponding actuator zone; G, known as the interaction matrix, is a n by n square matrix and each column of G can be formed with the CD response shape of its corresponding actuator; g(zl) is a scalar-valued process time-response model and can usually be described as a first order, delayed process transfer function. In most industrial CD control systems, control moves are determined by treating each actuator as an independent element, and then the subsequent control moves are adjusted taking into account spatial interactions using some decoupling techniques. For tractability of the tuning algorithm this thesis assumes that the actuator responses do not overlap (in (2.1) G is an identity matrix) and that each of actuators has the same dynamic response which can be described with a first order plus delay model. It is further assumed that the disturbance dynamics for each of actuator control zones are identical and there is no disturbance correlation between adjacent control zones. 9 Chapter 2: Models and Assumptions Targets Actuator resolution profile Profile processing JL High resolution profile Scanner system Profile transformation Paper machine Mapped controller Actuator interface Figure 2.1: A simplified paper machine CD control system 2.2 Feedback Loop in Measurex CD Control System Under the above CD control assumption, the system is simplified to n parallel independent feedback loops, which are identical. One of such loops is shown in Figure 2.2. It consists of a first order plus delay process model, a Dahlin controller, a first order control filter, a display filter (minislice filter) and a disturbance model. In this loop, the tuning of the Dahlin controller and the control filter needs to be determined. This tuning depends on the dynamics of both the process and the disturbances. The following gives a brief description for each component in the feedback control loop. 10 Chapter 2: Models and Assumptions Noise Shaping Filter sp Dahlin Controller u Process + y—• r V Control Filter Display Filter Figure 2.2: One of simplified closed-loops of Measurex CD control system Process time-response model The process model relating actuator setpoint to the scanned measurement is described by a first order plus time delay transfer function of the form kp(l-a)z 1-oz -1 -d-\ (2.2) where z is the backward shift operator; kp is the process gain, d is the discrete process time delay and a defines the process time constant; d and a can be determined through continuous to discrete parameter conversion: d=mt(Td/Ts), -r.ir. a-e (2.2*) where int denotes rounding a real towards to nearest integer, Ts is the scan time, Td is the continuous process time delay, and Tr is the continuous process time constant. 11 Chapter 2: Models and Assumptions Dahlin feedback controller The feedback controller used in the Honeywell-Measurex CD control system is a Dahlin controller. Its tuning requires a desired closed-loop time constant a being chosen and a process model being identified. From Figure 2.2, it can be seen that the 'process model' used in design of the Dahlin controller should be a first order model whose dynamics is equivalent to that of the process model plus the first order filter. Thus, the discrete-time transfer function of the Dahlin controller can be written as: kc(\-rz-x) D ( Z ) "(l-z- 1 ) [ l + (l- a)z- 1+...+(l-a)z-' /]' ( 2 3 ) where a defines the desired discrete closed loop time constant; d is estimate of the discrete process time delay; x is the equivalent time constant of the filtered process (after the control filter). The controller gain, kc is defined as: where kp is estimate of the process gain. The discrete-time parameters a and r can be obtained from respective continuous-time parameters -T IT a = e hn°, T = e~T',T', (2.4*) T T=T- ln(l-/?)' 12 Chapter 2: Models and Assumptions where Ta is desired closed-loop time constant (in seconds), Tz is the filtered process time constant (in seconds), also called the control time constant, T, is the scan time and ft is the control filter factor. The Dahlin controller is a widely used process delay compensator. It has a simple structure and includes a single tuning parameter. This tuning parameter a specifies the response speed of the closed-loop system to a setpoint change. The smaller it is, the faster the process response is. However overshoot and oscillation may result due to a too small a. The process disturbance characteristics also affect the optimal tuning of a. As a rule of thumb, a is detuned if the controlled process is influenced by large high frequency disturbances. For a first order process perturbed by a first-order integrated-moving-average (IMA) noise, [4] proves that if the tuning parameter a equals to Ci (the only parameter in the noise model) the Dahlin controller is equivalent to a minimum variance controller. Control filter For HMX CD control system, the scanned measurement profile is filtered by a control filter to remove high frequency noise from each of the profile databoxes. The control filter is an exponential filter, whose transfer function can be written as follows: where B is the control filter factor and is one of the tuning parameters. 13 Chapter 2: Models and Assumptions Display filter The display filter in Measurex CD control system affects only the displayed and measured process profiles and has no effects on the actuator moves (see Figure 2.2). Unless the display filter is disabled, it would filter the logged high resolution process profiles used by the identification algorithms, so its effect on the identification of the process and disturbance models should be taken into account. The display filter is also an exponential filter, whose transfer function can be written as follows: where <p is display filter factor. 2.3 Tuning Problem In the feedback control loop of Honeywell-Measurex CD control system in Figure 2.2, the structure of the feedback controller including the Dahlin controller and control filter is given. The main tuning problem considered in this work is to identify the process time-response model, identify the process disturbance model, and determine the two tuning parameters: the desired closed-loop time constant a of the Dahlin controller and the control filter factor B. In general, identification of the process time-response model includes data collection in the identification experiment, model structure selection and use of an identification techniques to determine the model parameters. In practice, identification experiment is performed as an actuator bump test - by applying a step change to actuators and logging the process response data. These data are used for fitting a chosen model, usually a first order with time delay model. This approach is currently standard for the industry. Herein it is assumed that the data collection procedures of 14 Chapter 2: Models and Assumptions the bump test and the process model structure are given and the focus will be on the development of an algorithm for identification of the process time-response model. The detailed problem statement and developed algorithms for identification of the process time-response model are presented in Chapter 3. In order to tune the controller and filter optimally, a disturbance model should be established. It is highly desirable to determine this model from the same set of bump test data as used for the process model identification, without resorting to additional data collection. The problem of identifying the disturbance model from bump test data can be formulated as follows. The data available for identification of the process disturbance model are the residuals of the CD process model identification, obtained by subtracting the response of the process model to the actuator bump from the measured bump response i(t) = Kt)- g(z-l)Gubump(t), (2.7) where the bump test time index t = l,2,...,m (scan number); ubump(t) is the applied actuator setpoint vector at time t; y(j) is the measured bump response at time t. The process time-response model gfz1) can be identified using the method presented in Chapter 3. The spatial response model G can be identified with an available method [7]. Chapter 4 will discuss the disturbance model structure and how to use 2-dimensional disturbance realization £ (t), t = 1,2,..., m, to identify the disturbance model. Based on the feedback loop component models represented by (2.2), (2.3), (2.5), (2.6) and a disturbance model to be developed in Chapter 4, the Dahlin controller and the control filter are tuned to rninimize the following performance index: 15 Chapter 2: Models and Assumptions J(a,B) =E{Y2 +pAu2} -> min, (2.8) where a is the desired closed-loop time constant of the Dahlin controller; B is the control filter factor; E(.) denotes mathematical expectation of a random variable; Efy2) is the process output variance and E(Au2) is the incremental control move variance; p is the control weighting factor that penalizes excessive actuator moves. The tuning objective is to determine optimal a and /? values which rninimize the process variations while keeping the control action within acceptable bounds. Chapter 5 will discuss how to evaluate the performance index (2.8) and obtain the optimal solution. 16 Chapter 3: Identification of the Process Time-response Model Chapter 3 Identification of the Process Time-response Model The Dahlin controller in the feedback loop of the Honeywell-Measurex CD control system (see Figure 2.2) is model based and its tuning requires to set the parameters of the process such as gain, time delay and time constant. In this chapter, an identification method for determination of those process parameters is developed. This method combines a global search algorithm with Levenberg-Marquardt refinement of the parameter estimate to achieve reliable, fast and robust estimation. 3.1 Formulation of the Process Identification Figure 3.1 illustrates the dataflow in the time-response identification from CD bump test data. Time response CD bump test data HMDI *| identification tool Process identification algorithm -> Process gain -> Time delay -> Time constant _^ Estimated MD response Figure 3.1: Dataflow of the process time-response identification In Figure 3.1 a Honeywell-Measurex Devron (HMDI) identification tool is shown to provide an estimate of the process time-response from a set of bump test data. This estimate is based on a CD response model identified by the tool. The developed process identification algorithm uses the time 17 Chapter 3: Identification of the Process Time-response Model response to estimate process gain, delay and time constant needed for control. The developed algorithm can be also used for identification of the MD actuators. In this case, the MD response of the process is directly supplied to the algorithm. There are some difficulties and issues related to the identification of the time-response model for a paper machine process. (1) The bump test data used in the identification are very noisy due to process variation and disturbances. Process variations are typically caused by fluctuations of volume and consistency of stock flow and local variability of the headbox consistency. The developed algorithm should be robust to process variations and noise and work well with various data sets. (2) The scan time may be not uniform due to scanner standardization which occurs every several scans. The developed algorithm should be able to use uneven scan time to perform the identification. (3) The collected bump test data usually are filtered by a low-pass display filter that removes high frequency noise. This should be taken into account during the development of the algorithm. The problem of identifying the time-response of a paper machine can be formulated as follows. The bump test time-response data available to the identification algorithm include two arrays: Y array of the process response samples and T array of the sampling times (not necessarily uniformly spaced). r (3.1) T-\t t t 1 18 Chapter 3: Identification of the Process Time-response Model where y j is the process sample at the scan time t} and N is the number of collected scans in the bump test. The test data also include time tb at which the actuator bump was applied. The amplitude £7 of the bump is further assumed to be unity. This can be always achieved by scaling the data. The collected data are matched against a first order with time delay model. The response Y generated by the model has the form: f b, for t}<.tb+td ^ ^ ^ • • • • M ^ = 1 ^ . ^ ^ ] + ^ for t j > h + t d (3-2) where kp is the process gain, tb is the bump start time, td is the process time delay, tr is the process time constant and b is the response baseline average. The objective is to find the optimal estimate of the parameter vector 9 — [kp td tr b] which minimizes the loss function J = \\Y-Y(d)\\2=\\Z(e)\\\ (3.3) where || • || denotes the Euclidean norm of a vector and Z(0) is a model fit error. Minimization of the loss function (3.3) is a typical nonlinear least squares problem that can be solved by using standard computational methods, e.g., see [11]. In the developed algorithm a global search and Levenberg-Marquardt refinement of the parameter estimate are combined to achieve reliable, fast and robust estimation. The global search algorithm is described in Section 3.2. The Levenberg-Marquardt local optimization is discussed in Section 3.3. Section 3.4 describes a modification of the algorithm needed to take MD pre-filtering of the data into account. Section 3.5 contains comparison of the performance of the developed 19 Chapter 3: Identification of the Process Time-response Model algorithm vs. the Linear Least Squares method. Finally, conclusions for this chapter are drawn in Section 3.6. 3.2 Global Optimization of the Loss Index The basic idea of global optimization is to calculate the loss function (3.3) for a range of different combinations of the parameters kp, td, tr and b, and then find their optimal values by locating the niinimum of the loss function. The purpose of the global optimization is to provide a good initial guess for the iterative Levenberg-Marquardt algorithm and prevent Levenberg-Marquardt iteration from converging to a local minimum, which is possible when a bad initial guess is chosen. Choice of various combinations of the parameters kp, tj, tn and b in the global search affects the accuracy and computational load of the global optimization. The larger the number of combinations of the parameters tried is, the more accurate the estimate is, and the larger the computing effort is. The results of the global optimization will be passed to the Levenberg-Marquardt algorithm as an initial guess and further refined by the Levenberg-Marquardt iterations, therefore at this stage the accuracy of the results is not a major concern and our priority, instead, is to reduce the computational load, while ensuring that a vicinity of the global niinimum is attained. Note that the parameters b and kp appear in a linear way in the modeled response (3.2), thus, the estimates of these parameters can be obtained for each given pair of the parameters td and tr. In particular, in (3.2) b can be calculated by averaging the baseline of the process response. (3.4) 20 Chapter 3: Identification of the Process Time-response Model where yt is the value of the time response at ti, and n is the number of the baseline scans. Since (3.2) is linear with respect to kp, (3.3) can be rewritten as For CD control of paper machine the process time delay is usually between 1 to 10 scans. Therefore, 0 and 10-7^ can be taken as the lower and upper boundaries of the global search domain for the time delay parameter td, where Ts is the scan time. The number of search points taken from this search domain determines the accuracy and the computing effort in searching the optimal time delay. Five to ten search points would usually satisfy the requirement. For the global search of the time constant, the initial value and final value can be chosen as 0.1- Ts and 40- Ts. The number of the search points can be chosen as 6 to 12. In order to customize the algorithm to identification of different processes, the initial value, step length and final value of time delay and time constant in global search can all be adjusted. The above mentioned settings for the global search of time delay and time constant have proved to be robust and applicable to most cases of the studied paper machine bump test data (about 20 sets). In some special cases the settings for the global search can be modified accordingly, for example the number of the search points for the time delay td and time constant tr may need to be increased in the presence of strong process noise. J=\\Y,-kpY2\\\ (3.5) where Yx = Y- band Y2 = A least squares estimate for the gain kp can be calculated from (3.5) as (3.6) 21 Chapter 3: Identification of the Process Time-response Model 10 Figure 3.2: Loss surface is plotted versus delay and logarithm of time constant. There are too many search points for the time delay and time constant. One can get a feel on choosing the number of the search points by examining surface plots of the loss function (3.5) dependence on the time delay fa and time constant tr. In Figure 3.2, the loss function has only one minimum (global minimum). This is an ideal case, where the signal/noise ratio is large. The choice of 20 search points for each of the time delay and time constant seems unnecessary. On contrast, Figure 3.3 shows that although there are only 5 search points for the time delay and 7 for the time constant the global minimum still can be spotted, while the computational load is much less than for the case shown in Figure 3.2. Yet for some bad data sets the number of the search points need to be large enough (Figure 3.4). 22 Chapter 3: Identification of the Process Time-response Model 6^ 0 10 10 Time dela^ (in scans ) -10 -10 lii( Time const ) Figure 3.3: For some ideal cases the search points can be reduced without missing the global minimum. This decreases computational load in the global search. 10 Time delay Ln( time const ) Figure 3.4: In this noisy data example loss surface has many local minima. It is possible to miss the global minimum in global search if the search points are too few. 23 Chapter 3: Identification of the Process Time-response Model 3.3 Local Optimization Using Levenberg-Marquardt Method The global search discussed in Section 3.2 provides only a coarse guess of the optimum location. An accurate computation of the optimal parameters can be done by one of the available standard nonlinear optimization methods. Levenberg-Marquardt method is one of the well acknowledged solutions. Other well known methods, such as Gauss-Newton or Gradient Descent, can be considered as special cases of the Levenberg-Marquardt scheme obtained for special values of the scheme parameters. 3.3.1 Levenberg-Marquardt Optimization Scheme The Levenberg-Marquardt scheme can be derived as follows. Consider a fit error vector Z[G]=Y-Y{d) (3.7) We expand Z [0] in a first order Taylor series about 9^ as Z [ 0 ] * Z w + | k 0 - 0 w ) , 0.8) where # w and Z ( f c ) are the model parameter vector 0 and the model fit error vector Z at the dZ rt, iteration k, respectively. The gradient matrix should be evaluated at 0K . By substituting (3.8) into the problem (3.3), we obtain a typical quadratic optimization problem that is solvable analytically. The Lenvenberg-Marquardt method imposes a limit on the optimization step length of the form || 2 <d2, where d is the maximal allowed step length. By using the Lagrange multiplier method to take into account this step limitation, we obtain the Lenvenberg-Marquardt update of the optimal solution in the form: 24 Chapter 3: Identification of the Process Time-response Model fll + (dZ\ dO \dO) (3.9) where JJ. is a scalar Lagrange multiplier for the step limitation and 7 is a unity matrix. To dZ dY implement the algorithm (3.9) we need to compute the Jacobian matrix ~^Q~=~ ~~^Q~ f° r m e error vector. This is discussed in the next section. 3.3.2 Computation of the Jacobian dZ Given the process response model (3.2), the Jacobian matrix can be calculated as dZ dY d6~~ d9~ dY dY dY dY dkp dtd dtr db (3.10) where dY dk„ if t^t„+td (3.11) dY dtA 0, if t^tb+td k p t - ^ - « - ' ^ , if t]>tb+td (3.12) dY dt 'f tj<tb+td ^ { - K ^ r h - t y ^ - ^ ' , if tj>tb+td (3.13) dY_ db = 1. j=\,2,...,N (3.14) 25 Chapter 3: Identification of the Process Time-response Model Equations (3.11)-(3.14) allow us to compute the gradient of the loss function analytically. The Jacobian matrix (3.10) plays an important role in the Levenberg-Marquardt algorithm (3.9). The introduce a method of improving the condition number of the Hessian matrix by scaling the parameters in the process model. 3.3.3 Scaling of the Parameters Nonlinear scaling of the parameters in the process model can improve convergence and impose constraints on the values of the parameters. In our approach, the process time constant tr is transformed into a new variable trX called scaled time constant: trl = — ln tr. Note that for any real value of r,, the time constant remains positive. In this way the constraint tr > 0 on the time constant is automatically enforced. With the parameter transformation (3.15), the process model becomes: dzT dZ algorithm will not work if the Hessian ee ee —— is close to a singular matrix. Next section will (3.15) (3.16) dY The derivative can be found by using the chain differentiation rule as dY _ dY dtr (3.17) dtrX dtrdtrX 26 Chapter 3: Identification of the Process Time-response Model dY dY dY The derivatives -—— ——and - — can be obtained by substituting (3.15) into (3.11), (3.12) dg dtd db and (3.14). 0.5 SY/dk, -0.2 -0.4 -0.1 3Y/dti -0.2 10 —r— 30 40 50 10 20 30 40 50 10 20 30 40 50 10 20 30 40 Scans 50 60 60 60 60 70 70 70 70 80 80 80 80 Figure 3.5: Before scaling of time constant dY/dtr, the derivative with respect to the time constant, has smaller amplitude (<0.14) than other plots. On the other hand, the above scaling facilitates convergence of the estimates by improving conditioning of the Jacobian matrix • Scaling of time constant improves the conditioning by making the Jacobian matrix column vectors have more uniform norm. Figure 3.5 shows four plots that corresponds to four columns of the Jacobian matrix represented by the equations (3.11)-(3.14). As shown in Figures 3.5 and 3.6, before the scaling, the amplitude of the derivative of Fwith respect to the time constant tr is small (<0.14), so the loss function converges to the minimum slower in the time constant direction than in the delay direction. After the scaling, the 27 Chapter 3: Identification of the Process Time-response Model 0 o Time const, in Scans Figure 3.6: Before scaling of time constant loss surface converges to the rninimum unevenly 1 0.5 SY/aic,, o ( 0 -0.2 3Y/&d -0.4 1 0.4 l # 1 1 1 1 1 1 ) 10 20 30 40 50 60 70 80 V I ) 10 20 30 40 50 60 70 80 I 1 1 1 1 0.2 SY/cXr 10 20 30 40 50 60 70 80 1 ' 1 1 ' 1 1 8Y/db 0 ( i < i i i i i ) 10 20 30 40 Scans 50 60 70 80 Figure 3.7: After scaling of time constant, the amplitude of dYldtx increases to 0.38 compared with the previous value 0.14 shown in Figure 3.5 28 Chapter 3: Identification of the Process Time-response Model amplitude of the derivative of Y with respect to trX is up to 0.34 and the loss function changes more evenly in the two directions. See Figures 3.7 and 3.8. Figure 3.8: After scaling of time constant, the loss surface converges to the minimum more evenly than in the case shown in Figure 3.6 3.3.4 Identification of Integrator Process An integrator process is characterized by continuously increasing output in response to a step change in the input. The control of caliper by an inductive heating actuator can be considered as a case where the process behaves as an integrator. An integrator process presents a difficult problem for the time-response identification. For example, an unstable model may be obtained if a standard Linear Least Squares method is used. For an integrator process, the step response never levels off (no steady state is reached). 2 9 Chapter 3: Identification of the Process Time-response Model 2 i 1 1 1 1 1 1.5 -1 -0.5 -0 --0.5 ' 1 1 1 1 1 0 1 0 20 Scans 30 40 50 Figure 3.9: Identification of an integrator process: thin line - Process response; thick line -Identified model response; the data is collected from a caliper process 1 0 . Delay in Scans 0 -20 -Ln(Time Const.) Figure 3.10: Loss function (model fitting error) of identification of an integrator process Our approach is based on the assumption that the integrator process can be approximated by a first order model with a large time constant. The identification method is still the global search plus 30 Chapter 3: Identification of the Process Time-response Model Levenberg-Marquardt update. Note that this method guarantees that a stable model will be obtained. With the use of time constant scaling in (3.15), Levenberg-Marquardt method works especially well for the estimation of process with a large time constant. Figures 3.9 and 3.10 illustrate that the developed tool works well for the identification of an integrator process. 3.4 Identification for Dynamically Filtered Data In the Honeywell-Measurex control systems, the process response data collected in a bump test may be filtered by the MD display filter. The transfer function of the display filter is given by (2.6). Using the filtered profile to perform the identification will result in a model including u Y J Process Display filter » Model 7 J Display filter Figure 3.11: Time-response identification with a Display filter the dynamics of the process as well as the display filter. Clearly this model is not what we want. One possible solutions is to restore 7 by inversely filtering Yf (use the inverse display filter), but this will amplify the process noise and may cause problems in the parameter estimation. 31 Chapter 3: Identification of the Process Time-response Model Our method is illustrated in Figure 3.11: A composite model (first order model + display filter) is used to simulate the filtered process output Yf. The parameters of the first order model are computed by minimizing the modified loss function of the form j(9)=\\Zf(0)\\2=\\Yf-Yf(9)\\2 ->min, (3.17) where Yf is the filtered time-response of the process Yf=[yf„yf,2,...,yfM]T (3.18) and Yf is the filtered model response (3.2) Yf=[yfA,yf,2,...,yfMf (3.19) Both filters in Figure 3.11 should have the same filter factor q>. The previously developed algorithm can be modified to perform minimization of the loss function (3.17) instead of (3.3) in the following way. Both the global search and local update algorithms of Section 3.2 and 3.3 should use a filtered model response (3.19) instead of 9(9), where the components of the filtered response can be found as yfJ=(\-q>)yf^+9yi (3.20) The filtering is a linear operation which can be represented as Yf = FY, where F is an dlf dl appropriate matrix. Thus, in the Levenberg-Marquardt update we obtain —— = F~— 32 Chapter 3: Identification of the Process Time-response Model k l , d k r . (i-XdYf~ dkf_ ' dY~ j <P) j-+ <p -i dkp (i-d t * _ ' dY j <P) j-+ <p l k l ~dYf~ (i-' df <P) j-+ <p l kJ, db (i- db 'df <P) + <p db j j--i (3.21) (3.22) (3.23) (3.24) dY where the components of the Jacobian matrix —— are given by (3.10). Gain Time delay Time constant Estimated model 0.98 3.95 9.97 Simulation model 1 4 10 Table 3.1: Simulation model parameters and their estimates. The responses are shown in Figure 3.12. The display filter factor is 0.7 and the noise amplitude = 0.3 Simulation shows that the developed algorithm works well in the presence of a display filter. In Figure 3.12 thin solid line represents the filtered measurementsYf that are generated by a simulated process(a first order model plus a exponential filter). The dashed line shows the filtered mode response Yf. The thick solid line stand for the response of the estimated first order model (Y) whose parameters are given in table 3.1. Figure 3.12 shows that the filtered model response closely follows the filtered measurements. Table 3.1 shows that the estimated model gives good 33 Chapter 3: Identification of the Process Time-response Model match with the simulation model. All three estimated parameters are close to their true values. Each is within a bound of 2% of its true value. 1.2 0.81-0.6 0.4 0.2 0 -0.2 Model response before filtering Filtered processi response Filtered model Response 20 40 Scans 60 80 100 Figure 3.12: Identification with filtered data. Simulation model and estimated model shown in Table 3.1.. 3.5 Comparison of the Developed Algorithm and Linear Least Squares This section provides performance comparison between the developed identification algorithm and a standard Linear Least Squares (LLS) algorithm. Linear Least Squares is a widely used method for the parameter estimation. When the controlled process dynamics are linear and process noise is uncorrelated, LLS will give a good estimate. Unfortunately, real papermaking processes do not provide such ideal conditions. Therefore LLS is not a best choice for solving our problem. Figure 13 compares the performance 34 Chapter 3: Identification of the Process Time-response Model between LLS and the developed identification algorithm. The test data were collected from a paper mill (heavy grade paper). 0.8 - / Y S , 0.6 -0.4 -0.2 V 0 •0.2 -0.4 0 10 20 30 40 50 60 70 80 Scans Figure 3.13: Comparison of the developed algorithm and LLS. Thick solid - the developed algorithm model response; dashed - LLS model response. The controlled property is basis weight of paper sheet. LLS give a larger estimated value of the time constant than the developed identification algorithm does. As one can see, the identification algorithm model seems to describe the process response in a more accurate way. The reason is probably that the process variation is highly correlated. Note that unlike LLS the developed identification algorithm uses a nonlinear optimization technique and is able to provide a more realistic estimate. Figure 3.14 shows the performance comparison of the developed identification algorithm and LLS using the data generated by a digital model whose parameters are given in Table 3.2. Although we assume that the user of the LLS algorithm knows the exact Time Delay (1 scan in this 35 Chapter 3: Identification of the Process Time-response Model case), this does not seem to help the LLS in the estimation of the other parameters (the Gain and Time Constant). time response-thin solid; proposed method-thick solid; LLS-dashed tj/ iv/ if il IB if ill ill // WW J I I I I I I " O 10 20 30 40 50 60 70 80 Time in Scans Figure 3.14: Comparison of the developed algorithm and LLS using simulated data Gain Time constant(scans) Time delay(scans) Simulation model 1 6 1 LLS 0.9614 3.8495 1 (selected by user) The developed algorithm 0.9871 5.3464 .1.1082 Table 3.2: Estimates obtained by LLS and by the developed identification algorithm; Noise to signal ratio: 0.2 Table 3.3 summarizes the comparative performance and features of the above described identification algorithms. Clearly the developed identification algorithm is superior to LLS in many aspects. 1 0.8 0.6 0.4 0.2 0 36 Chapter 3: Identification of the Process Time-response Model Linear Least-Squares The developed algorithm First-order process fair good Integrator process identified dynamics can be unstable fair Fractional delay identification no yes Works for uneven scan time no yes MD filtering taken into account no yes Robustness to noise fair good Table 3.3: comparison of performance and features between the developed algorithm and LLS 3.6 Conclusion This chapter has described a process time-response identification algorithm. The identification algorithm is capable of reliable and fast identification of the time-response model in the presence of strong noise. It has enhanced functionality, robustness and performance compared with the linear least squares method. It can be used for first order processes as well as integrator processes (inductive heating control of caliper). It is designed with consideration of scanner standardization, so uneven scan time will not be a problem. It can identify the process dynamics by using either the MD-filtered scanner data or the raw data. The developed algorithm has been tested on many sets of mill bump test data. The results are quite positive and promising. 37 Chapter 4: Identification of the Disturbance Model Chapter 4 Identification of the Disturbance Model A key part of the developed feedback controller tuning tool is the disturbance identification module. The disturbance identification plays an important role in prediction of the process variation and tuning of the Dahlin controller and control filter. In this chapter, the scheme for identification of the paper machine disturbance model is presented. First an Integrated Moving Average (IMA) disturbance model is introduced. This is followed by discussion of how to use the entire 2-D disturbance realization array extracted from CD bump test data to identify the disturbance model. Finally the Recursive Extended Least Squares method used in the disturbance identification module is described. 4.1 Disturbance model structure The disturbances associated with the papermaking process have complicated dynamics and various forms. Besides step-wise load disturbances, such as grade change, stock volume and consistency change, there are substantial periodic disturbances caused by mechanical vibration, hydraulic pulsation and periodic variations in raw material. The periodic disturbances can further be classified as short term (wavelength: 0.1-10 meters), medium term (wavelength: 10 - 2000 m) and long term (wavelength: > 2000 m) components. The control action can eliminate the medium and long term variations and does little to improve the short term variation. Another characteristic of the disturbances is randomness. The actual disturbances can be characterized as a combination of step, periodic, and random components. A suitable disturbance model should reflect generic characteristics of paper machine disturbances and should also be tractable for the automatic controller tuning. For a sheet forming process, such as a paper machine or a plastic film line, modelling the process disturbance in the form of an Integrated Moving 38 Chapter 4: Identification of the Disturbance Model Average (IMA) model [8] proved to be successful. According to our assumption that the disturbance dynamics is the same for all actuator zones and there is no disturbance correlation between adjacent control zones, the IMA model for the disturbance £,(t) in (2.1) is assumed to have the form i(t)=N(z-l)e(t), (4.1) N(z-')=^Q, (4.1a) l — z C(z-l)=l+clz-1+...+c,z-', (4.1b) where e(t) is a vector containing independent white noise elements with zero mean and variance a2 and C(z_1) is a MONIC polynomial of order /. MONIC means that the first coefficient of C(z_1) is 1. In what follows, Cfz'1) will be identified from the measured data. Note that in [4] it is proved analytically that a Dahlin controller is a minimum variance controller for the processes of the form (2.2) if the disturbance structure has the form (4.1a), where / = 1. Furthermore the optimal tuning in this case is a = cx and no filtering (J3 = 1). This allows for a verification of the tuning methods to be developed further. 4.2 Identification using CD residuals The CD identification residuals (2.7) £(t) = [gl(t),...,gn(t)]T will be used for the disturbance model identification. Experience with mill data processing shows that due to limited duration (around 20 to 60 scans) of a common industrial bump test the residual sequences for individual zones £t (t) can be too short to give a consistent estimate of disturbance dynamics. If for all zones, the disturbance sequences can be described by the same dynamics and are uncorrelated, as assumed in (4.1), then each sequence can be 39 Chapter 4: Identification of the Disturbance Model regarded as a different realization of the same random process. In our disturbance identification scheme, the zone disturbance sequences are cascaded to form a single vector where C(z'1) is the monic polynomial (4.1b); e = \e(l), e(mn)]T, is a scalar sequence of white noise with zero mean and unit variance; cr is standard deviation of the white noise. The model parameters to be identified include the coefficients of polynomial C(z'') and a. By using the extended sequence (4.2), it is possible to obtain much better disturbance identification results and more accurate process variation prediction than by identifying process disturbance dynamics for each zone separately. Notice that the polynomial C(z_1) and the sequence e appear in a bi-linear way in equation (4.3), so a suitable identification method solving the nonlinear parameter estimation problem should be used. In this case a Recursive Extended Least Squares (RELS) identification method [9] is used to estimate simultaneously the white noise sequence e and the coefficients of the polynomial C(z~x) in (4.1b). The following outlines the basic RELS algorithm. 4.3 Identification algorithm In accordance with (4.1b), the stochastic process in equation (4.3) can be represented by the following difference equation v = [^( l ) , . . .^ 1 ( /«)^ 2 ( l ) , . . . ,^( / W ) , . . . ,^( l ) , . . . ,^(m)f (4.2) and used to fit the following model: (4.3) v(r)-v(f-l) = e(t)+cle(t-\) + ...+c,e(t-r), (4.4) 40 Chapter 4: Identification of the Disturbance Model where v(t), t = I, 2, ..., mn, is the disturbance sequence given in (4.2), e(t) is the white noise sequence to be estimated as in (4.2) and Cj, c2, Cu are the coefficients of the polynomial C(z~') in the disturbance model (4.3). The equation (4.4) can be represented in the regression form as Av(t) = xT(t)0+e(t), (4.5) where the differentiated disturbance sequence is given by Av(t)=v(t)-v(t-I), (4.6) the parameter vector to be estimated is 0 = [ c 1 , c 2 , . . . , c , ] r , . (4.7) and the regressor vector is x(0 = [e(t- 1),e(t -2),...,e(t~ Of (4-8) The difficulty of the identification problem (4.5)-(4.8) is related to the fact that the variables e(t-l), ... , e(t-l) in the regressor vector x(t) (4.8) are unknown and have to be estimated jointly with the vector 6 (4.7). In this work, such estimation is performed by the Recursive Extended Least Squares (RELS) algorithm. The RELS algorithm can be derived as follows. Assume that for a particular t, the regressor vector x(t) (4.8) is known. Define the one-step ahead prediction error as e(t)=Av(t)-xT(t)0(t-l), (4.9) where 9{t) is an estimate of the vector 0 (4.7) at the time t. 41 Chapter 4: Identification of the Disturbance Model By replacing e(t-l), ... ,e(t-l) in (4.8) with the prediction error, we can obtain the estimate of x(t) at any time: x<t)=[e(t-l),e(t-2),...,e(t-l)f (4.10) Using (4.10) as the regressor vector in a recursive Least Squares estimation scheme instead of (4.8) and propagating the estimate of the residuals eft) (4.9) and a least squares estimate of parameter vector 6 (4.7) forward in time yields RELS algorithm, which has the form where P(t) is the covariance matrix of the recursive least squares estimator of 6 and K(t) is the prediction error gain vector at step t. This method, also known as Pseudolinear Regressions (PLR), combine the estimation of the parameter vector and unobserved components in the regressor. The initial conditions of the algorithm (4.11) are set as follows: the parameter vector 6= 0; the covariance matrix P = cl, where c is a large positive constant and / is identity matrix. The correct selection of the order of the polynomial C(z_1) (4.1b) depends on the characteristics of the disturbances present on the process. Normally assuming the order of C(z~') to be between 1 and 10 is adequate to capture the frequency contents of interest in the disturbance signal. Validation of the developed disturbance identification algorithms is presented in Chapter 6. 6(t) = 6(t - \)+K(t)[Av(t)-xT(t)9(t -1)], K(t)=P(t)x(t) I [1 + xT{t)P(t - l)xT(t)], (4.11) 42 Chapter 5: Tuning of Dahlin Controller and Control Filter Chapter 5 Tuning of Dahlin Controller and Control Filter Having obtained the process time-response model and disturbance model through the methods described in the Chapter 3 and Chapter 4, the tuning problem to be considered in this Chapter is to determine the closed-loop time constant a of the Dahlin controller and the control filter factor ft based on those models. In this Chapter, the tuning strategy based on minimization of a given LQG performance index is introduced. In Section 5.1 the procedure for evaluation of the performance index is described. Section 5.2 gives a straightforward method for minimization of the performance index. Section 5.3 presents a modified tuning scheme that targets on minimizing the worst impact from a set of disturbances. Finally, the developed feedback controller tuning tool and its usage are described in Section 5.4. 5.1 Evaluation of the performance index Assume that the process time-response model (2.2), the controller and filter parameters (2.3 - 2.6), and the stochastic disturbance model (4.1) are known. Then the tuning problem is to determine the closed-loop time constant a of the Dahlin controller and the control filter factor /? which minimize the following performance index: J(a,B) =E{y2 + pAu2} -> min, (5.1) where a is the desired closed-loop time constant of the Dahlin controller; (3 is the control filter factor; E(.) denotes mathematical expectation of a random variable; E(y2) is the process output variance and E(Au2) is the increment control move variance; p is the control weighting factor that penalizes excessive actuator moves. The performance index (5.1) can be evaluated in a straightforward way once the closed-loop 43 Chapter 5: Tuning of Dahlin Controller and Control Filter transfer functions from the white noise e(t) in (4.1) to the loop output y(t) and the actuator moves Au(t) are formed. Evaluation of the performance index (5.1) can be divided into two steps (1) Evaluation of E(y2) (2) Evaluation of E( Au2) The procedures for computing (1) and (2) are similar. We will explain (1) in detail and computations for (2) can be performed in completely similar manner. From the closed-loop system block diagram shown in Figure 2.2, we can derive the transfer function relating the process outputy(t) to the white noise input eft): where g(z~l) is the process model given by (2.2), D(z~') is the Dahlin controller transfer function given by (2.3), Ffz'1) is the control filter transfer function in (2.5), Mfz'1) is the display filter transfer function in (2.6), and Nfz'1) is the disturbance model in (4.1a); A(z~') and Bfz'1) are polynomials with real coefficients dependent on the coefficients in g(z), D(z), Ffz) and N(z). K J ° 1 " (5.3) A(z 1) = a0+a1z l+...+anz " Since eft) is the unit variance white noise, the evaluation of E(y2) can be done according to [16] as 44 Chapter 5: Tuning of Dahlin Controller and Control Filter E(y2) = j§Hy(z)Hy(z-l)z-'dz, (5.4) where <j> denotes the integral along the unit circle in the complex plane, computed in the counterclockwise direction. The integral in (5.4) can be presented in the form A efficient numerical method is given in [16] for computing (5.5). This method requires to build the following table ax .. b0 Z>, .. • K-a n an-x • .. ax ao .. ax i n - l t « - l b0 bx . l n-l A- l an-\ .or _M-1 _M-1 an-\ an-2-.a0 K bl < a\ a\ (5.6) a0 bQ The first row in the table (5.6) is given by the polynomial coefficients in (5.3). The first half of each even row is obtained by writing the coefficients of the first half of the preceding row in reverse order, and the second half of each even row is same as the first half row. The coefficients of the odd rows are obtained using the following formula £-1 k k k i k ai =a*-/ > Y k = al I ao (5.7) 45 Chapter 5: Tuning of Dahlin Controller and Control Filter where a" = at and b" = bi,for / = 0,1,...,n. Once the table (5.6) is built, the integral (5.5) can be found as (5.8) The numerical procedure for computing an integral of the form (5.4) was implemented as a MATLAB program EVALOSS.M. Accordingly, the steps for evaluating E(Au2) are outlined as below: 1) Based on the closed-loop block diagram in Figure 2.2, find the transfer function Hu(z) relating Au to where A = 1-z'1 and Au(z) and Bu(z) are polynomials with real coefficients computed from (5.9), (2.2), (2.3), (2.5) and (4.1a). 2) Make a table similar to (5.6) with the coefficients of the polynomials Bu (z-1) and Au (z-1) in (5.9), and compute the required table coefficients by using formula (5.7). 3) Substitute the coefficients of ,bf obtained from step 2) into formula (5.8) to compute E(Au2) . 5.2 Minimization of the performance index Having obtained the formulae for evaluation of E(y2) and E(Au2) it is straightforward to compute a value of the performance index (5.1) against a chosen pair of the tuning parameters a and B. The optimal the white noise input eft) in (4.3) Bu(z) AD(z)F(z)N(z) (5.9) Au(z) l + g(z)D(z)F(zY 46 Chapter 5: Tuning of Dahlin Controller and Control Filter values of a and 8 can be obtained by minimization of the performance index (5.1). A direct global search method is used in our rrdnimization scheme. It includes computing the values of the performance index (5.1) for different combinations of the tuning parameters a and B, locating the minimum value of the performance index, and finding the corresponding optimal values for a and 8. P 0 IK (1, 1) 1 a Figure 5.1: A search domain with arbitrary grid size for finding optimal tuning parameters a and B The above method requires determination of a two dimensional search domain projected by a and B. Notice that the desired closed-loop time constant Ta (in seconds) is mapped into the parameter a (2.4*), where a e[0, 1] and admissible values for the control filter factor /? in (2.5) are also within [0, 1]. Therefore we chose a unit square with vertex coordinates (0, 0), (0, 1), (1, 1) and (1, 0) as the search domain (Figure 5.1). A grid node in the search domain in Figure (5.1) represents a combination of a and /? tested in the search. Increasing the number of tested combinations of a and B improves the accuracy and increases computational load of the global search. Usually 10 to 20 values of the parameters for each of a and P have to be tested to reach a reasonable vicinity of the minimum of the performance index (5.1). Only approximate optimal values for a and /? can be obtained in the described scheme, more accurate optimal values can be obtained by refining these values using an iterative optimization technique. In our numerical experiments using both mill data and simulated model data, it is found that the performance index surface with respect to a and P is very flat close to the minimum. This suggests that from practical control performance viewpoint the approximate optimal tuning parameters may work as well as the more accurate 4 7 Chapter 5: Tuning of Dahlin Controller and Control Filter ones. Thus it is not necessary to sacrifice computational speed of the tuning algorithms in order to pursue absolutely accurate values of the optimal tuning parameters. 5.3 Performance index modification for robust tuning As mentioned earlier, the disturbance identification uses a particular realization of the stochastic disturbance. Some characteristics of the disturbance may not be captured by the resulting model. Besides, the characteristics of the disturbance may vary with time and with operating conditions. The uncertainties in disturbance model might lead to sub-optimal control performance if not taken into account in the controller tuning. This section presents a modified tuning scheme where the controller and filter are tuned in such a way that a performance lower bound is guaranteed for a wide class of disturbance characteristics. Borrowing the idea of Paganini [10] we define a disturbance set and a worst case gain that characterizes quantitatively the impact of the worst disturbance in the set on the controlled variable, then a new performance index is formed by incorporating the worst case gain into the old one, and again the optimal tuning parameters are obtained by searching for the rninimum of the performance index. Note that the achieved tuning performance is guaranteed not for a particular disturbance but for a set of disturbances. 5.3.1 Notation and assumption The feedback control loop in Honeywell-Measurex CD control system shown in Figure 2.2 is redrawn in Figure 5.2 with the white noise shown as a loop input. The transfer function relating the white noise to the process output was obtained in (5.2) and can be written in the form of an impulse response model OO (5.10) t=o 48 Chapter 5: Tuning of Dahlin Controller and Control Filter where z is the forward time shift operator and h(t) for t = 0, 1, ... are the closed-loop impulse response coefficients. Figure 5.2: A block diagram of feedback control loop in Measurex CD control system Define the autocorrelation function corresponding to the transfer function HY in (5.10) as ry(T): = Zh(t + T)h(t). t=o (5.11) In (5.11) it is assumed that h{t) and r (r) e so that the series (5.11) ry(r) converges for any t and 5.3.2 Definition of a disturbance set Consider a realization of the process variation as a sequence e(t) with length of N and its correlation function: 49 Chapter 5: Tuning of Dahlin Controller and Control Filter r.(r)=5>(f + *X0, T = 0,...,N-1 (5.12) f=0 \rM'r.(0)\*r, r = \,2,...,N-\ (5.12a) Figure 5.3 shows the normalized correlation function re (r) / re (0) of a disturbance that belongs to a disturbance set specified with a parameter The disturbances whose correlation function are within y bound are considered as the set members. 1 Autocorrelation function of the disturbance © I I | j 0.5 Y 0 ~y -0.5 ( • ) 5 10 15 X 20 25 30 Figure 5.3: Gamma specifies a disturbance set Definition: The set of signals of length N which are white in the time domain sense with accuracy y, up to lag T, is defined by WN,rJ: = {ezRN:\re(T)\<r.re(0), r = \,...,T) (5.13) The parameter y indicates the degree to which the disturbances in the set WN T ^ c l ° s e t o white noise; y=0 means that the set contains only pure white noise; y = 1 suggests that the signals in the set can 50 Chapter 5: Tuning of Dahlin Controller and Control Filter be of various forms (e.g., steps, sinusoidal and stochastic signals); 0 < y< 1 indicates that the signals are some versions of white noise (filtered white noise, etc.). In order to express the worst impact of a set of disturbances on the system output we now define the worst gain of the system under signals in WN<rj: | | ^ I L w : = 8 u p { M : ^ = J5r y e,eeFr w , | |e |kO> (5.14) where || y || denotes the Euclidean norm for a finite sequence y(t). According to [10], the following theorem holds Theorem: Let Hybe a LTI, SISO system with h(t) e /,. If e(t) e WN.r.T ^ s 3 1 1 input signal to the system with the transfer function Hy and y = Hye is the corresponding steady-state output, then \\Hy\\l<\\Hy\\iNrT<\\Hy\\\(\-y)+y 2>,(r)|, (5-15) r = - J where \\Hy\\22 istheH 2 norm of the system transfer function (5.2); 1 Hy \\Wji r is the squared worst gain of Hy under signals in WNyT; y is the set-defining parameter as in (5.13); ry(r)is the system autocorrelation function defined by (5.11). The parameter y plays an important role in parameterizing the freedom allowed in the disturbance signal model and results in a worst-case gain which varies between the H 2 norm for y= 0 and the V\x norm for y= 1. 5.3.3 A robust tuning performance index Having obtained a useful expression (5.15) for the worst-case gain, we can now form a new tuning performance index that targets minimization of the worst impact under a set of disturbances on the process 51 Chapter 5: Tuning of Dahlin Controller and Control Filter output. This performance index is obtained by modifying the index (5.1) in accordance with the guaranteed variation estimation (5.15) and has the form j=(i-r)\\Hy\\l+r hry(T)\+p{(l-r)\\HJ\l+rhru(T)\}, (5-i6) r=-r T=-T where Hu denotes the closed-loop transfer function given in (5.9) relating the increment control action to the white noise; \\Hu\\l is H2 norm of Hu; ru(t) is the autocorrelation function of Hu as in (5.11); the control weighting factor p is used to penalize excessive control moves. In the performance index (5.16) y= 0 and p = 0 gives a H2 optimal problem, and y = 1 and p = 0 results in a H«, optimal problem, while y = 0 corresponds to a LQG problem. When y and p are chosen properly, the tuning based on (5.16) is robust. For our tuning scheme, y is computed in the process of the disturbance identification of Chapter 4 by using the following formula: y =Max (rjr)) U(0)J forr=i, 2, ...,T, (5.17) where re (r) is the autocorrelation function of the prewhitened process disturbance given in (5.12) and T is the maximum lag of the autocorrelation function and is set 30 as default value. 5.4 Feedback controller tuning tool Based on the above described tuning methods, a feedback controller tuning tool FCTune was designed for tuning the Dahlin controller and control filter in Honeywell-Measurex CD control system. The tool is written in MATLAB and some of its key algorithms have been compiled into MATLAB MEX files for fast execution of the software. One can type in FCTune in the command line in MATLAB window to invoke the main tuning screen of the tool. 52 Chapter 5: Tuning of Dahlin Controller and Control Filter The main tuning screen shown in Figure 5.4 is divided into two sections. The upper section allows to perform the disturbance identification. One can adjust the disturbance model order and observe the identification results shown on another screen. The value of GAMMA on the upper-right corner of the tuning screen shows parameter y (5.12a), which indicates closeness of the identification residual to white noise. This value, together with the autocorrelation function (5.11) of the identification residual helps to determine a suitable disturbance model order. The field # of ZONE displays number of valid control zones. The lower section of the tuning screen is designed for tuning the closed-loop time constant Ta and the control filter factor /?. One may follow the following procedures for tuning Ta and p. • Enter the current controller parameters into the CURRENT PARAMETERS column: closed-loop time constant Ta, process gain kp, process time delay Td, control time constant TT and control filter factor /?. The predicted performance for the current controller, represented by the profile 2 a and actuator 2 cr, is shown below the CURRENT PARAMETERS column on the tuning screen. • The suggested tuning parameters Ta and P depend on the value of CONTROL PENALTY (the control weighting factor p in the performance index (5.1)). For the initial value p = 0.1 of CONTROL PENALTY, the suggested tuning parameter Ta, kp, Td, TT and P are shown in the SUGGESTED PARAMETERS column on the tuning screen. The predicted profile 2 a and actuator 2 cr for the suggested controller are also displayed on the tuning screen. • A user can adjust the CONTROL PENALTY parameter and compare the predicted performance of the suggested controller with that of the current controller. Keep doing this till the predicted performance trade off for the suggested controller is superior to that for the current controller. • Implement the suggested tuning parameters into the control system. 53 Chapter 5: Tuning of Dahlin Controller and Control Filter The accuracy of the predicted profile 2 cand actuator 2 crfor a particular controller can be verified by collecting and processing steady state data obtained with this operating controller. 54 Chapter 6: Validation of the Tuning Algorithms Chapter 6 Validation of the Tuning Algorithms This chapter is devoted to validation of algorithms and methods for tuning CD controller feedback loop as described in the Chapter 3 to 5. The key parts of the developed algorithms, the performance of identification of the disturbance model, as well as prediction of the process variation Ety2) and incremental control variation E(Au2) critically influence practical applicability of the algorithms, therefore these parts need careful testing and validation on real-life data. Simulation models are used mainly for software verification, while mill data are used for checking the process variation prediction and proving the proposed tuning concept in the real-life conditions. 6.1 Validation Approach Before proceeding to validation of the disturbance identification and process variation prediction, some useful references and formulae are given. Calculation of variation using the process data For validation of predicted process variation, the closed-loop variation need to be computed from the steady state measurements. Traditionally the total process variance is divided into three independent components: a CD-only, a MD-only and a residual. °'2=C!'cD+a'MD+aR- (61) The total 2 cr, CD 2 cr, MD 2cr and residual 2 a can be computed using an available method, for example, the TAPPI recommended formula [14]. 55 Chapter 6: Validation of the Tuning Algorithms Prediction of variation using the model Based on the identified process and disturbance models, the process variance and incremental actuator variance can be predicted. From the closed-loop system block diagram shown in Figure 2.2, one can derive the transfer function relating the process output y(t) to the white noise input eft): where g(z~') is the process model given by (2.2); D(z1) is the Dahlin controller transfer function given by (2.3); Ffz'1) is the control filter transfer function in (2.5); Nfz'1) is the disturbance model in (4.1a) and M(z'2) is the display filter transfer function in (2.6). In a similar way, the transfer function relating the incremental control move uft) to the white noise input eft) can be obtained II A D ( z - 1 ) F ( z - 1 ) A T ( z - ) \ +g{z-x)D{z-x)F{z-') where A = 1-z'1. The process and actuator variance predictionE(y2) and E(Au2) can be evaluated according to [16] E(y2)=-§Hy(z)Hy(z-')z-*dz, (6.4) 56 Chapter 6: Validation of the Tuning Algorithms (6.5) where <J> denotes the integral along the unit circle in the complex plane, computed in the positive direction. According to variance partition, (6.4) represents MD and residual variance of the process output and (6.5) denotes MD and residual variance of the incremental actuator move. Estimation of variation spectrum using the data MATLAB function SPECTRUM.M is used to estimate the process and actuator spectrum from steady state data. The spectrum can be used as a benchmark to check the predicted process and actuator spectrum, which is based on the identified process and disturbance models. Prediction of variation spectrum using the model The disturbance model (4.3) can be written as where A£ (z) is Z transformation of the incremental process variation sequence, e (z) is Z transformation of the white noise sequence. The power spectrum of the predicted incremental process variation is A£(z) = C(Z-1)<T, (6.6) e(z) (6.7) The closed-loop predicted process power spectrum can be obtained directly from (6.2) 57 Chapter 6: Validation of the Tuning Algorithms In \Hy(e~n\2 (6.8) The predicted incremental actuator power spectrum can be obtained from (6.3) (6.9) 6.2 Simulation Validation As mentioned earlier, the feedback controller tuning scheme crucially depends on quality of prediction of the process and actuator variation. This requires to have a good process model and especially a fairly accurate disturbance model. Therefore it is necessary and important to validate the disturbance identification and prediction of the process and actuator variation. Example 1: Verification of the disturbance identification algorithm The block diagram for the disturbance identification for a simulated process is shown in Figure 6.1. The parameters of the simulated first order process model and the integrated moving average disturbance model are given in Table 6.1. The number of control zones is set to 1. For simplicity the display filter factor is set to 1 (no filtering). The standard deviation of the white noise input e is 0.1, while the excitation signal u is a unit step. All simulated sequences, e, u, and y, have 100 samples. The scan time is 30 seconds. Figure 6.2 shows the process model identification results. In Table 1 the estimated parameters of the first order model are compared with the nominal model parameters which were used in simulation. The estimation error is due to the correlated disturbance (random walk). 58 Chapter 6: Validation of the Tuning Algorithms Noise model u Process model + + Identification algorithm Display filter Figure 6.1: Identification of disturbance model for a simulated process Gain Delay Time constant Disturbance standard kp Td Tr model (9) deviation cr Simulation 1 50 sec. 130 sec. l-0.5z_1 0.1 model 1-z"1 Identified 0.913 63.34 sec. 131.97 sec. 1-0.41Z"1 0.1109 model 1-z"1 Table 6.1: Comparison of the simulation model and estimated model For check of the disturbance model identification accuracy, the autocorrelation function of the identification residual is plotted against that of the disturbance in Figure 6.3. The former is much closer to the white noise autocorrelation function than the latter. This indicates that a 59 Chapter 6: Validation of the Tuning Algorith 1.2 Model:g=0.913; Tdel=63.34; Trise= 131.9660 1 1 i -1 0.8 0.6 0.4 0.2 A U -0.2 --0.4 ( 1 I 3 500 1000 Time in Seconds; 1500 MD response; 2000 2500 — estimated response 3000 Figure 6.2: The process model identification using simulated model Figure 6.3: Disturbance model identification using simulated model 60 Chapter 6: Validation of the Tuning Algorithms plausible disturbance model is obtained. Table 1 shows the estimates for the noise model. In general for the simulated sequence of the length of 100 samples, the estimation error for the noise model coefficient is within 7-18 %. Example 2: Optimalitv of the Dahlin controller Dumont [4] proved that if the process is a first order with dead time and the process noise is described by a first-order integrated moving-average (IMA) process, then the Dahlin controller is a rninimum variance controller provided that the closed-loop time constant a is set to Cj (the only coefficient of the polynomial C(z~l) (4.1b) in the noise model). In the following example, this analytical result is used to verify the developed tuning algorithm. A feedback control loop of a traditional CD control system, same as one shown in Figure 2.2, was simulated The parameters of the transfer function for each component in the feedback control loop are given in Table 6.2. The scan time is set to 30 seconds. For verification of the Dahlin Trans, function g(z'1) - process model D(z~l) - Dahlin controller Parameter kp Td Ta kp Td Tt Value 1 60 sec. 120 sec. * 1 60 sec. 120 sec. Trans, function F(z"1) - control filter - display filter Nfz'1) - noise model Parameter fi 9 Cl a Value 1 1 0.75 1 Table 6.2: Parameter values of each component transfer function used in the example 2 tuning constant a using the above results, the control filter and the display filter in the feedback loop are disabled (their filter factors are set to 1). The process model is assumed to match the real 61 Chapter 6: Validation of the Tuning Algorithms plant exactly. By continuous to discrete parameter conversion in (2.2*) and (2.4*), the considered process is described by Process variance vs. closed-loop lime constant 2.51 1 1 1 1 1 1 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Optimal alpha=0.750; P_var = 1.12 Figure 6.4: Optimality of the Dahlin controller 0.22z~3 1-0.75Z"1 y= u + -1-0.78Z"1 1-z"1 (6.10) The Dahlin controller for the process (6.10) is given by 0.55(1 -aXl-O.VSz-1) V(z ) = l-az'1 -(l-a)z -3 ' (6.11) where a defines the closed-loop time constant Ta in (2.4*). The tuning performance index (5.1), where the control weighting factor p was set to 0 for purpose of verification of the tuning algorithms, was computed in accordance with (6.4) for different values of the tuning parameter a. Figure 6.4 shows the performance index versus the 62 Chapter 6: Validation of the Tuning Algorithms tuning parameter a. It can be seen in Figure 6.4 that the optimal a is 0.75 that is the same as the noise model parameter Ci. This confirms the result of [4] and also verifies the consistency of the tuning algorithm and software. 6.3 Paper M i l l Data Validation 6.3.1 Case 1: Steambox actuator, moisture process In this section, a set of bump test data and a set of closed-loop steady state data from a paper mill in Quebec are used for validation of the developed disturbance identification and process variation prediction algorithm. The bump test data are used to identify the process and disturbance models. This allows to predict the process and actuator variation which is expected in closed-loop for the identified models and a given feedback controller. The prediction is then compared against the measured variation values which are computed for actual data logged during a steady state operation of paper machine with the given controller. The startup data were collected on a paper machine producing 48.8 gr. newsprint paper. The moisture content is controlled by a 70 zones (101.6 mm) Devronizer. To do initial alignment and tuning after the Devronizer installation, a moisture bump test was performed and the feedback controller parameters were determined with the HMDI software ProfileView (a process identification tool other than what is described here). The implemented controller parameters are shown in Table 6.3. After the new controller parameters were implemented and the process reached steady operation, 26 scans of closed-loop operational data were logged. Model identification Using the collected bump test data, the process time response model was identified using the method developed in Chapter 3. The estimated process parameters are shown in Figure 6.5. Figure 63 Chapter 6: Validation of the Tuning Algorithms 6.5 also shows the fitting curves of the process time response and estimated time response, which indicate that the process identification quality is satisfactory. The process disturbance model was chosen according to (4.1), the order of C(z'J) is selected as 1 and the number of valid actuator zones is 67, which determines how many zone disturbance sequences are put into a series and used for the disturbance model identification. The identified disturbance model is as follows , l v l-0.87z-1 JVYz_1) = — , a = 0.539, y = 0.1033, (6.12) 1 — z where cr is standard deviation of the residual and represents the process variation intensity; Y (0 < y < 1) > defined in (5.12), characterizes the disturbance identification accuracy. 1.6 Process parameters: Gain = -0.056; Tdel= 28.54; Trise = 181.0376 i i i f l • i i r 1.4 1.2 1 0.8 // // 0.6 // // 0.4 // 0.2 0 -0.2 ( -\ 1 1 • i • i i i 100 200 300 400 500 600 700 800 900 1000 time in sec; MD response - solid; estimated response - dashed Figure 6.5: Process time response model identification 64 Chapter 6: Validation of the Tuning Algorithms Controller Gain Time delay Ctrl time alpha DAT filter parameter Td constant TT Ta factor p value -0.046 43.6 104.13 150 0.2 Table 6.3: Implemented controller parameters Figure 6.6 shows the identified disturbance model parameters and the autocorrelation functions of the process variation and the noise identification residual. In Figure 6.6, a small y (0.103) shows that the noise identification residual is close to a white noise sequence and the estimated disturbance model has good applicability. Check of predicted process and actuator variation The disturbance model (6.12) was validated by checking how good it is for predicting the closed-loop process and actuator variations. As comparative references, the "actual" process 2 cr and incremental control 2 crwere computed using 26 scans of logged steady state data. The Covf. of disturbance - Solid; Covf. of residual - Dashed 1 1 i 1 i — 1 i — • ; • ; • i i i ' 0 S 10 15 20 25 30 C(1)= 1.00;C(2)=-0.87; Figure 6.6: Disturbance model identification 65 Chapter 6: Validation of the Tuning Algorithms predicted process 2crand incremental control 2 crwere computed based on above identified process model in Figure 6.5, disturbance models (6.12) and used controller in Table 6.3. Then the predicted process and incremental control 2 a were compared against the actual variations to validate the identified models and the tuning scheme. The process 2 a and incremental actuator 2 cr were predicted using equations (6.4) and (6.5) based on the above identified process model (Figure 6.5), disturbance models (6.12) and used controller (Table 6.3). In order to check the prediction, the actual process and actuator 2 cr were computed using the logged steady state data [14]. In Table 6.3, the predicted process/actuator 2 cr are compared with the actual process/actuator 2 cr. Note that both predicted and actual process 2 cr consists of two independent terms: MD 2 cr and residual 2 cr. In Table 6.4, the prediction of incremental actuator 2 a has 120% error (3.57 against 1.604), while the prediction of process 2 cr has 94% error (0.4918 vs. 0.253). The probable reasons for relatively large prediction error are the short sequence of the available steady state data (26 scans). Profile 2-sigma Actuator 2-sigma Measured value* 0.253 1.604 Predicted value** 0.4918 3.57 * computed from the steady state data with TAPPI recommended \ ** evaluated with (6.4) and (6.5) formula [14]. Table 6.4: Comparison between measured process/actuator variation and predicted process /actuator variation Check of predicted process and actuator spectrum The identified process model and disturbance model were validated by checking accuracy of the power spectrum prediction for the closed-loop operation. The "actual" process power spectrum was computed using the logged closed-loop steady state data, and the predicted process power 66 Chapter 6: Validation of the Tuning Algorithms spectrum was computed based on the identified models and used controller. The predicted spectrum was compared with the measurement spectrum so that the identified models were validated. Open-loop process spectrum check For the identified disturbance model (6.12), the predicted incremental process spectrum can be obtained using equation (6.7) 0.5392 In |(l-0.87e-y<a)|2 (6.13) Figure 6.7: Spectral check for the disturbance model: Open-loop spectrum for the model - solid; spectrum computed directly from the data - dashed The measured process spectrum was computed from the open-loop disturbance sequences using MATLAB function SPECTRUM.M, and the predicted process spectrum (6.12) was computed using MATLAB function DBODE.M. Here again the zone disturbance sequences are cascaded in 67 Chapter 6: Validation of the Tuning Algorithms order to increase data available for the spectrum estimation. The predicted process spectrum (6.12) is compared with the measured process spectrum in Figure 6.7. Figure 6.7 shows that the predicted spectrum matches the measurement spectrum fairly well over low frequency range. This indicates that the obtained disturbance model has good applicability. Closed-loop process spectrum check The closed-loop transfer function relating the process output y(t) to the white noise input eft) was obtained according to (6.2). The parameter values in each component of the transfer function are listed in Table 6.5. These values are obtained from Figure 6.5, Table 6.3 and (6.12), based on the above mentioned process and disturbance identification results. Closed-loop process spectrum validation was performed for the controller with parameters shown in Table 6.3. Trans, function g(z~l) - process model D(z~l) - Dahlin controller Parameter k, Td Ta K Td TT Value -0.056 28.54 181.0 150 -0.046 43.6 104.13 Trans, function F(z']) - control filter Mfz1) -display filter N(z~l) - noise model Parameter P 9 C l Value 0.2 0.2 -0.87 0.539 Table 6.5: Parameter values of transfer functions used in the spectrum prediction validation The closed-loop predicted process power spectrum was obtained by substituting the parameter values in Table 6.5 into (6.8). The predicted process spectrum is plotted in Figure 6.8 (solid line). In order to validate the predicted process spectrum, MATLAB function SPECTRUM.M was used 68 Chapter 6: Validation of the Tuning Algorithms to compute the measurement spectrum using the steady state operational data. For comparison, the measurement spectrum is also plotted against the predicted spectrum in Figure 6.8 (dashed line). Figure 6.8 shows that although the measured spectrum is below the predicted spectrum, the shapes of the two spectrum curves are similar. This means that although there are some errors in the prediction of the signal's total energy, the relative energy for each frequency component is predicted correctly. . Closed-loop model spectrum- solid; spectrum estimate - dashed 10P c 1 1 < 1 1 1 1 1 1 1 1 F F J| 1 , , 1 1 L _ _ 1 , 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 ' Frequency rad/sec Figure 6.8: Closed-loop process spectral check: Process spectrum predicted using the closed-loop model - solid; spectrum computed directly from the closed-loop data - dashed Closed-loop actuator spectrum check The predicted actuator spectrum was obtained by substituting the parameter values in Table 6.5 into equation (6.8), and the measured actuator spectrum was calculated using MATLAB function SPECTRUM.M based on the steady state data. The predicted actuator spectrum and the measured actuator spectrum are plotted in Figure 6.9 in order to verify applicability of the identified process and disturbance models. Mismatch between these two is clearly seen. Since the shapes of these two 69 Chapter 6: Validation of the Tuning Algorithms curves are similar, the relative energy for each frequency component in the signal is predicted correctly. Figure 6.8, Figure 6.9 and Table 6.4 all show that the predicted process and actuator variations are larger than the measured process and actuator variations. , Closed-loop model spectrum- solid; spectrum estimate - dashed 102 r . 1 1 1 . 1 1 1 1 , 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 Frequency rad/sec Figure 6.9: Closed-loop actuator spectral check: Actuator spectrum predicted using the closed-loop model - solid; spectrum computed directly from the closed-loop data - dashed There are many factors affecting accuracy of the spectrum prediction and the 2 cr prediction mentioned above, such as insufficient and/or inaccurate models, changing disturbance dynamics, actuation nonlinearity, short data sequences used for the model identification or the statistic calculations and so on. In this case, since only 26 scans of steady state data were logged and used for computing the "actual" process spectrum and 2 cr, the resulting spectrum and 2 cr may not represent the characteristics of the process variation. 70 Chapter 6: Validation of the Tuning Algorithms 6.3.2 Case 2: Slice lip actuator, weight process In order to validate the developed disturbance identification and process variation prediction algorithm, a mill trial was conducted in a BC paper mill on July 17, 1997. The paper machine produces 40.5 g. newsprint. The basis weight is controlled by motorized slice lip actuators and moisture by water spray actuators. The average scan time is 26 seconds. Model identification During the mill trial, a weight bump test was performed and 38 scan data were collected. The developed software was used to process the bump test data and identified the process time response model. The identified process parameters together with the model fitting curve are shown in Figure 6.10. 1 . 4 M o d e l : g = 0 . 2 0 4 ; T d e l = 5 6 . 6 7 ; T r i s e = 1 3 . 3 5 9 8 • I I I T I I I I 1 . 2 1 tfi \T./v" \~~ "7™ \7 ~ 0 . 8 0 . 6 -0 . 4 -0 . 2 - A J 0 - /\ - 0 . 2 - 0 . 4 ( • • i i • i • i i ) 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0 0 T i m e i n S e c ; p r o c e s s t i m e r e s p o n s e ; - e s t i m a t e d r e s p o n s e Figure 6.10: Identification of the process time-response model An integrated-moving-average model in the form of (4.1) was used for the process disturbance identification. The degree of the noise shaping filter C(z_1) in (4.1b) is chosen as 1 and the 71 Chapter 6: Validation of the Tuning Algorithms number of the valid control zones is 68. Figure 6.11 shows the identified parameters and the autocorrelation functions of the process variation and the residual. As mentioned earlier, y illustrates how close the autocorrelation of the residual (dashed line) is to that of an white noise. A small obtained y (0.098) indicates that the identification result is credible. Since the original disturbance is close to an white noise (see Figure 6.11), it is adequate to choose the degree of C(z_1) as 1. The identified disturbance model is as follows: JV(z~1)=1"a8^ <r=0.1709, ^ = 0:098 (6.14) 1 — z where cr is standard deviation of the residual and represents the process variation intensity. 1 0.8 0.6 0.4 0.2 0 -0.2 C Covf. of disturbance - Solid; Covf. of residual - Dashed • i : • i ) 5 10 15 20 25 30 C(1)= 1.00; C(2)=-0.81; Figure 6.11: Identification of the process disturbance model Check of predicted process and actuator variation In order to validate the above identified models, two different controllers with parameters shown in Table 6.6 were implemented and 55 scans of steady state closed-loop data for each of them were logged. Controller #1 is used at the mill, the controller #2 was implemented during the trial with the 72 Chapter 6: Validation of the Tuning Algorithms purpose of the prediction verification. In Table 6.6, Ta represents the continuous desired closed-loop time constant of the Dahlin controller and /? stands for the control filter factor (0 < B < 1, 1 -no filtering). The filtering of the control filter for the controller #2 was smaller than that for the controller #1 by more than 100%, so more aggressive actuator moves for the controller #2 could be expected (see Table 6.8). The prediction of the process and actuator variation based on the identified models and the used two controllers was evaluated through (6.4) and (6.5). The predicted variation was compared with the actual variation in Table 6.6 and 6.7. Controller one Controller two Process gain kp 0.2174 0.1952 Process time delay Td 85.09 74.77 Control time constant TT 129.35 119.27 Closed-loop time constant Ta 258.70 268.4 Control filter factor B 0.2 0.49 Display filter factor <p 1 1 Table 6.6: Parameters of two tested controllers Table 6.7 shows the measured and predicted 2 sigma of process variation for the two tested controllers. The prediction is very accurate and the error is less than 5%. Table 6.8 shows the predicted 2-sigma and measured 2-sigma for incremental actuator move. Although there were some prediction errors, the change direction of the actuator variation when tuning the controller was predicted correctly. When the controller parameters was changed from the setting one to setting two, the predicted actuator 2a increased from 0.1601 to 0.4072, which was in satisfactory correspondence with the change in real actuator variation ( increased from 0.1657 to 0.4448). The main reasons for the prediction errors are insufficient and/or inaccurate models and varying process and disturbance dynamics. For tractability of the algorithm it is assumed that the process 73 Chapter 6: Validation of the Tuning Algorithms disturbances in adjacent zones are independent, the interaction of the neighboring zone disturbances can cause the prediction errors. These issues need to be addressed in order to improve further the disturbance identification and prediction. Controller one Controller two Predicted profile 2-sigma(6.4) 0.3471 0.3439 Measured profile 2-sigma [14] 0.3421* 0.3323* * consists of MD and residual component Table 6.7: Comparison of predicted 2-sigma and measured 2-sigma for the process variation Controller one Controller two Predicted incremental control move 2-sigma (6.5) 0.1601 0.4072 Measured incremental control move 2-sigma [14] 0.1657 0.4448 Table 6.8: Comparison of predicted 2-sigma and measured 2-sigma for incremental actuator move Check of predicted process and actuator spectrum The above process model and disturbance model identified from the open-loop data have been verified by checking accuracy of the power spectrum prediction for the closed-loop operation. The predicted closed-loop process and actuator spectrum based on the identified models and used controller are obtained with (6.8) and (6.9), while the measured power spectrums are computed from the process steady state measurement data using MATLAB function SPECTRUM.M. The predicted spectrum is compared with the measurement spectrum so that the identified models can be validated. Open-loop process spectrum check 74 Chapter 6: Validation of the Tuning Algorithms For the identified disturbance model (6.14), the predicted incremental process spectrum can be obtained using equation (6.7) 0.17092 2n |(l-0.81e-;<n)|2 0.12 0.1 0.08 0.06 0.04 0.02 -Noise model spectrum- solid; spectrum estimate - dashed 0.06 0.08 Frequency rad/sec 0.12 0.14 (6.15) Figure 6.12: Comparison between predicted spectrum and measurement spectrum of the open-loop process variation / This predicted spectrum can be checked with the measurement spectrum computed from the open-loop disturbance sequences using MATLAB function SPECTRUM.M. Here again the zone disturbance sequences are cascaded in order to increase amount of data used for the spectrum estimation. Figure 6.12 shows that the predicted spectrum matches the measurement spectrum fairly well. This illustrates a good applicability of the identified disturbance model. Since the same data set (the weight bump test data) was used for the disturbance model identification and for the measurement spectrum estimation, this check can only prove the disturbance model is adequate and 75 Chapter 6: Validation of the Tuning Algorithms accurate for this data set. It will be shown in the following section that the identified process and disturbance models can be used for predicting the process variation based on another data set. Closed-loop process spectrum check The closed-loop transfer function relating the process output y(t) to the white noise input eft) is given by (6.2), where the parameter values in each component are listed in Table 6.9. These values are obtained from Figure 6.10, Table 6.6 and (6.15), based on the above mentioned process and disturbance identification results. Controller #1 was used for the closed-loop process spectrum validation. Trans, function gfz'1) - process model Dfz'1) - Dahlin controller Parameter K Td Tr Ta kP Td Value 0.204 56.67 13.36 258.7 0.217 85.09 129.4 Trans, function Ffz'1) - control filter Mfz'1) - display filter Nfz'1) - noise model Parameter fi 9 C l Value 0.2 1 -0.81 0.1709 Table 6.9: Paramer value of each component transfer function used in the spectrum check The closed-loop predicted process power spectrum is obtained by substituting the parameter values in Table 6.9 into (6.8) and is plotted in Figure 6.13 (solid line). MATLAB function SPECTRUM.M was used to compute the measurement spectrum using the steady state operational data. For comparison, the measurement spectrum is also plotted against the predicted spectrum in Figure 6.13 (dashed line). Good match between these two spectra shows not only the process model but also the disturbance model are applicable to the real process variation prediction. Li Figure 6.13 the mismatch over very low frequency range is likely caused by the way in which the 76 Chapter 6: Validation of the Tuning Algorithms 0.07 Model based predicted spectrum - solid; spectrum - dashed i i i i i i 0.06 0.05 / -0.04 J V i / 0.03 0.02 1 I f 0.01 ( • • • i i i ) 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Frequency rad/sec Figure 6.13: Comparison between predicted process spectrum and measurement spectrum 10° Closed-loop model spectrum- solid; spectrum estimate - dashed • i 1 1 1 1 10' 10"2 \\ v\ \ \ V \ < \ : 10° • 1 i i i i 3 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Frequency rad/sec Figure 6.14: Comparison between predicted spectrum and measurement spectrum for actuator variation disturbance sequences are handled. Putting the zone disturbance sequences into one series will generally add some false low frequency component into the signal, so the information over very low 77 Chapter 6: Validation of the Tuning Algorithms frequency range (period > the bump test duration or duration of steady state data collection) should be disregarded. Closed-loop actuator spectrum check The predicted actuator spectrum was obtained by substituting the parameter values in Table 6.9 into equation (6.9). The measurement spectrum was calculated using MATLAB function SPECTRUM.M based on the steady state data. The predicted actuator spectrum and the measurement spectrum are plotted in Figure 6.14. Good match between these two curves suggests that the identified process and disturbance models can be used for predicting the process variation and for tuning the feedback controller. 6.3.3 Case 3: Dilution actuator, weight process In this section, a set of bump test data and a set of closed-loop steady state data collected from a dilution-actuator-controlled basis-weight process in a US paper mill were used for validation of the developed feedback loop tuning algorithms. The controlled process contains 84 dilution actuators and 168 databoxes. Average scan time is 42 seconds. Both the bump test data and the steady-state data contain 20 scans. The bump test data were used to identify the process and disturbance models. The predicted values for process variation and for actuator variation were obtained based on the identified models and implemented feedback controller. The prediction was then compared with the measured variation values which were computed for data logged during a steady state operation of paper machine with the given controller. Model identification 78 Chapter 6: Validation of the Tuning Algorithms The process time response model was identified using the method discussed in Chapter 3. Figure 6.15 shows the estimated process parameters and the fitting curves of the process time response and estimated time response. A small process time constant (10.2 sec.) is in correspondence with the real situation, where the effect of the dilution water on the basis weight is instant. An 3rd order integrated-moving-average model (4.1a) was chosen to describe the process variation. The number of the valid actuator zones is 75. It means that 75 zone disturbance sequences were put into a single data vector and used for the disturbance model identification. The identified disturbance model is as follows: 1.2 Model: g=-0.017; Tdel=106.67; Trise= 10.2042 1 0.8 / / /> // 0.6 0.4 0.2 0 - " -0.2( i t i l l 3 100 200 300 Time in Seconds; 400 500 600 700 MD response; — estimated response 800 Figure 6.15: Process time response model identification for Case 3 , l-0.47z"1-0.36z-2-0.06z"3 N(z~1)= Tl , o- = 0.4585, ^ = 0.1006, 1 — z (6.16) where o"is standard deviation of the identification residual and represents the noise input intensity; Y (0 -1), defined in (5.12a), indicates the quality of the disturbance identification. Figure 6.16 79 Chapter 6: Validation of the Tuning Algorithms shows comparison of the autocorrelation function of the disturbance with that of the identification residual. Since y (0.1006) is small, the identified model is applicable to the process variation prediction. 0.8 0.6 0.4 0.2 -0.2 1 ! \ l\ \\ \\ \.\ \ \ \ \ \ \ — — J ~~ " -» ^ \ ^ —-• 2 4 6 8 10 C(1)= 1.00; C(2)=-0.47; C(3)=-0.36; C(4)=-0.06; 12 Figure 6.16: Disturbance model identification for Case 3 Check of predicted process and actuator variation The process model shown in Figure 6.15, the disturbance model (6.16) and the implemented controller in Table 6.10 were used for predicting the closed-loop process and incremental actuator 2 a. Then the predicted process and incremental actuator 2 a were compared with the actual process and incremental actuator 2 o"to validate the identified models and the tuning scheme. Controller Gain Time delay Ctrl time Alpha Ctrl filter parameter k, Td constant TT Ta factor p value -0.0052 252 189.46 378.91 0.2 Table 6.10: Implemented control er parameters 80 Chapter 6: Validation of the Tuning Algorithms Profile 2-sigma Actuator 2-sigma Measured value* 1.421 20.5 Predicted value** 0.523 27.16 * computed from the steady state data with TAPPI recommended \ ** evaluated with (6.4) and (6.5) brmula [14]. Table 6.11: Comparison between measured process/actuator variation and predicted process/actuator variation for Case 3 The process 2 cr and incremental actuator 2 cr were predicted using equation (6.4) and (6.5) based on the above identified process model (Figure 6.15), disturbance models (6.16) and used controller (Table 6.10). In order to check the prediction accuracy, the actual process and actuator 2 cr were computed using the logged closed-loop steady state data [14]. In Table 6.11 the predicted process/actuator 2 cr are compared with the actual process/actuator 2 cr . In Table 6.11, the prediction of incremental actuator 2 a has 32% error (27.16 against 20.5), while the prediction of process 2 cr has 63% error (0.523 vs. 1.421). The reason for large prediction error is probably due to the short sequence of the available steady state data (20 scans). Check of predicted process and actuator spectrum The identified process model and disturbance model were validated through check of the power spectral prediction. The "actual" process power spectrum was computed using the logged closed-loop steady state data, and the prediction of the process power spectrum was performed based on the identified models and used controller. The predicted spectrum was compared with the measurement spectrum so that the identified models were validated. Open-loop process spectrum check 81 Chapter 6: Validation of the Tuning Algorithms For the identified disturbance model (6.16), the predicted incremental process spectrum can be obtained using equation (6.7) 0.45852 2n •|(1 - 0.47e-'* - 0.36e-2jVD - 0.06e-3j<D)| -3ja> \|2 (6.17) 0.5 N o i s e m o d e l s p e c t r u m - so l id ; s p e c t r u m e s t i m a t e - d a s h e d 0.01 0.02 0.03 0.04 0.05 Frequency rad/sec 0.06 0.07 0.08 Figure 6.17: Spectral check for the disturbance model The measured process spectrum was computed from the open-loop disturbance sequences using MATLAB function SPECTRUM.M, and the predicted process spectrum (6.17) was computed using MATLAB function DBODE.M. Again the zone disturbance sequences were cascaded in order to increase data available for the spectrum estimation. The predicted process spectrum (6.17) is compared with the measured process spectrum in Figure 6.17. Figure 6.17 shows that the predicted spectrum match the measurement spectrum fairly well. This indicates that the obtained disturbance model has good applicability. Closed-loop process spectrum check 82 Chapter 6: Validation of the Tuning Algorithms The closed-loop transfer function relating the process output y(t) to the white noise input e(t) was obtained according to (6.2). The parameter values in each component of the transfer function are listed in Table 6.12. These values are obtained from Figure 6.15, Table 6.10 and (6.16), based on the above mentioned process and disturbance identification results. The controller with parameters shown in Table 6.10 was used for the closed-loop process spectrum validation. 101 Closed-loop model spectrum- solid; spectrum estimate - dashed 10° \ \ 10"' \ N» \ v» 10"2 10 s ( i i i i i i i ) 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Frequency rad/sec Figure 6.18: Closed-loop process spectral check for Case 3 The closed-loop predicted process power spectrum was obtained by substituting the parameter values in Table 6.12 into (6.8). The predicted process spectrum is plotted in Figure 6.18 (solid line). In order to validate the predicted process spectrum, MATLAB function SPECTRUM.M was used to compute the measurement spectrum using the steady state operational data. For comparison, the measurement spectrum is also plotted against the predicted spectrum in Figure 6.18 (dashed line). Figure 6.18 shows that the averaged measured spectrum is larger than the averaged predicted spectrum (the dashed curve is above the solid curve). This means that there are 83 Chapter 6: Validation of the Tuning Algorithms some errors in the prediction of the signal's average energy. This is confirmed by the results in Table 6.11, where the prediction error for the process 2 a is 63%. Trans, function P(z~') - process model Dfz'1) - Dahlin controller Parameter Td Tr Ta kp , Td TV Value -0.017 106.6 198 378.91 -0.0052 252 189.42 Trans, function Ffz'1) - control filter Mfz'1) -display filter Nfz'1) - noise model Parameter P 9 Cl c2 cr Value 0.2 0.2 -0.47 -0.36 -0.06 0.4585 Table 6.12: Paramer value of each component transfer function used in the spectrum check Closed-loop actuator spectrum check The predicted actuator spectrum was obtained by substituting the parameter values in Table 6.12 into equation (6.9), and the measured actuator spectrum was calculated using MATLAB function SPECTRUM.M based on the steady state data. The predicted actuator spectrum and the measured actuator spectrum are plotted in Figure 6.19 to check applicability of the identified process and disturbance models. The prediction error for the averaged actuator spectrum is about 32% (Table 6.11) and the prediction error mainly distributed over low frequency and high frequency (co < 0.02 and co > 0.06). In this case, the Nyquist frequency is 0.074. The errors of predicting the process and actuator spectrum can be caused by various factors, such as insufficient and/or inaccurate models, changing disturbance dynamics, actuation nonlinearity, unidentified components in the feedback loop, short data sequences used for the model identification, and so on. In this case, since only 20 scans of steady state data were logged and used for computing the "actual" process spectrum and 2 a, the resulting spectrum and 2 a may not 84 Chapter 6: Validation of the Tuning Algorithms represent the characteristics of the process variation very well. Despite all these difficulties, a reasonable qualitative agreement between the prediction and the measured spectra is obtained. 104 Closed-loop model spectrum- solid; spectrum estimate - dashed 10s --101 7 10° l l i i l l l } 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Frequency rad/sec 0.08 Figure 19: Closed-loop actuator spectral check 85 Chapter 7: Conclusions Chapter 7 Conclusions Controller automatic tuning is a highly desirable and useful feature for an industrial control system. The feedback controller tuning tool developed in this thesis is designed to work with industrial Honeywell-Measurex CD control systems and aims at assisting the field personnel in tuning and maintaining the feedback loop of the CD control systems. This would reduce the production loss due to intervention of the controller manual adjustment and improving paper product quality. The main ideas and results presented in this thesis are siunmarized as follows. The tuning of the Dahlin controller and control filter in the feedback loop of the Honeywell-Measurex CD control system requires the process time-response model being identified. The developed algorithm for identification of the process time-response model has been integrated into a control product currently on the market and successfully operated in different paper mills for 9 months. Successful identification was obtained for a wide range of processes, such as basis weight, moisture and caliper. Mill operation of the developed process identification algorithm has shown that it is robust in the presence of strong process noise; it allows reliable and fast identification of a process model using nonuniformly sampled bump test data; it is able to identify the process dynamics by using either the MD-filtered scanner data or the raw data and can identify fractional time delay. In order to tune the Dahlin controller and control filter for optimal disturbance rejection, it is necessary to establish a fairly accurate model of the paper machine process disturbance. To increase data available for identification of the disturbance model, all zone disturbance sequences from the bump test CD residual were cascaded to form a longer disturbance sequence. By doing so, overall accuracy of the disturbance identification was improved due to the significant increase of 86 Chapter 7: Conclusions the data available for the disturbance identification. Since optimality of Dahlin controller for a first order process and first order Integrated Moving Average disturbance has been proved analytically, the used Integrated Moving Average disturbance model allows to verify the developed tuning algorithms by comparing the tuning results to the analytical solution in the literature. It was shown by simulation and mill data tests that the used Recursive Extended Least Squares method had good convergence properties and applicability. For a few cases, when the disturbance data were not sufficiently rich or the order of the disturbance model was improperly chosen (too large), the RELS failed to give a correct estimate. It is possible to correct the "overfitting" by introducing regularization into the basic RELS algorithm, but the estimation failure caused by insufficient data can only be fixed by repeating bump test and extending collection of the bump test data. The tuning of the closed-loop time constant a in the Dahlin controller and control filter factor is based on minimization of the given quadratic performance index. Although the obtained a and J3 are approximate optimal values, the accuracy satisfies the most applications of paper machine CD control. The purpose of introducing the modified tuning performance" index is to cope with uncertainty in the process disturbance dynamics. Theoretically, the modified tuning algorithm guarantees a performance lower bound for a set of disturbances, however, in some cases the modified tuning algorithm overestimates the process and actuator variances and results in too conservative a tuning. The developed feedback controller tuning tool was successfully tested through simulation model, HMDI hardware-in-the-loop paper machine simulator and many sets of mill data. Simulation has shown that the implemented identification algorithms were capable of identifying the process model as well as the disturbance model with satisfactory accuracy. For mill data tests it was shown that the identified disturbance model captured basic characteristics of the process variation and overall the process variation was predicted with a satisfactory degree of accuracy. In 87 Chapter 7: Conclusions order to validate the identified process and disturbance models, the predicted closed-loop process spectrum was compared with the measured spectrum. Fairly good match between the predicted spectrum and the true spectrum indicates that the identified models are applicable to prediction of the real process and actuator variation. In some cases of the spectral check, the mismatch over very low frequency range was observed. This was likely caused by the way in which the disturbance sequences were handled. Putting the zone disturbance sequences into one series would generally add some false low frequency component into the signal, so the information over very low frequency range in the spectrum of the process variation should be disregarded. The observed process variation prediction errors can be attributed to insufficient and/or inaccurate models, changing disturbance dynamics, actuation nonlinearity, short data sequence for the identification of the process and disturbance models and so on. For better tractability of the problem, in the developed algorithm it is assumed that the disturbances in adjacent zones are independent and the actuator CD response is narrow (about 1 or 2 zones wide). In most cases, the developed algorithms give reasonable tuning parameters as long as these assumptions hold. For possible future research, one might reconsider the assumptions made in development of the tuning algorithms and taking the process model uncertainty into account. For basis weight control of heavy grade paper products, the assumptions that the actuator CD response is narrow and the disturbances in adjacent zones are independent are usually not valid, so some kinds of modification of the process and disturbance model structures are needed. 88 ,1 References [1] Chen, S.-C, and Wilhelm, Jr., R. G., "Optimal control of cross-machine direction web profile with constraints on the control effort, Proceedings of the American Control Conference, pp. 1409 - 1415, 1986. [2] Ziegler, J.G. and Nichols, N. B., "Optimal settings for automatic controller", Trans. ASME, 64, 759, 1942. [3] Dahlin, E.B., "Designing and tuning digital controllers", Instruments and Control Systems, Vol. 41, pp 77-83, 1968. [4] Dumont, G. A. "Analysis of the design and sensitivity of the Dahlin regulator". Internal report, Pulp and Paper Research Institute of Canada, 1982. [5] Cluett, W.R., Wang, L. "New Tuning Rules for PID Control", Control System '96, 75 - 80. [6] Astrom, K. J. and T. Hagglund. "Automatic tuning of simple regulators with specifications on phase and amplitude margins". Automatica, Vol. 20, 645 -651, 1984. [7] Gorinevsky, D., Heaven, M., Hagart-Alexander, C, Kean, M and Morgan, S., "New Algorithms for Intelligent Identification of Paper Alignment and Nonlinear Shrinkage", Control Systems '96, Halifax, Nova Scotia. [8] Duncan, S.R. & Corrscadden, K.W., "Minimising the Range of Cross-directional Variations in Basis Weight on a Paper Machine", Proceedings of the IEEE international Conference on Control Applications, 1996. 89 [9] Ljung, L., and Soderstrom, T., Theory and Practice of Recursive Identification, The MIT Press, 1983 [10] Pagariini, F., "A Set-Based Approach for White Noise Modelling", IEEE Trans. Automat. Contr., Vol. 41, No. 10, pp. 1453 - 1465, 1996 [11] Dennis, J.E., and Schnabel, R.B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice Hall, Englewood Cliffs, 1983. [12] Elnaggar, A., Dumont, G.A., and Elshafei, A.L. Adaptive Control with Direct Delay Estimation, Pulp & Paper Canada, T149-154, 1994. [13] Ljung, L. System Identification: Theory for the User. Prentice-Hall, 1987. [14] TAPPI, TIP 1101-01, 1996 [15] Heaven, E. M., et al. "Recent advances in cross-machine profile control," IEEE Control System Magazine, pp. 36 - 46, Oct. 1994. [16] Astrom, K. J. " Introduction to Stochastic Control Theory", Academic Press, Inc., New York, 1970. [17] Ljung, L., Soderstrom, T. and Gustavsson, I. "Counterexamples to General Convergence of a Commonly Used Recursive Identification Method" , IEEE Trans Auto. Cont. Vol AC-20, No. 5, 1975, 643-652. [18] Zhang, M. and Gorinevsky, D., "Design and Development of IntelliTune", Internal report, Measurex Devron Inc., Sept. 1996. 90 [19] Astrom, K. J., Hagglund, T., Hang, C.C. and Ho, W.K., "Automatic Tuning and Adaptation for PID Controllers - A Survey", Control Eng. Practice, Vol. 1, No. 4, pp. 699 - 714, 1993 [20] Zhang, M., "Validation of the Disturbance Identification and Control Loop Variation Prediction", Internal report, Honeywell-Measurex Devron Inc., October, 1997 [21] Zhang, M.,Gorinevsky, D., "Robust tuning of CD controller Feedback Loop", Internal report, Honeywell-Measurex Devron Inc., October 1997 [22] A.P Kjaer, W.P. Heath, and P.E. Wellstead. Identification of cross-directioin behavior in web production: techniques and experience. Control Engineering Practice; 3(1): 21-30, 1995 [23] Astrom, K. J. and Wittenmark, B., Computer Controlled Systems, Prentice-Hall, New Jersey, 1984. [24] Wang, X. G., Dumont, G. A. and Davies, M. S. (1993a), "Estimation in paper machine control", IEEE Control Systems 13(8): 34-43 [25] Astrom, K. J. and Wittenmark, B., Adaptive Control, Second Edition, Addison-Wesley Publishing Company, Inc., 1995. [26] Box,, G. and Jenkins, G., Time Series Analysis, Holden-Day, San Francisco, California, 1976. [27] Oppenheim, A. V., and Schafer, R. W., Digital Signal Processing, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1975. 91 [28] Braatz, R. D., Ogunnaike, B. A. and Featherstone, A. P., "Identification, Estimation, and Control of Sheet and Film Processes", 13th Triennial World Congress, IF AC, San Francisco, USA, 1996. [29] Ghofraniha, J., Davies, M. S., and Dumont, G. A., "CD Response Modelling for Control of a Paper Machine", The 4th IEEE Conference on Control Applications, Albany, N.Y. Sept. 1995. [30] Dumont, G. A., "Application of Advanced Control Methods in the Pulp Paper Industry - A Survey" Automatica 22: 143 - 153, 1986. [31] Heath, W. P. and Wellstead, P. E., "Self-tuning Prediction and Control for Two-dimensional Processes Part 1: Fixed Parameter Algorithms", Int. J. of Control 62: 65 - 107, 1995. 92
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Automated tuning of paper machine cross-directional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Automated tuning of paper machine cross-directional feedback controllers Zhang, Ming 1997
pdf
Page Metadata
Item Metadata
Title | Automated tuning of paper machine cross-directional feedback controllers |
Creator |
Zhang, Ming |
Date Issued | 1997 |
Description | Automatic tuning of the feedback controllers in a commercial paper machine CD control system is considered in this thesis. A feedback controller tuning tool has been designed and developed for tuning the Dahlin controller and control filter in the Honeywell-Measurex CD control system. This tool includes the process identification module, the disturbance identification module and the controller and filter tuning module. The tool is coded in MATLAB and some algorithms have been implemented as "C" modules embedded into a product currently on the market. The tuning strategy used in this tool is based on modern system identification and controller parameter optimization techniques. A global search method combined with the Levenberg- Marquardt refinement is used for identification of the paper machine time-response model (first order plus delay). To obtain a longer disturbance realization sequence and increase data available for identification of the process disturbance, CD residuals, extracted from bump test data, are used for identifying an Integrated Moving Average disturbance model. The identification method is based on the Recursive Extended Least Squares. Based on the identified process and disturbance models, the Dahlin controller and control filter are tuned to minimize a LQG based performance index. Since disturbance dynamics for a paper machine vary with time and with the operating conditions, the possible changes of the disturbance characteristics should be taken into account in the developed tuning algorithms. In the modified tuning scheme, the controller and filter are tuned in such a way that a performance lower bound is guaranteed for all disturbances in a certain set. The feedback controller tuning algorithms were tested extensively using simulation models, Honeywell-Measurex Devron hardware-in-the-loop paper machine simulator and many sets of mill data. A mill trial was carried out for validation of the disturbance identification and control loop variation prediction. Simulation and mill data tests showed that the developed identification algorithms were applicable to real-life mill data and the process variations were predicted with satisfactory accuracy. For better tractability of problem, in the developed algorithms it was assumed that the disturbances in adjacent zones were independent and the actuator response was narrow (about 1 or 2 zones wide). In most cases, the developed algorithms gave reasonable tuning parameters as long as these assumptions held. Validation using mill data also showed that the modified tuning algorithm sometimes overestimated the process and actuator variances and resulted in too conservative a tuning. |
Extent | 6570045 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-29 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065212 |
URI | http://hdl.handle.net/2429/7692 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1998-0076.pdf [ 6.27MB ]
- Metadata
- JSON: 831-1.0065212.json
- JSON-LD: 831-1.0065212-ld.json
- RDF/XML (Pretty): 831-1.0065212-rdf.xml
- RDF/JSON: 831-1.0065212-rdf.json
- Turtle: 831-1.0065212-turtle.txt
- N-Triples: 831-1.0065212-rdf-ntriples.txt
- Original Record: 831-1.0065212-source.json
- Full Text
- 831-1.0065212-fulltext.txt
- Citation
- 831-1.0065212.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065212/manifest