UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Cross directional response modeling, identification and control of dry weight profile on paper machines Ghofraniha, Jahangir 1997

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

Item Metadata


831-ubc_1997-250571.pdf [ 9.17MB ]
JSON: 831-1.0065124.json
JSON-LD: 831-1.0065124-ld.json
RDF/XML (Pretty): 831-1.0065124-rdf.xml
RDF/JSON: 831-1.0065124-rdf.json
Turtle: 831-1.0065124-turtle.txt
N-Triples: 831-1.0065124-rdf-ntriples.txt
Original Record: 831-1.0065124-source.json
Full Text

Full Text

Cross Directional Response Modeling, Identification and Control of Dry Weight Profile on Paper Machines.  by Jahangir  Ghofraniha  B.Sc. Sharif University of Technology, Tehran, 1985. M.A.Sc. University of British Columbia, Vancouver, 1990.  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  in  THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING  We accept this thesis as conforming to the required standard  i  THE UNIVERSITY OF BRITISH COLUMBIA July 1997  © Jahangir Ghofraniha, 1997  In  presenting this  degree  at the  thesis  in  partial  University of  fulfilment  of  of  department  this thesis for or  by  his  or  requirements  British Columbia, I agree that the  freely available for reference and study. I further copying  the  representatives.  an advanced  Library shall make  it  agree that permission for extensive  scholarly purposes may be her  for  It  is  granted  by the  understood  that  head of copying  my or  publication of this thesis for financial gain shall not be allowed without my written permission.  Department  &L€C7£It/lL  of  The University of British Columbia Vancouver, Canada  Date  DE-6  (2/88)  fag  2V,  111?.  ^Glfi/EERlh)G>  Abstract  In this thesis, dry weight response modelling, identification and control of the cross directional (CD) profile of a paper sheet is investigated. To date, the available modelling strategies are either very detailed and require computational fluid dynamics techniques or are based on purely input/output identification schemes. The developed model is based on physical principles and is simple enough to be used for C D control purposes. The model is developed as a dispersive wave at the surface of slurry.  The combination of surface wave theory and drainage theory results in a closed form  model that can simulate different scenarios on the paper machine. Based on the developed model, an identification scheme is developed to identify parameters. The identification scheme makes use of continuous wavelet transform. The identification was applied to simulated data and industrial data. The continuous wavelet transform led to a new C D control algorithm that eliminates the mapping problem which is a draw back of C D control systems due to ill-conditioning and lack of inverse at all times. The control algorithm is demonstrated by application to industrial data.  ii  Table of Contents  Abstract  1  List of Tables  vii  List of Figures  viii  Dedication  xii  Acknowledgments  xiii  Nomenclature  xiv  Introduction  1  1.1  1  Background  1.1.1  Definition of Papermaking Terms  3  1.2  Motivation and Objectives  4  1.3  Contributions  5  1.4  Previous Work  6  1.4.1  C D Response Model Survey  6  1.4.2  Drainage Survey  9  1.4.3  C D Control Survey  1.5 2  ii  11  Thesis Overview  14  Cross Directional Response Model  15  2.1  Introduction  15  2.2  Dispersive Wave Model  15  2.2.1  Equation of Motion  17  2.2.2  Boundary Conditions  20  2.2.3  Solution for exponentially damped waves  22  iii  3  2.3  Drainage Equation and Mat Formation  26  2.4  Pressure Profiles on Drainage Elements  29  2.5  Mat Compaction  34  2.6  Superposition and Response Model at the Edge  35  2.7  Simulation Results  37  2.8  Summary  39  Testing and Experimental Methods  46  3.1  46  3.1.1  Purpose of the Experiment  46  3.1.2  Scope of the Experiment  47  3.1.3  Data Collection Scheme  47  3.1.4  Background  47  3.2  4  Introduction  Multirate Filtering  48  3.2.1  Aliasing  49  3.2.2  Filter Design  50  3.3  C D Profile Calculation  52  3.4  Calibration of Basis Weight Sensor  53  3.4.1  Test Specifications  53  3.4.2  Analysis Scheme  54  3.5  Results  55  3.6  Summary  58  System Identification  65  4.1  Introduction  65  4.2  System Identification Approach  65  4.2.1  The Problem  65  4.2.2  Analogy  67 iv  4.3  68  4.3.1  Introduction  68  4.3.2  Definitions and Terminology  70  4.3.3  Continuous Wavelet Transform (CWT)  74  4.3.4  Correlation Analysis . .  76  4.3.5  Proposed System Identification Approach  77  4.3.6  Alignment  80  4.4  Parameter Estimation  82  4.4.1  Estimation with Simulated Data  82  4.4.2  Estimation with Industrial Data  84  4.4.3  Model Validation  85  4.4.4  Scale Bound Considerations for C W T  88  4.4.5  Choice of Mother Wavelet Parameters  90  4.5 5  Wavelet Theory and Correlation Analysis  Summary  90  Cross Directional Profile Control  96  5.1  Introduction  96  5.2  C D Profile Control with Wavelets  97  5.3  Overview of the Control Scheme  99  5.4  Control Algorithm  100  5.4.1  Slice Lip Deflection Model  101  5.4.2  C D Profile Decomposition with C W T  103  5.4.3  Thresholding  104  5.4.4  Reconstruction of C D Profile for Control  105  5.4.5  Change of Scale (Profile Compression)  106  5.4.6  Dahlin Control  107  5.5  Simulation of C D Control Algorithm with Industrial Data  108  5.6  Summary  110 v  6  Conclusions and Future Work  113  6.1  Conclusions  113  6.2  Future Work  115  A  116 A.l  Solution of Laplace's Equation by Separation of Variables  A.2 Signal Processing Scheme of the Measurex Scanner  116 118  A.2.1  A / D Conversion in V F C  118  A.2.2  Data Compression  119  A.2.3  Basis Weight Calibration .  120  A.2.4  Moisture Calibration  122  A.3 Definition of Norm  122  A.4 Calculation of Admissibility Function  124  Bibliography  127  vi  List of Tables  2.1 Machine parameters used for simulation of 5 scenarios mentioned in the list  38  3.2 Results of manual measurements of basis weight across the sample sheet  55  3.3 Comparison of the Measurex profile with different profiles in the squared error sense. . . . 57 3.4 Comparison of hand measurement of basis weight with different profiles in the squared error sense  58  4.5 Mother wavelet parameters used in Figure (4.36)  82  4.6 Estimated parameters for the example in Figure (4.36)  83  4.7 Mother wavelet parameters used in Figure (4.37)  84  4.8 Estimated parameters for the example in Figure (4.37)  84  4.9 Mother wavelet parameters used in Figures (4.38 and 4.41)  87  4.10 Estimated parameters for the example in Figure (4.38)  87  4.11 Estimated parameters for the example in Figure (4.41)  87  A. 12 Possible solutions of Laplace's equation depending on the value of separation constant from Reference  117  vii  List of Figures  1.1 A simplified schematic of a paper machine  2  2.2 Wave crests generated by a ship on a calm sea  16  2.3 Schematic of wave generated by slice deflection on a porous bed (mat)  18  2.4 Surface S surrounding the fluid  19  2.5 Total drainage resistance vs. basis weight ( from reference )  27  2.6 Pressure profile over a table roll for a wire speed of 1500 f.p.m (from reference )  31  2.7 Effect of wire speed on pressure profile over a table roll (from reference )  31  2.8 Schematic of some foils used in modern Fourdrinier machines (from and  32  2.9 Definition of the parameters used in Taylor's formulation ( from reference )  33  2.10 Typical pressure profiles over a foil and a foil with a vacuum pump  33  2.11 Headbox induced edge waves due to misalignment of cheeking elements  35  2.12 Proper use of edge bleeds reduces the edge wave effect (from reference )  36  2.13 C D basis weight response to an adjustment of a slice lip actuator for a heavy grade paper (All C D and M D position units are in samples. Dry line is an average location of the region)  41  2.14 C D basis weight response to an adjustment of a slice lip actuator for a light grade paper (All C D and M D position units are in samples. Dry line is an average location of the region)  42  2.15 Basis weight response at the edge for different grades of paper (All C D and M D position units are in samples)  43 viii  2.16 Effect of changing the drainage elements on a Fourdrinier machine and the corresponding shift of the dry line location (All C D and M D position units are in samples)  44  2.17 Superposition of surface profiles in two different cases. Figures (a), (b) and (c) refer to two adjacent bumps with the same magnitude and the cases (d), (e) and (f) with different amplitudes (All C D and M D position units are in samples)  45  3.18 A system diagram for Measurex scanner  48  3.19 Multirate system building blocks, M-fold decimator L-fold expander  49  3.20 Sampling rate reduction without prefiltering results in aliasing  50  3.21 Minimum requirement for a typical decimation filter when subsampling by M  51  3.22 A schematic of the sample used for calibration. The circles show the position of puncher for manual measurement of the basis weight  54  3.23 Calibration data collected by the data acquisition system before processing  57  3.24 C D profile calculated based on a Rectangular window averaging scheme  60  3.25 C D profile calculated based on a Hamming window averaging scheme  61  3.26 C D profile calculated based on a Hanning window averaging scheme  62  3.27 Comparison of the Measurex profile with profiles obtained from different averaging methods  63  3.28 Comparison of different basis weight profiles with that of the hand measurements of the sample sheet  64  4.29 Distortion in the Fourier analysis for frequencies in adjacent intervals and improvement by Short Time Fourier Transform (STFT). Improvement of Wavelet Transform over STFT by detecting the discontinuity at sample 128 (t=l sec)  69  4.30 Multiresolution wavelet decomposition scheme  74  4.31 Wavelet packet algorithm  75 ix  4.32 A schematic of the identification procedure with Continuous Wavelet Transform  81  4.33 Determination of misalignment using C W T  82  4.34 Results of a bump test at three different zones across the lip (Courtesy of Measurex Devron Inc) 4.35 Bump test on a liner board machine from reference  86 showing a typical M-shape  response  86  4.36 Identification of a simulated response on a typical light grade with Continuous Wavelet Transform  91  4.37 Identification of a simulated response on a liner board machine with Continuous Wavelet Transform  92  4.38 Identification of a bump test response on a 81 lb/3300 sqft grade with Continuous Wavelet Transform  93  4.39 Normalized autocorrelation of estimation error with 95% confidence interval for the example shown in Figure (4.36)  94  4.40 Normalized autocorrelation of the estimation error and 95% confidence interval for the example shown in Figure (4.38)  94  4.41 Identification of a bump test response on a liner board machine (Industrial data taken from ) with Continuous Wavelet Transform  95  5.42 Schematic of C D / M D control of dry weight profile with Continuous Wavelet Transform  100  5.43 A schematic of the slice lip as a beam with loading  102  5.44 Schematic of the Hard and Soft thresholding algorithms  105  5.45 Different stages of C D control algorithm shown through plots (a) to (h)  Ill  5.46 Residuals after the implementation of C D control and control action  112  A.47 Schematic diagram of a voltage to frequency converter  119  A.48 Schematic of the averaging process viewed as a mapping of many-to-one process.  . . .  120  A.49 Straight averaging viewed as a decimation block  120  A.50 A typical transmission curve used in Measurex calibration scheme  121  A.51 A typical calibration curve used in Measurex scanners  123  xi  Dedication  I dedicate this thesis to Bahar, Zahra and Maryam and my parents Taghi and Mina Ghofraniha.  xii  Acknowledgments  A Ph.D. thesis might not qualify as a team effort but it will definitely require the contribution of many individuals. I would like to take this opportunity to thank those who played an important role in this work. First, I would like to thank my supervisors, Dr. Michael Davies and Dr. Guy Dumont for their support and guidance and friendship during all those years at U B C . They taught me how to think critically and independently. The efforts of my committee members Drs. Farrokh Sassani, Ezra Kowk, Chris Ma and Avron Soudack are greatly acknowledged. I would like to thank my past and present fellow students at the Pulp and Paper Centre Control Group, Ashraf Elneggar, Ye Fu, Abdoul-Latif Elshafei, Kristinn Kristinsson, Scott Morgan, Navid Mohsenian, Lahoucine Ettaleb, Roger Shirt and Zoran Nesic for their friendship and heated discussions that eased the graduate school experience. The help of the Pulp and Paper Centre staff, especially Rita Penco for all her help with library research is greatly appreciated. I would like to pay special tribute to my family for all they have done for me. M y heart goes to Bahar, Zahra and Maryam. This thesis has affected them more than anyone else. Their continuous love and support carried me through the hardship of graduate school. I am indebted to my parents, brothers, sister and in-laws for their support and love and believing in me. I would like to thank Sandra Davies for her friendship and specially for proof reading the thesis. Special thanks go to all my Iranian friends who made the student life for me and my family an enjoyable experience. The financial support of British Columbia Science Council, Paprican and the N C E Mechanical WoodPulp Network is greatly appreciated.  xiu  Nomenclature  Headbox consistency (%)  c  s  Filtrate consistency (%) Mat consistency (%) Ah  Change in free surface of slurry (m)  v  f  Filtrate velocity (m/s)  Lit)  Mat thickness (m)  Bw  Dry weight (gr/m )  Ap  Pressure drop across mat (Kpa)  V-  Viscosity of the slurry (kg/m.sec)  K  Permeability constant (1/m)  rt  First-pass retention (%)  rt  Initial retention (%)  p  0  2  Initial filtrate consistency (%) Constant in retention equation Ps  Density of slurry (kg/m )  A, B, D  Constants  k  Spatial frequency (cycle/m)  CO  Angular frequency (Hz)  3  Damping factor in wave equation V  Wave amplitude (m) Velocity potential  Pb  Pressure function at the mat (Kpa) Amplitude in wave equation (m)  PS  Pressure in the fluid (Kpa)  u  Velocity field (m/s)  N  Complex wave number  q  Velocity term in Bernoulli's equation (m/s)  CD  Cross direction  MD  Machine direction  xiv  D  Amplitude in C D model (m)  ARMAX  Auto Regressive Moving Average time series with Exogeneous term  MIMO  Multi-input Multi-output system  V xf  Curl operator, curl of function /  V/  Divergence operator, divergence of function /  D/Dt  Fluid mechanics notation for total derivative  V /  Laplacian o f /  A  m  Cross section area in Darcy's Law (m )  L  m  Filter thickness in Darcy's Law (m)  0  2  2  R  Flow resistance in Taylor's table roll model  rb t  Table roll radius (m)  U  Wire speed (m/min)  K'  Constant in table roll equation  M, Ni  Constants in mat compaction imperical formula  kdj  Initial damping in dispersive wave equation  c  Wave phase velocity (m/s)  SFR  Specific Alteration resistance (m/kg in SI or cm/gr)  R,  Total mat resistance (1/m)  R  Wire resistance (1/m)  Kr85  Krypton 85  D(n,k)  Moisture content in Lindeborg's model (%)  Mr  Reference moisture in M D (%)  p"  C D moisture profile (%)  B  Nonlinear interaction term between M D / C D in Lindeborg's model  CWT  Continuous Wavelet Transform  DWT  Discrete Wavelet Transform  STFT  Short Time Fourier Transform  t/j j,  Wavelet function in continuous domain, a = scale & b = shift  w  a  Admissibilty function W^f  Wavelet coefficients in continuous domain  ip „  Wavelet function in discrete domain  ao, bo  Mother wavelet scale and shift parameters  A', B'  Frame bounds in wavelet definition of Frames  MW  Mother wavelet function  m>  xv  WB  Wavelet bases  a  iP> Pip  Wavelet coefficients in discrete case  y(x), u(x)  Dry weight and actuator profiles (gr/m & microns)  G  Interaction matrix  Pi  Radius of curvature (m)  M(x)  Bending moment in slice lip (Nm)  EI  Flexural rigidity (N.m )  2  2  Spring constant in flexible actuator rods (N/m) GDC  Dahlin controller pulse transfer function  A  Desired closed loop time constant in Dahlin controller (sees) System gain in Dahlin controller formula (gsm/microns)  7  Constant representing the area in drainage equation (m ) 2  XVI  Chapter I:  Introduction  Chapter 1 Introduction  1.1 Background The art of papermaking dates back at least to 105 A . D . when the Chinese inventor Ts'ai Lun described a technique to drain a mixture of bamboo fibres and water through a woven screen [43]. Handmade papermaking was the only available method until 1804 when the first commercial paper machine was designed by the Fourdrinier brothers [81].  Although paper machines have evolved  since, the principles have not changed dramatically. The main constituent of a sheet of paper is pulp. Wood structure, in simple terms, consists of cellulose fibres, lignin and resin materials. The purpose of the pulping process is to extract the fibre from the wood structure by the lignin material that binds the fibres together. Two main pulping methods, chemical and mechanical, are employed to produce pulp each with different application and end results.  Chemical pulping produces fibres that are  long, flexible and strong and normally used in products for which strength is the main characteristic, such as linerboard and sack paper. One of the main consumers of unbleached chemical pulp is the packaging industry, while bleached chemical pulp is used widely in high quality and high end printing. Mechanical pulp has a higher yield compared to chemical pulp (90-95% compared to 40-50%), but does not have the strength and color stability of chemical pulp [82]. The modern papermaking process starts with dissolving the pulp in a large storage tank known as the machine chest to a consistency of 3-4 % for a fine paper. The slurry is then diluted with a low consistency mixture of water and pulp and passes through stages of cleaning and screening before entering the headbox. Although the Fourdrinier headbox is in widespread use other designs are in increasing use, and avoid the one sided initial drainage [92]. A simplified schematic of a Fourdrinier type paper machine is shown in Figure (1.1). The mixture of fibre and water is delivered to the headbox. The headbox is a high pressure container that delivers the mixture onto a mesh conveyor (wire). The headbox is responsible for providing a uniform flow onto the wire by reducing fluctuations in the consistency of stock (water/pulp mixture) and irregularities across the width of the machine. The headbox also 1  Chapter I:  Introduction  Figure 1.1: A simplified schematic of a paper machine. prevents the fibres from binding together to form floes. The primary dewatering takes place on the wire when water drains from the slurry as it passes over the dewatering elements, and as a result, the fibre sheet forms towards the end of Fourdrinier section. The second stage of the dewatering is carried out in the press section where excess water is squeezed out of the sheet. The moisture of the sheet reaches its final target (roughly 4%) in the dryer section of the paper machine. The thickness of the sheet becomes more uniform when the sheet passes through a series of rolls called the calender stack (in some installations the calendering is an off-machine operation). 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 adjustment of properties across the sheet are done through a series of actuators mounted across the machine. The weight actuators are mainly the slice screws at the headbox. For consistency profiling headboxes, the weight actuators are a series of valves (linear or globe valves) that control the local consistency across the headbox. The control of moisture profile across the sheet is carried out by steam boxes, normally mounted close to the wet end of the paper machine (close to press section), and steam showers installed at the dry section of the paper machine. The function of the steam boxes is mainly to increase the local temperature of the sheet such that water can be more easily squeezed out of the sheet at the press section. Steam showers remove moisture by locally  Chapter I:  Introduction  exposing the sheet to dry hot steam. Water sprays are used to increase the moisture at different C D locations. The combination of steam showers and steam boxes with water sprays provides a wide control range for the moisture across the sheet. Most actuators that control the properties of the sheet in the direction of sheet movement are before the headbox [82] except for drying cylinders (can dryers) and Infrared dryers which are mounted at the dryer section of the paper machine.  1.1.1 Definition of Papermaking Terms Although some terms have already been introduced, a more organized definition of some of the papermaking terms used in this thesis are presented in this section. The sheet properties that are most commonly used in the paper industry are basis weight, moisture and caliper. Basis weight is defined as the total sheet weight per unit area, and dry weight (moisture removed) as the weight of the dry solids per unit area of the sheet. Moisture is the percentage by weight of the water content to the total sheet, and caliper is a measure of mat thickness under specific surface squeeze pressure. During manufacturing, process fluctuations in the paper machines' approach system tend to cause disturbances that effect the entire sheet. Lack of uniformity on the machine itself tends to cause irregularities across the sheet as it is produced. It is therefore customary to decompose property variations such as weight and moisture of the sheet into components in cross direction (CD), machine direction (MD) and residual effects [11].  In the same manner, actuators that affect the properties  of the sheet in C D and M D are referred to as C D and M D actuators. The number of Databoxes represents the number of measurement points at the scanner resolution. A C D profile with 400 data boxes is equivalent to a profile with 400 samples. This resolution is normally reduced to the lower resolution of actuators for C D control calculations. C D weight control on a Fourdrinier headbox is normally carried out by means of a series of long screws mounted on a beam of a small cross section at the opening of the headbox.  This opening is at the bottom of the headbox facing the  wire (mesh conveyor) and is the orifice that delivers the stock onto the wire at a desired angle and velocity under headbox pressure. The beam is referred to as the slice lip and serves as the upper boundary of the orifice. The long rods that are attached to the slice lip are called spindle rods and are either motorized or actuated by thermal expansion. The lower part of the opening of the headbox  3  Chapter 1:  Introduction  is called the apron. After leaving the headbox, the stock passes over the drainage elements and fibre is deposited on the mesh conveyor, and some water carrying a small amount of fibre is drained. This process continues as the stock moves with the wire towards the press section of the paper machine. After the free moisture drains, the reflectance of the sheet surface changes over a narrow region, and a distinct line can be seen across the wire. This line is called the dry line and is an indication of where the sheet changes from a flowing suspension to a fibre matrix. Machine operators normally use the location of the dry line as a visual indicator of the drainage quality on the Fourdrinier. Recently, techniques for controlling C D basis weight by adjusting stock consistency across the slice have been introduced. Such "dilution headboxes" allow additional water to be injected differently across the sheet, so reducing fibre content for fixed flow and a corresponding local reduction in basis weight [92]. The measuring device that is mounted on a frame and measures the weight, moisture and caliper in a zigzag manner is referred to as the scanner [82].  1.2 Motivation and Objectives The reduction of variability is one of the keys to successful paper making. In Canada, the quality of fibres and pulp gives an edge to Canadian pulp and paper industry in the world market [7]. Having a good quality raw material does not guarantee a good finished product unless quality is maintained in the production stage. In order to meet the increasing requirements for quality in the paper industry, reduction of variability has become the focus of paper makers. As the Canadian industry moves to higher value-added grades, the importance of close control of the manufacturing process increases.  The nature of variations in paper has been investigated in the past [22].  The  variations can be assigned to different sources depending on their frequency. The paper machine control objective is to reduce those variations that are within the bandwidth of the actuators. With a proper control strategy, not only are the stringent quality requirements met, but better runnability and smooth operation of the paper machine are also achieved. Inadequate control performance may result from poor tuning, but is often caused by a lack of proper models. In the C D control of weight and moisture, a good response model that formulates the interaction between adjacent actuators will be a key element in the reduction of C D variations.  Chapter I:  Introduction  In the past, the mathematical modelling of the wet end of Fourdrinier paper machines emphasized the filtration dynamics.  In recent years, some researchers have also considered fluid dynamics  of the Fourdrinier wire and headbox [98].  These models are very complicated and require finite  element numerical techniques to solve the partial differential equations with the boundary conditions. Unfortunately, these modelling techniques are often of little use in improving C D control, and control engineers have resorted to simple models of the relationship between C D actuators movements and measured response at the scanner. In this work, the objective is to develop a model that relies on physical principles at the wet end of the paper machine and is simple enough to provide a parametric description suitable for C D control. If we show the complexity of modelling techniques on an axis with current physical models at the one extreme and the models used in the current C D control schemes at the other extreme, the current work will be somewhere in the middle of the complexity axis. The second objective is to use experiments on an operating paper machine in order to identify the parameters of the developed model. The final goal is to use the model in a C D control strategy to demonstrate the usefulness of this approach.  1.3 Contributions  The main contribution of this thesis is the development of a model based on the fluid mechanics of the wet end of the paper machine, combining surface wave theory to the surface of the stock with drainage theory. Previously, drainage theory has been treated as a M D modelling technique. The approach presented in this thesis extends the drainage analysis over the width of the sheet, resulting in a two-dimensional analysis of sheet formation. The resulting model is simple enough to be used for C D control purposes. Currently, in the C D control context, the C D response requires many samples to reach a reasonable accuracy. The developed parametric model is a closed form solution of the Laplace equation that represents the C D response of an actuator and requires four parameters, each representing a physical variable. In the experimental part of the thesis, the signal processing scheme applied to the raw scanner data makes use of Multirate Filtering Technique. The mapping of the profile from high resolution to 5  Chapter 1:  Introduction  actuator resolution has been interpreted as a change in sampling rate. This approach to the mapping problem is new. The technique that is used to identify and decompose the variations in C D profile is a Continuous Wavelet Transform ( C W T ), a method that has been borrowed from the field of communication and radar signal detection and processing. Identification of model parameters with C W T is one of the new applications of wavelet theory in the field of system identification and control. By making use of C W T , a novel C D control algorithm is proposed. This approach to C D control eliminates the mapping problem from scanner resolution to actuator resolution and its associated issues such as ill-conditioning of the interaction matrix. The use of this technique in the context of C D response identification and control is unique to this thesis to the best of author's knowledge.  1.4 Previous Work The subject of this thesis covers several areas in engineering. Fluid mechanics, drainage theory and surface wave theory are prominent in the modelling section of the thesis, and multirate filtering and wavelet methods are used in the system identification and control chapters. A n overview of all these fundamental areas is not necessary in the context of a thesis. The literature survey has therefore been concentrated on the modelling component including the three major areas of C D response modelling, drainage theory and C D control of weight profile. A brief introduction to other topics is included in the appropriate chapters.  1.4.1 CD Response Model Survey There is an abundance of research related to cross directional control of sheet profiles in the literature, but specific results on paper machine C D response modelling is very limited. The C D response model refers to the shape of the response to a single movement of a slice lip actuator, normally a screw at the slice lip.  Treatment of the wet end modelling of the paper machine  has traditionally been a part of the C D control problem, or in the drainage modelling of the wet end.  C D response modelling schemes can be divided into three main categories.  The first group  deals with the physical aspects of the process and starts from the underlying principles. 6  These  Chapter I:  Introduction  models include complete hydrodynamic behaviour and leads to nonlinear partial differential equations, normally solvable by numerical schemes. This approach provides a good understanding of the process and is effective when investigating different design scenarios. The second method is an empirical input/output approach and using an estimated matrix to represent the interaction in C D and normally a linear model with dead-time to represent the M D variations.  These methods are suitable for  implementation of the C D profile control. The third method which could also be considered as a special case of the second, treats the C D variations as a superposition of polynomials with the M D dynamics assumed as in the second case. Westermeyer [98] carried out a detailed analysis. He addressed the wet end modelling problem as a separate issue and investigated the fluid mechanics of the Fourdrinier wire. Westermeyer used a control volume approach to set up the model and analyzed the flow characteristics using a threedimensional flow field with boundary conditions. His approach treated the stock above the Fourdrinier wire and the fibre mat as a single system. The partial differential equations describing the flow are solved for the velocity field using a finite element scheme. The method calculates cross flow velocities resulting in weight variations in cross direction. The approach has been used to look at the flow characteristics of the jet at slice lip with different nozzle geometries. An important conclusion of this work is the influence of rush/drag ratio (difference in jet and machine speed) on the amplification of C D variations. This phenomenon was first pointed out by Cuffey [20] and later investigated in more details by Parker et <a/[71]. Bergh and MacGregor's technique [6] is a linear model approach with a diagonal interaction matrix. They introduced the idea of spatial time series in the sheet forming process. In particular, they used a spatial time series (a series as a function of both time and space) to represent the disturbance model. The disturbance model is formulated as a multivariate low order A R I M A (Auto Regressive Integrating Moving Average) model resulting in a sparse symmetric matrix. The time dependence or M D variations is represented by the M A part of the model. This method can also describe the interaction of C D and M D variations. The response model is treated as a steady state C D gain matrix with low order M D model (order < 2) to account for actuator dynamics. Because of the scanning measurements, the analysis is treated as an identification problem with missing data.  7  Chapter 1:  Introduction  The extrapolation of the results to a real paper machine, where there are up to 400 data points at high resolution measurement and up to 120 points in the low resolution case, could be challenging as a result of the dimensionality of the system. Eriksson et al [29] and Karlsson et al [47] at STFI (Swedish Pulp and Paper Institute) have mentioned the use of surface waves on the Fourdrinier as a method of C D response modelling without disclosing the algorithm. Karlsson has pointed out that the implementation of such models in C D profile control at a specific mill in Sweden has resulted in reduction of average moisture by 0.5% and on heavy grade paper there has been an improvement of up to 60% in quality. Wellstead et al [97] model the wet end of the paper machine as a two-dimensional time series. This method treats the process variables at a particular point as a function of two Cartesian variables. The resulting system is a function of discrete indices, m and n, each variable representing cross direction/space and machine direction/time variations. Implementation of the two-dimensional theory to sheet forming process assumes an infinite domain which needs to be restricted for a finite width paper and modified near the edges of the sheet. The domain localization and a recursive parameter estimation technique have been reported [41]. The choice of system order (number of parameters to represent the whole sheet) and the trade off between accuracy and increase in dimensionality are difficult tasks in most two dimensional modelling schemes. The C D response model normally deals with the steady state response of the actuator while with a successful implementation of such models (2—D time series) the dynamic variations could also be controlled. A n important advantage of the two dimensional modelling scheme is the ability to deal with the C D / M D interaction in the same context. Chen [14] has used bump tests on a paper machine to estimate the C D response model parameters. Changes in the dry weight are correlated with changes in the slice lip position and a least-squares algorithm is used for parameter estimation. The estimated parameters are then used to set up the band diagonal elements in the interaction matrix. Duncan [28] approached the C D response identification using a closed loop estimation scheme. In this method a perturbation signal was added to the set points of the actuators to overcome difficulties associated with the identification of parameters in the closed loop. Due to inherent delay between the time the control is applied and the measurement, the open loop response of actuators was extracted  Chapter I:  Introduction  from the closed loop estimates. Acknowledging the dimensionality of the problem if the complete interaction matrix were to be estimated, the perturbations were only added to few actuators spread across the machine and few columns of the interaction matrix were identified. This methodology was implemented on a plastic film extrusion line. The use of orthogonal polynomials to represent the C D profile has been popular in steel industry and has been used to control C D variations of a steel sheet [37]. Application to the paper industry has been pointed out by Keyes [51] in using terms in a Fourier Series to represent variations of C D moisture profile. The term "Modal Control" has been used to describe this technique. A more recent application of a similar approach to C D profile representation and control can be found in [54] where Gram orthogonal polynomials are used. In this approach the interaction matrix is parameterized in terms of the Gram polynomials (the basis functions to represent the profile and the slice lip shape) and the C D control problem reduces to adjusting the coefficients of these polynomials to reach the desired performance. Rigopoulos et al [76] used a Karhunen-Loeve expansion to identify the significant modes of paper machine disturbance profile. The modelling of the eigen-vectors for the expansion was done with an auto-regressive (AR) process. The potential of such models for use in a model predictive control strategy was considered.  1.4.2  Drainage Survey Mathematical modelling of drainage phenomenon on the Fourdrinier machines was developed  in the decade following 1960. There has been extensive research on the drainage modelling in the context of flow through porous media andfiltrationtheory, and an exhaustive overview of the literature is not attempted here. Selected samples of the more relevant articles in this field are reviewed. Early work in the area of paper machine drainage theory is described in the article by Burkhard and Wrist [12]. The research investigates the correlation of the machine speed with the shape of the pressure distribution over a table roll. The results of the experiment on a paper machine with a speed range of 500 — 2500 fpm are compared with predictions by Taylor [87]. The relation of instability of the stock with jet landing location on the table roll and turbulence in the jet is investigated. The amplitude of  Chapter 1:  Introduction  the kick-up (instability in the stock) over a table roll was correlated with the speed of the machine. Meyer's work on the drainage theory at the wet end of a paper machine is one of the cornerstones [62]. The nature of the mat formation on a Fourdrinier machine with variable thickness has prompted Meyer to use a multi-body drag model instead of the capillary model of Kozeny-Carman [16]. These two methods differ in the way the constant in Darcy's Law is calculated. In Meyer's approach the variability of the mat thickness with time is represented with a variable resistance term in Darcy's Law, resulting in a variable consistency in both the mat and the drained white water. A careful mass balance on the wire side results in the calculation of basis weight in terms of mat thickness and consistency. The derived model is used to investigate two cases, namely constant rate filtration and constant pressure filtration. Meyer has pointed out limitations in the resulting model, for example neglecting the initial acceleration in the direction of flow (perpendicular to sheet surface) and flow resistance of the wire. Ingmanson [44] looked at the effect of mat compressibility on the variations of the velocity through the mat. The velocity is the greatest at the wire side and the lowest at the top of the slurry. The inclusion of this velocity variation in the hydraulic flow resistance equation will result in a nonlinear partial differential equation. A n important conclusion arrived at by solving the equation is that increasing slurry consistency gives a higher sheet porosity distribution. In the constant rate filtration case the model is in good agreement with the experiment. Nelson [65] develops the approach in three aspects. Flow resistance has been included when the solid phase is in motion. In addition, the most significant inertial forces have been included and, in a special case, it has been shown that the problem can be formulated with a single partial differential equation. The extension of the model is used to describe retention of small particles in a porous medium. Initial retention has been calculated by different numerical methods. Han et al [39] tried to simplify the complicated formulation of the drainage on the Fourdrinier for applications in paper machines. Fof low consistency stock (<0.5%) the effect of differing velocity between water and fibre in drainage may be neglected. Solutions for rapid drainage and slow drainage with complete fibre retention are presented. The importance of the wire resistance in forming thin sheets, and the importance of viscous forces in the formation of thick mats are described through examples. Drainage resistance measured with an off-line apparatus is compared with that measured on suction formers directly and found to be in disagreement. The  10  Chapter I:  Introduction  discrepancy has been related to the mat elasticity and compression on the machine. A n excellent survey of the drainage theory up to 1978 is by Branion et al [10]. A detailed critical review of the models by Meyer, Ingmanson and others is included in this monograph. Pires et al describe a computational model for drainage on the Fourdrinier in [73]. The model is based on a simplified version of Darcy's equation and Taylor's model on table rolls. The model is used to predict the location of the dry line.  The specific resistance of the mat (slope of the  increase in mat resistance along the paper machine) is calculated in the laboratory. The apparatus measures the dynamic specific resistance, which differs from that obtained under constant pressure and constant rate filtration conditions. The results indicate an agreement with model prediction to within 10%. Fielding et al [30] describe the methodology used to establish drainage curves on forming tables using a back scatter gauge. Although this article does not contain any modelling scheme, the experimental determination of the drainage curve is critical in the evaluation of performance under different operating conditions. For example, the impact of an increase in speed of the machine or changes in consistency of the headbox on wire drainage can be estimated.  Sayegh et al have  reconsidered the compressibility of the mat during drainage recently. They have found that specific filtration resistance of the mat is a nonlinear function of the applied pressure and of the time variation in pressure. The developed model to account for these effects is based on a Maxwell element which assumes the mat undergoes an elastic/plastic deformation. This approach results in three components in filtration resistance, namely, a zero-pressure term, a pressure-history term and an instantaneous pressure term.  1.4.3  CD Control Survey  The C D control problem is closely related to the resolution of the measuring instruments. If the C D profile measurements are widely spaced, or are heavily filtered, the controller cannot correct for short wavelength variations. In modern paper machines the measurement system normally provides the high resolution profile for control purposes. C D control performance is normally evaluated using the standard deviation of the C D profile. It is important that this standard deviation is calculated using the high resolution C D profile (standard deviation based on the low resolution profile is always 11  Chapter I:  Introduction  smaller than that based on high resolution profile due to filtering). Basis weight profile control is widely considered in the literature, using separate approaches to the two topics of C D control and M D control. More recent work has tried to combine C D and M D control techniques [38]. There have been surveys on the topic of C D profile control in recent years [52] and [93], with a more recent survey to be found in [54]. Although the focus of this survey is the C D control profile, the work by Astrom is briefly mentioned due to its historical value. Astrom's work [1] on paper machine M D control is one the early implementations of a computer control system in an industrial application. Dynamics were represented by an A R M A X model with the disturbance formulated as colored noise. Astrom used a maximum likelihood method to identify the disturbance and the plant dynamics and then implemented a minimum variance control law [1]. Early work on C D profile control was carried out by Beecher et al [4].  Beecher and Bareiss implemented the first successful application of an  analog controller to the basis weight profile. An interaction matrix was used, the elements of which were estimated by minimizing the sum of estimated errors. To simplify the control law, the interaction matrix had to be inverted. The authors claimed 10% to 50% improvement over manual control. A majority of other researchers have taken the same route, controlling the C D profile using different cost functions.  A Linear Quadratic Gaussian formulation with even more complex functions has  been pointed out by Boyle [9]. The major drawback in this technique is the numerical difficulty of calculating the inverse of a large matrix, and of detecting situations in which it is singular [9]. Boyle showed that the bump response that establishes the columns of the interaction matrix can usually be treated as uniform across the sheet, and so reduced the interaction matrix to a symmetric matrix. This simplification makes the inversion process simpler, but still the singularity issue remains unresolved.  Chen used the interaction matrix with a quadratic index [15] to derive an optimal control law. A penalty function to limit lip bending was included. Chen et al have also applied adaptive control laws based on Quadratic Programming to control the B W profile [83].  Laughlin and Morari [57] used a robust control method for the B W control problem. The model used is a first order to account for actuator dynamics, and a band-diagonal matrix with a delay. The C D control problem is formulated as a multivariable (MIMO) control system based on Internal 12  Chapter 1:  Introduction  Model Control (IMC) strategy. By calculating the bounds on the eigenvalues of the model, a robust controller is achieved.  Graser [35] analyzed the C D profile using a Fourier analysis approach. The controllability of variations in C D profile has been correlated with the width of the actuator response, and then the control problem is formulated as an optimization problem with penalty on the slice bending. The optimization algorithm confirms the well known fact that a narrow response results in better control of C D frequencies.  The effect of increasing or decreasing the penalty function on the smoothing  of the actuator positions has been demonstrated. It is concluded that Fourier analysis alone can be misleading when analyzing a C D profile, and a combination of the Fourier analysis and slice lip bending, together with the width of the response and headbox characteristics should be considered at the design level.  The idea of coordination of multiple actuators (MD/CD and CD/CD) has received attention in recent years. Hall [38] addressed the application of multiple C D actuators to control a single profile and showed that economic benefits may be obtained when both actuators are used in coordination as opposed to being used separately. He also investigated the improvements achieved by coordinating M D actuators with the C D actuators. He showed through simulations that this coordination allows for a better control of high frequency M D variations. Kastanakis et al [49] looked at C D / M D interaction problem on paper machines and proposed a procedure to identify this interaction. The decoupling parameters are identified using a frequency analysis approach which could be extended to investigate the achievable control and actuator performances. It is concluded that although standard PI and Smith controller have limited flexibility, a dedicated controller based on an on-line identification scheme is not always justifiable. Kwok et al [55] developed a static/dynamic model of C D / M D interaction based on consistency variations. The model describes variations of the C D dry weight profile in terms of changes in stock valve. The interaction model is identified with industrial data. A generalized predictive control strategy was shown to reduce variability when applied to the C D / M D interaction model in a simulated environment.  13  Chapter 1:  Introduction  1.5 Thesis Overview The development of the model using combined surface wave theory and filtration theory are discussed in chapter 2. Data collection and experimental work on a paper machine and multirate filtering of measurements are explained in chapter 3. Chapter 4 deals with the system identification of the parameters in the model based on wavelet transform methodology using simulated and industrial data. In chapter 5 the C D control problem based on the developed model and continuous wavelet transform is discussed, together with interpretation of the controllability of C D profile in terms of wavelet resolution. Finally, in chapter 6, a summary of the thesis is presented, and extensions and future directions for the research are pointed out.  14  Chapter 2: Cross Directional Response Model  Chapter 2 Cross Directional Response Model  2.1 Introduction A primary goal of the paper making process is to maintain a uniform distribution of fibres across the sheet (CD profile). Fibre distribution is usually controlled by adjusting the slice, an opening at the bottom of the headbox, with variable width across the machine . Variability.in fibre distribution 1  across the sheet can be attributed to disturbances before the headbox, in the headbox, and on the wire. The opening at each location can be manipulated by a long screw and provides the most popular means to control the variations across the sheet. A n adjustment at one point in the slice will result in changes in the neighboring points due to the bending of the lip, and to maintain a flat profile requires a procedure to take into account the coupling between slice actuators. The main focus of this chapter is to model the actuator coupling and fibre flow using a model of the wet end fibre distribution mechanism. The model is based on two main physical principles. Dispersive wave theory explains the fluid flow at the surface of the slurry on the wire and Darcy's law is used as the physical principle behind mat formation and drainage on the wire side. The surface wave theory limits the application of the model to Fourdrinier machines, and twin wire formers and multiformer are not considered. The combination of these theories and a careful treatment of the boundary conditions at the slurry/mat interface will result in a model that predicts dry weight profile and the location of the dry line as well as its shape. This chapter is organized as follows. After the introduction, the dispersing wave model is derived, the drainage theory is discussed and mat formation is analyzed. The use of the model to explain wet end phenomena such as edge effect is demonstrated by simulation. Limitations of the model are explained.  2.2 Dispersive Wave Model Before arriving at the dispersive wave model for the stock flow on the wire due to an irregular slice opening, alternative theories for the surface profile modelling were examined. Intuition suggested 1  For a schematic of a Paper Machine refer to Figure (1.1) in chapter 1.  15  Chapter 2: Cross Directional Response Model Wave crests prediction based on Kelvin Ship theory  Wave x—position (<  ship heading)  Figure 2.2: Wave crests generated by a ship on a calm sea.  a surface wave model similar to that used for ship dynamics theory. The problem of a deflected slice on the wire is very similar to the surface dynamics of a moving stream at the presence of an obstacle. This problem has been discussed in naval architecture literature as Kelvin ship theory. Kelvin ship theory models the group of waves generated by a boat on a calm sea [56]. Figure (2.2) shows the wave crests generated in the wake of a boat. The generated waves in the wake of a ship are formed to make a fixed angle of 19.5° with the source of the disturbance for all hull shapes and speeds. The lack of variation in this angle was the main reason for avoiding Kelvin ship theory for application in paper machine modelling since it is observed that the angle varies significantly with consistency and other paper machine parameters.  The second model which was intuitively appealing was the spread  of a fluid on a solid or a fluid layer. The results obtained with the application of this theory were in contrast to observations at the wet end of the paper machine and the measurements of C D dry weight profile. Dispersive wave theory, on the other hand, could explain the physics of the problem and also, the results matched the observations of the dry weight profiles. The assumptions made for the model development are based on an approximation of the wet end formation process at the macroscopic level (It is a well known among papermakers that for a good formation at the wet end of the paper machine, a degree of turbulence is required in the headbox). The purpose of the development in this section is mainly to provide a better model for C D dry 16  Chapter 2: Cross Directional Response Model  weight control. The adequacy of assumptions will eventually be verified with the final results and their agreement with the industrial data. The major assumptions for the model are : 1.  Incompressible fluid; this assumption is a very good approximation.  2.  In viscid fluid; viscosity is not significant in stock.  3.  Irrotational flow; although turbulence is known to be present for good formation, the macroscopic behaviour is dominated by irrotational flow effects.  4.  Small amplitude waves; the amplitude of surface waves generated at the surface of the slurry compared to their wavelength is small.  5.  The response model is symmetric; although the consistency is not required to be constant in adjacent actuator zones, it is assumed to be constant locally in the headbox in the vicinity of an actuator bump.  Under these assumptions, the development of the dispersive wave model from the Laplace equation with the proper boundary conditions is carried out in the next section [18], [24].  2.2.1 Equation of Motion Consider the disturbance generated by the deflection of the slice lip at the surface of the slurry. This disturbance is characterized as a surface wave [56]. We assume that the flow characteristics in machine direction (MD) are negligible and the flow can be characterized as a two-dimensional flow across the sheet (CD) and vertically through the wire. This assumption does not require the basis weight or other properties of the sheet to be constant in M D . In Figure (2.3), a schematic of a two dimensional wave is shown. The amplitude is a function of C D position x and time t. The origin of the coordinate system is at the undisturbed surface of the liquid in such a way that a disturbance lies at both sides of the origin. When projected in three dimension, the time axis will be perpendicular to the x-axis out of the page, so each snap shot of the wave represents the shape of the surface of the slurry at different M D locations. Figure (2.3) shows the slurry with the surface wave on the wire. The x-axis shows C D and the y-axis, the thickness of the slurry. The positive x-direction points to the right. When a disturbance is generated at the surface, the wave fronts move in both the positive and negative x-directions. The  17  Chapter 2: Cross Directional Response Model  it  y  y=o  A Cross Direction (CD)= x  T  •BpHBpMMHMiHBiB  Figure 2.3: Schematic of wave generated by slice deflection on a porous bed (mat), wire represents the porous bed, and the shape of the wave at initial time is an approximation of the shape of the lip. The amplitude of the wave is shown as 7] and the mean value as h. Formulation of the small-amplitude surface wave equation is posed as a boundary value problem. This approach defines the physical governing equation in the region of interest with boundary conditions. The wave equation normally has an infinite number of solutions with the unique solution that is specific to the problem depending on the boundary conditions. The governing equation in the formulation of the C D response model is the Laplace equation, which occurs in many engineering disciplines and represents a potential field. In the case of small-amplitude waves, the field represents the velocity potential. The velocity potential is a function the gradient of which in any direction defines the fluid velocity vector in that direction. In order to find the velocity vector at any point requires only the tangent to the potential field at that point. Calculating the total variations (integral of velocity vector along a certain path) of the velocity field, it is only necessary to calculate the velocities at the beginning and end point of the path and the line integral to be independent of the path. The existence of a velocity potential is solely dependent on the irrotationality assumption made earlier (this assumption physically requires that there are no eddies or small scale turbulence in the flow). In mathematical terms, the curl of the velocity vector must be zero: V x « = 0  \ dy  dz J  \dz  dx y  J  \ dx  dy )  It may be shown that in order for the line integral of vector u to be independent of the integration path the above condition must hold and the gradient of a function <f> (potential function) is equal  18  Chapter 2: Cross Directional Response Model  Figure 2.4: Surface S surrounding the fluid. to velocity vector u.  = u  (2.2)  The incompressibility assumption reduces the conservation of mass law to the form presented in the following. Assume a closed surface S fixed in space. The conservation of mass within this surface, if it surrounds the incompressible moving fluid, implies the continuity equation, Figure (2.4). If the incompressibility assumption holds, there is no accumulation of matter inside the surface S and the outflow and the inflow should be equal.  \7^\over S  (2.3)  —0  The above equation simply states that the sum of the velocity gradients in any direction must be zero. The inviscid fluid assumption implies that the only external forces acting on the fluid inside S are the pressure force and a gravity force (all the viscous terms in the Navier-Stokes equation will be zero [24]). Newton's second law of motion (conservation of momentum) without the viscous terms is applied to fluid inside S in direction / (direction / , can be defined as the desired direction of flow). For example, C D can be considered one such direction. / as a vector in three dimensions is defined as I = x i + y j + z k):  (2.4) v  s  V  19  Chapter 2: Cross Directional Response Model  The left hand side of this equation represents the sum of momentum variations in direction / and the right hand side represents the sum of external forces (pressure and gravity). For the first term on the right the divergence theorem gives: j pl.ndS  = j  S  Vp.IdV  (2.5)  V  and collecting the terms to one side, the above equation will reduce to:  J {^m l VP' ) pL  +  9  dv  =°  ->  (2  6  v Since the above equation should hold for any volume V and any direction / , it results in the following Euler formula. Du  In the above equations ^  1  is the total derivative and the y is the gradient operator . Equation (2.7) 2  provides the pressure distribution given the velocity field. As mentioned earlier, an irrotational flow results in a velocity potential cf> and equation (2.3) expressed in terms of the velocity potential gives the Laplace equation :  d4> u ,  — = — ox °2  U  X  dx  2  A  Uy,  —  dy 2  x  d2<f>  d(f>  d(j>  ,  x  d 4> ,02^  c, , 2  ' dy  2  '  U,  Z  °  z  dz  2  „ ox C2.8)  With the assumption of two-dimensional flow, equation (2.8) reduces to: d d> 2  86 2  The equation of motion has been presented in its general form. In order to solve for the solution specific to the surface wave problem, boundary conditions need to be defined.  2.2.2 Boundary Conditions There are two types of boundary conditions, kinematic and dynamic. The kinematic boundary condition is derived from the continuity equation and states that the velocity of surface and local 2  For definitions and notations in the formulas refer to the Nomenclature.  20  Chapter 2: Cross Directional Response Model  velocity of the fluid should be the same at the interface so that no fluid particle leaves the interface. The dynamic boundary condition is found by applying Bernoulli's equation on the surface of the fluid. This condition explains the response of the free surface to pressure variations. Defining S as the boundary that moves with the fluid so that no particles leave S, we have the following equation:  g  = 0  (2.10)  This equation expresses the condition at the surface of the slurry (surface of the slurry is denned as 77) where the particles do not leave the surface. Defining S as below, and applying equation (2.10), we derive the kinematic boundary condition at the surface: S = r)(x,z,t) - y = 0 DS  drj  Or,  8n  (2.11)  for the two dimensional case and after linearization around the undisturbed surface, the above equation reduces to:  y = »?0M) d — - u - 0  (2-12)  v  on  y  y = 0  The dynamic boundary condition is found by applying the integral of equation (2.7) to the free surface interface between water and air. A characteristic of the free surface is that it cannot support pressure variations across the interface and must respond by changing its elevation to keep the pressure uniform. The integral of equation (2.7) is the familiar Bernoulli equation which relates the pressure field to the kinematics. Applying the Bernoulli equation to the free surface will give the dynamic boundary condition: dd)  1  -TTT +  77  at  2  , 1  +fi"?H  p„ - = Const p  on  y = n  , „. (2.13)  q in Bernoulli equation represents the velocity of fluid particles (in two dimensional case q  2  u  2  x  + Uy).  The dynamic boundary condition is a nonlinear equation.  =  Using the small amplitude  wave assumption, equation (2.13) is linearized [24] around the undisturbed fluid surface (y=0). The linearization procedure is shown below.  i +^  d  2 +  ^  +  7)l-=(^  +  5*  d fd(f>  + 21  ^ 7)lv=o + +  2+  1  2  p„V  (2.14)  Chapter 2: Cross Directional Response Model  For infinitesimally small 77 the velocities and pressures are small (p = 0 is gauge pressure) and all v  the second order terms can be neglected. After linearization equation (2.13) reduces to: +g =  0  V  on  y =0  (2.15)  At the interface of the slurry and the bed, two conditions should hold which are specific to this problem. The conditions are that the vertical velocity in the fluid and in the permeable bed should be the same, and the pressure across the fluid/bed interface should match. In the following, subscripts f and b are referring to the fluid and the bed. d<f> 777  =  y  u  dt  (2.16)  Pf = Pb  at  V-  -h  The continuity equation should hold for the bed, since the flow in the bed is also assumed to be incompressible.  The velocity in the bed is governed by Darcy's law which relates the vertical  velocity to the pressure across the bed. Vu = 0 K  u= K  p  (2-17)  -vp  represents the permeability of the mat. Substituting for the velocity in the continuity equation  with a similar term in Darcy's law results in a second Laplace equation for the pressure inside the bed. V P6 = 0  (2.18)  2  The boundary condition at the slurry/mat interface requires that both the velocity and pressure at the interface be the same. Equations (2.8) and (2.18) can now be solved using kinematic and dynamic boundary conditions of (2.12), (2.15) and (2.16).  2.2.3 Solution for exponentially damped waves Equations (2.8) and (2.18) are solved by the separation of variables technique with the boundary conditions (2.12)-(2.16).  A complex exponential surface disturbance of the following form is  considered, Since the solution considered represents the amplitude of the surface wave it should be real valued therefore, the assumption of the complex solution is merely a mathematical convenience. rj = A eJ ~^  for  (Nx  0  22  x >0  (2.19)  Chapter 2: Cross Directional Response Model  It may be verified that equation (2.19) is the solution for potential fields d> and pb defined by the real parts of  3  <f>(x,y) = p (x,t)  A  f  f  l  e  ,  M  = De ^ ^e^ - ^ N  b  A  h+  Nx  (2.20)  wt  = [A cosh. N(h + y) + B sinh N(h + y)].  w  The choice of solution as in equation (2.19) is presented in [24]. The continuity of pressure and velocity at the slurry/mat interface is the starting point. For a continuous pressure across the interface we have: Pb = Ps (2.21)  d<t>. P \y=-k  = ~Ps  Tt  Substituting for the pressure term in equation (2.20) gives: -jcopA  =D  (2.22)  The continuity of velocity components across the interface requires: d<f> _  lip dp  s  dy  p  dy  at y  —  —h (2.23)  or KB  B  A third equation is required to solve for A, B and D. Making use of dynamic boundary condition yields: ldd> g dt — = -. —(A 9 J  3  cosh Nh + B sinh Nh)  ^ ~ )  (2.24)  tanh Nh  (2.25)  Nx  wt  e  = Ai e ^ ~ ^ 3  Nx  wt  Substituting A and B from equations (2.22) and (2.23) gives:  D = pgAi< cosh Nh  The complete solution of the dispersive wave equation is presented in the Appendix  23  Chapter 2: Cross Directional Response Model  The linear kinematic boundary condition gives the dispersion equation when velocity term is substituted with the potential <f>: dr\ _ d(p &  dy  juAi  (2.26)  - AN sinh Nh + BN cosh Nh  The dispersion equation defines the relationship between spatial and time frequencies. This formula shows that if a group of waves starts from a single point they will eventually disperse since each time frequency corresponds to a different spatial frequency. The dispersion equation for damped dispersive wave on a sliding bed (such as mud) has been reported and similar relationships to those in [24] can be used here. After substitution for A\, A and B in terms of D, we reach the dispersion equation:  (  u  2  u)K \ n  tanh Nh  = gNy,anhNh ( t a n h N h - ^^ A J  (2.27) u  2  - gN tanh Nh = - j (  ) (gN -  tanh Nh)  The solution of dispersion equation results in the formula for N. In this case N is complex and can be written as N = k + jkj. The real part represents the wave number (spatial frequency) and the imaginary part determines the spatial damping. Based on the formulation of N, the real part of the equation (2.19) can be written as:  A e~  kdX  V  0  cos(kx - ut)  for  x>0  which is a damped progressive wave of period ^ and amplitude A e~  (2.28)  kiX  0  when the wave front is  travelling in the positive x-direction. This indicates that overall boundary conditions will determine the shape of the surface of the stock, and the interface boundary condition at the bed in specific, i.e. drainage, induces an energy dissipating mechanism that appears as the damping factor e~  kd x  , in  equation (2.28). As mentioned at the beginning of this chapter, a disturbance travels in both positive and negative x-directions. This implies an even function for the dispersive wave. When the wave front is moving in the negative x-direction a similar equation to that of equation (2.28) exists. The following formulation describes the surface elevation of a wave moving in the negative x-direction 24  Chapter 2: Cross Directional Response Model  (at t=0 wave crest at x=0). | a; | = —x ??_ = A  e  x < 0  - " l l cos(-kx f c  0  when  - uit)  x  ??_ = A e~  kdX  a  for  cos(kx + tot)  x < 0  for  (2.29)  x > 0  The last line in equation (2.29) simply gives the wave amplitude when the wave is travelling in negative x-direction and since the sign of x has been incorporated in the equation, only positive values of x are considered.  The overall solution of dispersive wave with damping is shown in  equation (2.30). j] = A e~ ^  cos(kx ^ ut)  kd  0  for  - b\ < x < &i (2.30)  where  bi > 0  The solution of dispersion equation (2.27) has been outlined in [75] and [24]. Results from [75] and [24] for the shallow water are: 2" 1/2  1 - LQ h/g  I  2  1 -  4  (1 - u h/g) 2  k d  gK^  U  (uK /!f) p  fc  (2.31) (1 -  2h  =  k h)(u;Ii /£ di  p  2h  The first line in the above formulation shows that A: is a function of different parameters, especially a change in depth has a great impact on spatial frequency. The speed of a dispersive wave ( phase velocity ) in the simple case of standing waves with no damping is defined as: CO  D  k =  v gh.  => c = — k  /  or  c = yjgh  (2.32)  Comparing the above formula with equation (2.31) shows that the speed of a travelling wave in the presence of drainage is similar to the speed of a travelling wave with addition of a correction term shown as the bracketed term in equation (2.31). The second line in equation (2.31) shows variations of damping factor  with depth (h) and process variables such as permeability constant (K ) and p  viscosity (ji) and initial damping factor k^. Due to the complexity of the variations of wave number with depth and uncertainties with other parameters such as permeability in Darcy's law, model parameters will be estimated by a system identification scheme in chapter (4). 25  Chapter 2: Cross Directional Response Model  2.3 Drainage Equation and Mat Formation The previous section explained variations in the surface of the slurry with time due to a bump at the slice lip. This section explains the mat formation model and is a summary of work on the Fourdrinier wire drainage by a number of researchers [10] [40] [44] [62]. The mat in this analysis is considered as a continuum, and the flow within the mat as laminar flow. Consider the wire as a fine meshed screen carrying a diluted suspension of water and fibre with consistency of c and s  height H. The initial filtrate velocity is zero, but builds up to vj as the suction underneath the wire is increased. We assume that the applied pressure across the wire is much larger than the static pressure of the fluid. The volume of the filtrate for a cross sectional area of 7 and time t is found to be  t  Vo = 7 / Vf dt. Some of the fibre is left on the top of the wire and forms a thin mat of consistency 0  c, m  which is larger than suspension consistency c . s  Smaller fibres (fines) pass through the pores  of the wire and form the filtrate with the consistency of c . w  c. s  Filtrate consistency c  w  is smaller than  As the slurry moves further down the machine in M D , the mat becomes thicker and the pores  of the wire and the mat become smaller providing a greater resistance to flow across the mat. We expect to see c  w  decrease and c  m  increase as the slurry moves closer to the dry line. The dry line  is the point at which the mat formation is finalized and there is no further deposition of fibres. The dry line can be easily identified on the wire since there is a significant change in the opacity of the sheet. Darcy's equation and a careful mass and volume balance across the sheet are the basis for the development of the mat formation model. The initial acceleration of the flow within the mat is assumed to be negligible so the governing equation reduces to Darcy's law. Ap =  (2.33)  This equation simply states that the pressure drop across the sheet is proportional to filtrate velocity, filtrate viscosity and a term that represents the mat resistance (jf-).  In the above equation L  m  and K  represent filter cake thickness and permeability of the filter cake respectively. Darcy's equation in this form is analogous to Ohm's law. A variation of Darcy's equation for a porous bed is normally used for paper machine drainage. This new representation is specifically adapted to the Fourdrinier 26  p  Chapter 2: Cross Directional Response Model  2  3  5 x 10  Mat Weight gr / s q . c m  Figure 2.5: Total drainage resistance vs. basis weight ( from reference [39] ). conditions. The formulation is: APA  r  h Rt  (2.34)  l  Rt — SFR  B  w  + R  w  In the above, v/ is the filtrate flow rate and AP is the pressure across the mat, fi is the filtrate viscosity, A is the cross section area of the mat and R  t  is the sum of wire drainage resistance R  and mat drainage resistance. SFR in equation (2.34) refers  w  to specific filtration resistance and B  w  is the total mat resistance. The total mat resistance Rt  is the mat weight. The specific filtration resistance is normally  determined experimentally by making use of Darcy's equation. A typical result of an experiment reported in [39] is shown in Figure (2.5). The slope of the plot of the resistance vs. mat weight indicates the specific drainage resistance. Although the total pressure drop across the mat is normally known through off-line analysis of the drainage elements of the wire, the problem formulation is still incomplete and a knowledge of the overall mat thickness is required. This could be achieved through mass and volumetric balance across the sheet. The following development is a summary of Meyer's work [62]. The decrease in the free surface depth of the slurry due to drainage is Ah and related to filtrate velocity. t Ah = - 7  27  Jvfdt  (2.35)  Chapter 2: Cross Directional Response Model  7 is the proportionality constant which represents the surface area and is considered to be 1 without loss of generality. A mass balance of the fibre on the same basis gives: L  t  m  J  - c )dl = - j  (c  m  (c - c )vjdt  s  0  s  w  (2.36)  0  Equation (2.36) states that the fibre loss in the suspension, neglecting the fines in the filtrate, will equal the amount of fibres retained in the mat. The negative sign indicates that the right hand side of the equation is the loss of fibres in the suspension. Considering the instantaneous weight formula, = J c  B (t) w  m  dl  (2.37)  o Equation (2.36) reduces to the following. i  B {t) w  - j  = c L (t) s  m  {c -  c )vj  s  w  dt  (2.38)  o As the mat builds up, the size of the pores within the mat decreases and the resistance increases. This translates to a variable filtrate consistency in M D . Han has suggested an empirical relation for a varying retention [39].  =  In equation (2.39) c  w  and c  Wo  C-Wo £ ^  (2.39)  are instantaneous (local) and initial white water (filtrate) consistencies,  3 is a constant and Bw represents the dry weight. Total instantaneous retention (first-pass retention), the ratio of the weight of the fibre retained at the top of the wire to the weight of total fibre [82], can be found as follows :  C-s  Cs  rt  0  = 1 - ^1  (2.40)  C  s  rt = 1 - (1 -  rto)e " _/?s  The last equation that completes the model formulation is the incremental weight formula [73]. AB  W  = c p Vf rt At s  s  28  (2.41)  Chapter 2: Cross Directional Response Model  In the above equation p is the suspension density and At is the time step equivalent to an increment s  in M D . The value of At is calculated by dividing the M D distance increment by machine speed (Ait = ^j-). This concludes the formulation of the drainage and mat formation. Knowing all the parameters in the above equations we can simulate the mat formation with a variable thickness. After a disturbance is generated by the slice lip, the resulting wave disperses to the sides while the mat forms underneath. Apart from the interaction at the mat/slurry interface, the drainage equations and the surface profile have been developed independently. At the dry line both models will be closely correlated. It is common to define the dry line as the point at which the sheet reaches a certain consistency [73]. Here we have adopted a different criterion. In our approach the dry line is the first point in C D at which the slurry depth becomes zero. The assumption made is that as soon as the surface wave dries at one point, the shape of the surface stays constant while drainage continues elsewhere in the process until the last point dries. To test the model, the parameters were taken from a fine paper machine and C D weight profile and the dry line shape and location were estimated. The simulation results were in good agreement with the observed profile in the mill.  In the next section, pressure profiles over different drainage elements are discussed.  These  pressure profiles will eventually be used in equation (2.34).  2.4 Pressure Profiles on Drainage Elements  Most of the dewatering on a Fourdrinier machine is carried out at the wet end by drainage elements.  The model includes pressure profiles and the necessary dewatering mechanism due to  drainage elements. The water removal at the wet end of a Fourdrinier machine can be divided into four mechanisms [10].  1.  Suction is probably the most important mechanism in dewatering of the sheet. The vacuum is generated by the relative motion of the fabric over the drainage element or by an active external device below the Fourdrinier.  29  Chapter 2: Cross Directional Response Model  2.  Inertial momentum generated by an inclined slice jet at the beginning of the wire is also a dewatering mechanism. When the jet hits the wire at an angle it tends to continue its way downward. This is not necessarily desirable due to its negative effect on formation.  3.  When the slurry travels over cylindrical formers the generated Centrifugal forces can cause dewatering.  4.  Gravity plays an important role in the dewatering of very slow machines but is only comparable with other dewatering mechanisms at a slow wire speed.  Suction and the centrifugal mechanisms are employed in the drainage elements normally used in the Fourdrinier wire, for example Table rolls; Foils; Flat boxes; Wet suction boxes and Suction couch. Modelling of table roll behaviour is still based on the original work of Burkhard and Wrist in [12]. They used an experimental paper machine and a piezoelectric sensor to measure the pressure profile over a table roll.  Their work covered the Fourdrinier machine with speeds up to 2500 fpm and  showed the effect of wire speed and mat resistance on the pressure profile. A typical drainage profile on a table roll is shown in Figure (2.6).  The T D C in the figure stands for Top Dead Centre of the  roll and is used as a reference point. Variation of the pressure profile with the wire speed is shown in Figure (2.7). The profile changes and covers a wider area as the speed increases.  With a higher  mat resistance the pressure profile tends to shrink towards T D C while the magnitude stays the same. Mathematical analysis of drainage at the table roll has been presented in [87]. Taylor's formulation relates the drainage rate over a table roll to wire speed, radius of the roll and a resistance term [87]. The complete formula is as follows: Ii  R r b P• U 2  2  t  30  3  (2.42)  Chapter 2: Cross Directional Response Model  Figure 2.7: Effect of wire speed on pressure profile over a table roll (from reference [12]). In the above formula, vj is the drainage rate, R is the resistance to flow, r b is the radius of the roll, t  U is the wire speed, p is the fluid viscosity, p is the fluid density and Ii' is a constant. Although the positive pressure could be advantageous at times, the positive section of the pressure profile over a table roll is generally an undesirable feature and other drainage elements have been considered. Pressure foils provide an alternative drainage element based on the suction mechanism and do not have a pronounced positive pressure upward at the early stage of drainage. Foils are more popular on modern Fourdrinier machines and come in different forms. Schematic of cross sectional views of some foils are shown in Figure (2.8) [13]. 31  Chapter 2: Cross Directional Response Model  o Contoured Foil Adjustable Foil  Figure 2.8: Schematic of some foils used in modern Fourdrinier machines (from [13] and [10].  The foil angle determines the developed suction [95], and there is an optimum angle associated with each foil width (length in cross sectional area). Typical values for the optimum angle are between 2 and 4 degrees with a common foil angle set to 3 ° . Foils have been reported to affect other sheet properties such as smoothness, formation and white water consistency. Improvement on first-pass retention is one of the main reasons for using foils over table rolls. With the exception of recent work on pressure profiles over gap former by Zhao et al and Green et #/[100], [36], there has not been a significant improvement in model development for pressure profiles since the original work of Taylor [88]. Taylor's formulation as reported in [10] is: ,  2a AP(x)  =  -pU  (2.43)  2  i + iMr) 32  ;  Chapter 2: Cross Directional Response Model  W (filtrate)  wire  Figure 2.9: Definition of the parameters used in Taylor's formulation ( from reference [10]).  25  vacuum  Eo o  -25  -50  -75 Figure 2.10: Typical pressure profiles over a foil and a foil with a vacuum pump, where parameters are defined as: a—  a'p k U p  k = Constant representing mat resistance p  p — Fluid density U — Wire speed A P = Suction over foil The above formula states that the developed suction over a foil is proportional to the dynamic head of the slurry on the wire with a time varying coefficient that is a function of the distance from the edge of the foil, and of the fluid properties. A schematic showing the important parameters is shown in Figure (2.9).  A typical suction profile over a foil is shown in Figure (2.10). The reduction of the positive  upward pressure segment at the beginning of the foil pressure profile is an important advantage of foils over table rolls. When the foil is also equipped with a vacuum device the pressure pulse will converge to the vacuum pressure.  For. the simulation, typical values of the above parameters taken 33  Chapter 2: Cross Directional Response Model  from a fine paper machine have been used.  2.5 Mat Compaction We now consider the mat compressibility under a vacuum pressure, which is another step towards a realistic wet end model. Drainage over a foil is caused by two phenomena, the suction generated by the drainage elements described in the previous section, and the squeezing effect due to deformation of the mat under the vacuum pressure pulse.  The mat can be considered as a sponge which is  compressed by the vacuum pressure resulting in an additional amount of water being squeezed out of the mat. This problem has been extensively addressed in the literature. Some recent reports can be found in [10] and [78]. Research on the compressibility of the mat has shown the importance of two factors. The mat becomes compressed due to the applied pressure and this compressibility changes with time [78]. The formulation presented here is based on [10]. Mat compaction is deduced from static compression tests in a cylindrical container on top of a screen. The mat compression is studied under different loads, and the mat thickness and the pressure are recorded. The test is carried out in an equilibrium state, meaning that only the permanent (plastic) deformation is considered. By correlating the mat concentration and the applied pressure, the following expression is found.  c  In the above equation, c  m  m  - M  AP  (2.45)  Nl  is the mat concentration and AP is the applied pressure. M and N/ are  constants. Since the mat concentration is a non—zero value as AP —*• 0, equation (2.45) will not be valid at low pressures. To remedy this problem, a constant term has been added. The modified mat compression equation is:  c  m  c  m o  = c  + M AP  Nl  mo  (2.46)  in the above equation is the mat concentration in the uncompressed state. The values of M and  Nj are found through experiment and have been reported in the literature. A number of researchers have investigated the variation of the constants M and N/ due to bleaching and Kappa number. In order to incorporate the effect of mat compaction in the wet end model, one has to resort to Darcy's 34  Chapter 2: Cross Directional Response Model  o  correct alignment  incorrect alignment  Figure 2.11: Headbox induced edge waves due to misalignment of cheeking elements, equation. Equation (2.34) is considered. The specific drainage resistance SFR is related to the mat concentration through the following expression: 1  SFR  K  p  ^ p Cm  (2.47)  is the permeability constant. Combining Equations (2.47) and (2.46) and using the appropriate  coefficients, we can incorporate the effect of mat compression in the model.  2.6 Superposition and Response Model at the Edge The  development of the surface wave on the slurry was based on a linearized theory. The  assumption of small wave amplitude allowed the linearization of the nonlinear boundary conditions. In a linearized theory, solutions can be superimposed. The response of a single actuator can be extended to overall shape of the slice lip in order to achieve a complete C D response model. Superposition is also used in developing C D response at the edge. Edge waves can be attributed to different sources. A major cause of waves on the Fourdrinier is an incorrect installation of cheeking pieces inside the headbox [67]. A cheeking piece as defined in [82] is: A removable dam within a headbox which forms part of the side of the box and seals the end of the slice opening. A schematic showing the proper installation of a cheeking piece, and an incorrect installation is given in Figure (2.11) [67].  The  misalignment of the cheeking piece can initiate an edge wave in the headbox. Only waves generated by the slice lip and its reflection at the edge are considered. On a Fourdrinier machine the fibre  35  Chapter 2: Cross Directional Response Model  Dry Line  it Main and 2nd Bleed properly adjusted  No Bleed  Figure 2.12: Proper use of edge bleeds reduces the edge wave effect (from reference [81]).  suspension is contained at the edge by a deckle board. Deckle boards are attached to the headbox and sealed to the wire by flexible rubber strips. With the introduction of flexible wires it is possible to curl the wire at the edge to contain the stock. Implementation of either or both of these methods will reflect the lateral surface wave and unless an energy dissipating mechanism is used, the edge wave will persist. A n energy dissipation method to dampen the edge waves has been mentioned in [72]. A common solution to reduce the edge effect is to use an edge bleed. Edge bleeds are nozzles mounted at the side of the wire and their function is to inject a diluted fibre mixture onto the sheet to improve edge profiles. The lateral position of the bleed depends on the shape of the edge wave. The effect of an improper and a proper use of an edge bleed is shown in Figure (2.12) [81]. The void area at both sides of the sheet caused by edge waves is very weak and normally trimmed. A n estimate of the location of the edge wave is beneficial in adjusting the position and the discharge of bleeds.  In  the formulation of the edge wave, equation (2.28) is considered since the only factor changed at the edge is variations in the surface profile. The dispersive wave model as in equation (2.28) describes a plane progressive wave moving in the positive ^-direction with a spatial frequency (wave number) of k and time frequency of u. At the edge, assuming a frictionless impact, the wave is reflected with the same amplitude and velocity in the negative x-direction and a phase shift of 180 degrees. The reflected wave will have the same form as the dispersive wave travelling in the negative x-direction. The interference of these waves generates a standing wave pattern at the edges.  -hi \x\  f]r - A e 0  cos (kx + ojt)  for  x < 0  (2.48)  Applying the superposition theorem and replacing for |x| = -x when x is negative, the surface profile 36  Chapter 2: Cross Directional Response Model  at the edge can be found as: Vedge = V + ledge  Vr,  = A e~ °  cos (kx  kdX  0  $1 —  - u)t ) + A e~  COS (kx! +  kdXl  0  0  0  Uti),  ~f~ ^^edge  (2.49) X\ — x 4/- Ax g ed  e  Vedge  = A e~ " ° k  x  0  A  cos (kx  -  cos (kx  + ojt + k A t  0  e~ ° kdX  0  e  e  U)t )+ 0  e  e d g e  + uAx )  *  edge  - ^ ' <" kd  x d  e  In the above equation, to and %o refer to the time and distance that an incoming wave meets the reflected wave from origin in positive x-direction. The parameter At / ec  time of the wave to the edge and Ax /  ge  refers to the travelling  is the distance between the centre line of the generated  ec ge  disturbance and the edge. Parameters t and x refer to the time and distance from the edge until the e  e  reflected wave reaches the incoming wave. If the reflection at the edge is not frictionless, the surface profile can still be found using the above analysis. In this case, the reflected wave is multiplied by a reflection coefficient a, 0<a<l, representing the percentage of energy loss at the impact. The surface profile of the reflected wave with energy dissipation is: ledge = V + Tf , r  Vedge  = A e~ °  COS (kx  kdX  0  0  aA e~ '  - Ljt )+  cos (kx  kdX  edge  (2.50)  0  e  + ojt ). e  In the above formulation an alternative representation of equation (2.49) is used. Here instead of using the amplitude at the origin Ao we have used the amplitude of the wave at the edge. The next section presents the simulation results of the developed model for different operating conditions of a paper machine.  2.7 Simulation Results In this section simulation results for different operating conditions of a paper machine are presented.  Model parameters represent the operation of a fine paper machine with five different  scenarios.  The scenarios considered are:  1.  C D response of a slice bump for a heavy grade sheet. 37  Chapter 2: Cross Directional Response Model  Simulation scenarios  Machine speed  Headbox consistency  (m/min)  (%)  Heavy grade  400  0.9  17.0  Light grade  700  0.4  10.0  Edge response  400 & 700  0.9 & 0.4  17.0 & 10.0  Drainage scenario  700.0  0.4  10.0  Superposition  800  0.4  9.0  Slice opening (mm)  Table 2.1 Machine parameters used for simulation of 5 scenarios mentioned in the list. 2.  C D response of a slice bump for a light grade sheet.  3.  C D response at the edges.  4.  Effect of drainage profiles on the shape of response and dry line location.  5.  C D response for superposition of bumps at the slice.  In the simulation results the location of dry line is an average M D position of the dry line since the dry region represents a band rather than a single line. The amplitude of the bump in heavy grade case is 1.5 mm and in the case of light grade paper is 0.8 mm. Figure (2.13) shows the model response for a heavy grade of paper. The M-shape response of a slice bump on a heavy grade of paper is often encountered among paper machine operators. The first figure on the left is the surface profile of the slurry and the mat being formed from underneath while the mat is shown next. The shape of the dry line and surface profiles at different M D locations are shown in the second row. The profile averaged across the sheet (MD profile) and the average pressure across the mat are shown in the last row. Figure (2.14) shows the same element as in the previous figure when the parameters have been adjusted for a light grade paper. The difference in the C D response shape and the increase in the average pressure across the sheet due to a higher machine speed (wire speed) for a light grade sheet is noticeable. The edge effect has been shown in Figure (2.15). The first column shows the mat when the wave generated by the slice bump is reflected from the edge of the wire for a heavy and a light grade sheet with no energy loss at the impact. The last row shows the change when the wave is reflected with 40% energy loss at the edge for a heavy grade sheet. The second column represents the C D response with the location of the bump shown. The 38  Chapter 2: Cross Directional Response Model  magnitude of the bump is not to scale. Figure (2.16) illustrates the effect of different arrangements of the drainage elements on a paper machine and the corresponding results on the shape of the C D response and the location of the dry line. The drainage elements used are based on the models presented in section (2.4). Three different foil models have been used in the simulation, namely, the 3-degree foil, vacuum foils and vacuum boxes. The case shown in the left column of Figure (2.16) has 12 3-degree foil banks and 2 vacuum foils and 5 vacuum boxes. The number of 3-degree foils has been increased to 15 and vacuum foils to 3, while the number of vacuum boxes has been reduced to one in the case shown in the right column. The shift in the mean value of the dry line is obvious in Figures (c) and (d), and the change in the shape of the C D response can be seen in Figures (a) and (b). In order to extend the C D response model of a single bump to the overall shape of the lip, superposition of surface profiles should hold. The application of the superposition principle for the surface profiles has been illustrated in Figure (2.17), where two bumps with different magnitudes and bumps with the same magnitude are simulated.  2.8 Summary In this chapter the main ideas of the C D response model were developed. The model development started from continuity equations and resulted in a dispersive wave as a solution of the Laplace equation. Darcy's equation was critical in the development of the drainage equations and the boundary conditions of the wave model. The combined surface wave theory and drainage principles provided a model for cross direction response that is simple enough to be used for C D control purpose and rigorous enough not to compromise any of the important physical principles involved.  Some of  the operating conditions on a paper machine such as edge effect and modification of the drainage elements can be simulated with the model. The developed model requires less parameters to estimate when used for control purpose, and the parameters represent physical variables. The model has it limitations which will be briefly presented. One of the limitations of the developed model is the exclusion of the rush/drag phenomenon on the wire. On a real paper machine there is a speed difference of 4-7 fpm between the slurry and the wire. The reference in this terminology is with respect to the wire, when the slurry is faster 39  Chapter 2: Cross Directional Response Model  than the wire the case is referred to as rush and when the slurry is slower than the wire, the case is called drag. The rush/drag phenomenon is crucial to a good formation of the sheet and strongly related to drainage on the wire. This factor can only be incorporated in a two dimensional wave equation with a varying drainage parameter. Since the response model is assumed to be symmetric it cannot model disturbances that will affect response symmetry. An asymmetric response results from different sources. The asymmetry in response may be attributed to measurement noise, although localized consistency variations in the vicinity of a deflected slice lip can also impact response shape. The other element not considered in the model is the variation of the parameters in the surface wave equation due to changes in depth. There is a point beyond which the condition of the potential field will not necessarily hold. A variable depth results in a variable wave number or spatial frequency, k, in equation (2.28). Since the spatial frequency is a variable, the dispersion equation obtained by solving the Laplace equation f.p.mcannot be used. This problem is solved by employing the technique presented in chapter 4. The next chapter deals with the data processing scheme in scanners. This step is important since most of the industrial data is processed and goes through a similar procedure before being used.  40  Chapter 2:  Cross Directional Response Model  200  CD Position  Q o  200  M D Position  CD Position  Shape of dry line  x 10~  3  Q o  M D Position  Surface profile of slurry 1.6 1.4 1.2 1 0.8 0.6  40 C D Position  Mat thickness without compaction effect  50  100 M D Position  150  0  20  40 CD Position  60  80  Average applied pressure across mat  200  50  100 M D Position  150  200  Figure 2.13: CD basis weight response to an adjustment of a slice lip actuator for a heavy grade paper (All CD and MD position units are in samples. Dry line is an average location of the region).  41  Chapter 2:  Average applied pressure across mat  Mat thickness without compaction effect  20  40  60 80 M D Position  Cross Directional Response Model  600  100 120  40  60 80 M D Position  Figure 2.14: CD basis weight response to an adjustment of a slice lip actuator for a light grade paper (All CD and MD position units are in samples. Dry line is an average location of the region).  42  Chapter 2:  Cross Directional Response Model  CD response at the edge for a heavy grade paper  CD Position CD response at the edge for a light grade paper 0.5  C D Position  CD response with 40% energy dissipation at the edge 0.55 i • • • • 1—i  C D Position  Figure 2.15: Basis weight response at the edge for different grades of paper (All CD and MD position units are in samples).  43  Chapter 2:  Cross Directional Response Model  (a)  Average mat profile with pressure pulse drainage  0  20  Average mat profile with pressure pulse drainage  40 60 80 100 120 (c) M D Position  0  20  40 60 80 100 120 (d) M D Position  Figure 2.16: Effect of changing the drainage elements on a Fourdrinier machine and the corresponding shift of the dry line location (All CD and MD position units are in samples).  44  Chapter 2:  (a)  Cross Directional Response Model  (d)  200  C D Position  0 0  M D Position  CD Position  (b)  0  0  M D Position  0  0  M D Position  (e)  200  C D Position  0  0  MD Position  CD Position  Superposition of adjacent slice bumps  20  40 60 (c) C D Position  Superposition of slice bumps with different amplitudes  80  20  40 60 (f) C D Position  80  100  Figure 2.17: Superposition of surface profiles in two different cases. Figures (a), (b) and (c) refer to two adjacent bumps with the same magnitude and the cases (d), (e) and (0 with different amplitudes (All CD and MD position units are in samples).  45  Chapter 3: Testing and Experimental Methods  Chapter 3 Testing and Experimental Methods  3.1 Introduction Estimation of the parameters of the developed cross directional response model requires accurate measurements of the slice lip actuator settings and the basis weight and moisture profiles. In this chapter, a careful analysis of the data processing in the basis weight sensor is presented. This analysis is a necessary step before the slice lip bump test is carried out. The knowledge of data processing method in the measurement system is crucial to the parameter estimation scheme presented in chapter 4.  The basis weight sensor is a scanning gauge based on the principle of nuclear backscatter. It  correlates the weight of the sheet to the intensity of the scattered beam from a nuclear source (beta rays or high energy electrons for basis weight) measured across the sheet. Normally, the basis weight signal is heavily filtered in time due to a high noise-to-signal ratio. Sampling rate reduction of the measured signal also takes place. Heavy filtering and sampling rate reduction may distort the signal and mask the real measurable weight variations if they are not carried out properly. In this analysis, an independent data collection system is used and multirate system theory is applied to analyze the basis weight signal. This chapter is organized as follows. Preliminaries and a schematic of modules in the scanner is presented first. Second, a brief introduction to multirate system theory is explained. Next, the calibration method in the Measurex system along with the calibration results are introduced. Finally, a comparison of the independent data analysis with the readings from the scanner is made and some conclusions are drawn.  3.1.1 Purpose of the Experiment Rapid pressure variations in the headbox and the approach system will result in basis weight variations that are not always detected by the scanner and therefore are not relayed to the operator console. This lack of detection may be due to either a scanner with low resolution or to the existence of heavy filtering on the scanned data. The purpose of this chapter is to investigate the data processing 46  Chapter 3: Testing and Experimental Methods  scheme of the Measurex scanner and compare the scanner output with that of an independent data collection method. This investigation was critical for the correct analysis of C D bump tests performed later in the study.  3.1.2 Scope of the Experiment The basis weight signal used for the analysis varies between 0 — 8 volts and is calibrated using a standard scanner calibration scheme. The moisture signal was not used for the calibration. Profiles provided by the Measurex system were used for comparison with profiles calculated from the independent data source. Measurex scanner reference manuals and personal communication with the Measurex representative in the mill were the secondary sources of information.  3.1.3 Data Collection Scheme The data collection device was a data acquisition board with four differential inputs and a maximum sampling rate of 100 ksamples/sec installed into a 386-Toshiba portable computer. The sensor was a Model 2002 Measurex scanner on a fine paper machine.  A system diagram of the  Measurex scanner is presented in Figure (3.18) at point A. The voltage signal is sampled at the V F C (voltage to frequency converter) which accumulates the counts every 50 ms providing a sampling rate of 20 Hz. After calibration and some error corrections, the data is sent to the control computer via a serial link. The basis weight signal was initially lowpass filtered with a first-order passive filter and a cut off frequency of 350 Hz and then sampled at a rate of 1 KHz by the independent data acquisition system.  3.1.4 Background The scanner time resolution is limited to 50 ms which translates to an approximate spatial wavelength of 8 mm for a sheet of 5 meters wide and a scanning time of 30 sees. This wavelength is well beyond the slice lip control capabilities for CD. High resolution C D profile and high frequency M D variations require a higher sampling frequency which in turn means more storage space. In order to deal with the large storage requirement of a higher sampling rate, an array of multirate techniques have been developed [19]. A n important function of the multirate system is to change the sampling 47  Chapter 3: Testing and Experimental Methods  to A/D & laptop amplitude detector P C B i WW  weak current 0 - 8 volts  -350 VDC volt inverter  centre electrode  Ion chamber assembly Xenon Argon  count  compensator / noise piece  n  n  n  BW calculation done here  j sheet  serial link to computer control centre  Source Assembly Kr  Figure 3.18: A system diagram for Measurex scanner, rate of a signal while minimizing distortion. There are two approaches to obtaining a sampled signal. If the signal has significant energy up to a frequency fc, an analog filter (antialiasing filter) can be used to lowpass filter the signal up to this frequency. The signal can then be sampled and digitized. This approach requires a high order filter with a small transition band. The second approach is to use a low order antialiasing filter with a higher cutoff frequency of only half the sampling frequency. However, if the signal has been oversampled by at least twice the required frequency, it may then be digitally lowpass filtered and subsampled. The second approach with a simpler filter design has been used in the experiment.  3.2 Multirate Filtering In many signal processing applications the sampling rate of the discrete sequence needs to be changed. When the sampling rate has to be increased, the process is referred to as interpolation or 48  Chapter 3: Testing and Experimental Methods  Sample _ f rate  s  M-rold decimator  Sample rate  =  f  s  /  Sample _ f rate  M  s  downsampling byM  upsampling by L  Sample _ f * ^ rate s  L-fold expander  L o w p a s s  filterjng  byH'(z)  Figure 3.19: Multirate system building blocks, M-fold decimator L-fold expander. upsampling. In brief, we have a sequence and need points between the samples. One method of obtaining the intersample points is by simple averaging between the two samples and hence, the name interpolation. In the application at hand, the Measurex scanner data set, a sequence of about 600 data points is reduced to a sequence of 83 points. This problem can be viewed as a data compression process or sampling rate reduction (down sampling). This section covers the basic principles of multirate filtering theory. The multirate system theory is based on two main building blocks, an M-fold decimator and an L-fold expander, as shown in Figure (3.19). In order to subsample a sequence by a factor M (M-fold decimator), the sequence should be lowpass filtered before choosing every M sample of the sequence. The lowpass filtering is required to avoid any distortion or Aliasing.  The characteristics of the lowpass filters used in  multirate systems play an important role in the quality of the end results, and proper filter design is crucial in any multirate filtering application to avoid amplitude and phase distortion. Multirate filtering has found application in many branches of signal processing and communication including implementation of narrowband lowpass filters and data and image compression for telephony and fax applications.  3.2.1  Aliasing In order to avoid any distortion of the signal due to subsampling, the signal has to be band limited 49  Chapter 3: Testing and Experimental Methods  Original lowpass signal 1 rv-  Signal subsampled by 2  10 Samples Spectrum of original signal  Signal subsampled by 4  2 4 Samples  Samples Spectrum subsampled by 2  »-  Spectrum subsampled by 4 1 New  Nyq.  1 frequency  l 6 I I  \ \ \  1 ' 1  ^ 1 /  /  /  /  \ 1 \ 1 \ 1 \ 1  /  0  20 Frequency Hz  0  10 Frequency Hz  0  /  /  5 Frequency Hz -  \ 10  Figure 3.20: Sampling rate reduction without preftltering results in aliasing, before sampling rate reduction. This translates to lowpass filtering the signal with a corner frequency that guarantees no overlap of the frequency responses (antialiasing filter). The effect of subsampling without lowpass filtering is shown in Figure (3.20). A typical lowpass signal is considered. The second row shows the magnitude response of the signal when the sampling rate has been reduced by 1/2 and 1/4 respectively. The subsampled versions of the signal have maintained the same energy as the original signal. Without appropriate lowpass filtering, the aliasing increases as the sampling rate is further reduced from 1/2 to 1/4 due to increased overlapping of the magnitude response replicas.  3.2.2 Filter Design It is evident that to avoid aliasing, the signal must be lowpass filtered before subsampling. Every filtering process introduces some form of delay in the output. This delay is referred to as phase shift and can be adjusted with a careful design of the filter. If we consider a lowpass filter with a sinusoidal input within the pass band, the output will be a sine wave of the same frequency and the 50  Chapter 3: Testing and Experimental Methods  PL M  Nyquist frequency  0 w C  w M  Figure 3.21: Minimum requirement for a typical decimation filter when subsampling by M.  corresponding phase shift (delay with respect to input). If the amount of phase shift increases linearly when the input frequency of the signal is changed, the filter is said to have linear phase property. This is a very important attribute of the filter since the output shift with respect to the input can be compensated after the filtering by a time shift. In the context of scanner data analysis, this property plays an important role in the final C D profile since the phase shift translates to misalignment of the databoxes and eventual distortion of the profile. As a general guide line, a proper decimation filter should have two main properties:  •  Linear phase property.  •  Proper bandwidth cutoff to prevent aliasing.  The minimum bandwidth requirements for a decimation filter to avoid aliasing are shown in Figure (3.21).  More stringent requirements in the passband and the stopband can be easily incorporated in  the design of decimation filters. The filter implementation scheme is an important aspect of digital filter design. Finite impulse response (FIR) and infinite impulse response (IIR) are two methods of implementing lowpass decimation filters. The design of FIR and IIR filters in the context of multirate filtering can be found in [19]. The description of data processing scheme in the Measurex scanner is given in the appendix.  51  Chapter 3: Testing and Experimental Methods  3.3 CD Profile Calculation The nature of a traversing scanner results in zigzag measurements of the sheet properties. The measurement is decomposed into three components: machine direction variations (MD), cross direction variations (CD) and random variations. The raw measurement needs to be decomposed into C D and M D components to enable proper control action. The standard approach for decomposition of C D and M D variations is exponential filtering, a moving average filter over points in machine direction at a C D location. The main assumption in applying this algorithm is that the C D profile is not changing rapidly in comparison to M D variations and can be represented by deviations from a mean value of basis weight or moisture taken over one scan. The mean value is referred to as the scan average and is an M D variable. The output of the exponential filter is the C D value for the specific databox. The exponential filter is represented by: y (t)  = (l-a)y (t-l)  n  y {t) n  n  + ax (t)  (3.51)  n  is the estimated profile, n refers to the C D location or databox and t is the scan number.  The filter pole is at (1 — a) and the raw measurement is represented by x {t). n  At the time of the  experiment, the value of the filter pole was set to 0.8 whenever the exponential filter was used to extract a C D profile from raw measurements. This setting is typical of industrial usage. A second approach to C D and M D profile estimation is the method proposed by Natarajan. et al [27]. This method uses Lindeborg's model for filtering the raw moisture measurements [59]: D(n, k) = MR + p  n  + u (l k  + Bp ) n  + 0  where D(n,k)= M R =  Moisture content at position n and time k. Reference/mean value of moisture in each scan.  p =  Profile deviations from mean in C D .  Ufc =  M D variation at time k.  B =  Nonlinear interaction between C D / M D .  n  (3.52)  The assumption by Lindeborg that "CD variations are slow compared to M D and the M D disturbance is the same all across the sheet" is similar to that in exponential filtering. Based on 52  Chapter 3: Testing and Experimental Methods  equation (3.52), p  n  and B are estimated using a least-squares parameter estimator in C D , and Uk  is estimated with a Kalman filter formulation [27].  For better estimation results, a least-squares  algorithm with resetting and forgetting factor may be used. The same formulation is used for basis weight estimation by setting B equal to zero. This is justified since the bilinear term in the moisture formulation is a result of the drying process and is not expected to be present in the basis weight case. Other techniques for decomposition of C D and M D variations from diagonal measurements have been described by Tyler [89], Goodwin [34] and Nesic [66].  Tyler et al employ an input output  model to estimate basis weight variations taking into account the dynamics of variations in C D . The scanning process is modeled as a linear, periodic time-varying operator with period N combined with the input/output model. The estimation of states is carried out by a periodic Kalman Filter, and the Riccatti equation reached in the optimization process is solved with the lifting method. Goodwin et al have treated the zigzag measurements as a single dimensional stream of data and use a periodic Kalman Filter to estimate the variations. The model for C D / M D is a linear spatial system. Nesic has used a two-dimensional discrete wavelet transform with soft thresholding to extract C D / M D components from noisy measurement.  3.4 Calibration of Basis Weight Sensor The calibration experiment was carried out on a stationary paper sample cut out of a standard sheet with a nominal basis weight of 72 lb/3000sqft. Masking tape and wooden rulers were used to make two artificial edges at some distance along the sheet. The sheet was held stationary and the scanner traversed the sheet several times while the data acquisition and the Measurex system were collecting basis weight data at the same time.  3.4.1 Test Specifications The calibration test was performed while the machine was down for regular maintenance. A schematic of the sample used for the experiment with the location of the artificial edges and the dimensions of the sheet is shown in Figure (3.22).  Three different methods were used to measure  the basis weight of the sheet. The scanner traversed the sheet 10 times while the sheet was held 53  Chapter 3: Testing and Experimental Methods  All dimensions in millimeter edge  „ 1015  J  punched samples  REF  2345  J  L_2-52  Figure 3.22: A schematic of the sample used for calibration. The circles show the position of puncher for manual measurement of the basis weight. stationary underneath it, and basis weight data were collected simultaneously by the Measurex system and the data acquisition system separately. A manual measurement of the basis weight across the sheet was also carried out after the scanning portion of the experiment was completed. The data collected from the manual measurement of the sheet is tabulated in Table (3.2). The dimensions of the puncher that was used to extract samples were as follows: 1.  Diameter of the puncher = 4.5 in.  2.  Area of the punch = 15.90 sq in  The nominal grade of the sample was 72 lb (lb/3000sqft). The data acquisition sampling rate was set to 1 KHz, and the corner frequency of the antialiasing filter to 350 Hz.  Samples 1—11 were taken  between the edges scanned by the sensor, and the remainder were outside that region.  3.4.2 Analysis Scheme The data processing in the Measurex system is a decimation process, and, in every decimation process, the key element is the lowpass filter. The aim of our analysis was to focus on the effect of different filter characteristics on the estimated C D profile. Three different FIR filters were chosen for this purpose. The first choice was a rectangular window since this is the filter used by the Measurex system. The other filters were a Hamming and a Hanning window, each commonly used in the signal processing application. The collected basis weight data were first split into air-scan and back-to-front 54  Chapter 3: Testing and Experimental Methods  Basis weight / sample area  Distance from Reference  (gr/cm2)  (cm)  1  1.2205  7.5  2  1.2199  21.5  3  1.2191  40.0  4  1.2225  54.0  5  1.2144  80.0  6  1.2209  112.0  7  1.2117  134.0  8  1.2137  155.5  9  1.2145  174.5  10  1.2150  195.5  11  1.2171  226.0  12  1.2215  248.6  13  '1.2189  260.0  14  1.2189  271.4  15  1.2168  282.8  16  1.2195  294.2  17  1.2170  305.6  18  1.2217  317.0  19  1.2221  328.4  Sample  Table 3.2: Results of manual measurements of basis weight across the sample sheet. and front-to-back half-scans and each sequence was filtered and subsampled. The final sequence obtained had the same number of samples as in the Measurex profile. The filtered sequences were calibrated, and the C D profile was calculated using an exponential filtering scheme. The profiles were finally compared with the Measurex profile. The basis weight of the sample sheet was also manually measured at different locations between the edges, and the results of the filtered profiles and the hand measurements of the basis weight were compared. Next, the results are presented.  3.5 Results Figure (3.23) shows the raw data from the scanner. 55  The scanner is off for the first 30000  Chapter 3: Testing and Experimental Methods  samples. It is then turned on and a standardization test is performed. Samples with a mean value of approximately 4.25 volts are the off-sheet measurements.  As the scanner moved on-sheet, it  generated a spike of approximately 1.6 volts. The presence of the spike could not be correlated with regular operation of the scanner and its cause was not clear. One possible explanation could be a forward differentiation process that generates a spike in the presence of a sudden change of signal (a step change). Spikes are also caused by the artificial edges introduced as marking points. Each complete scan consists of:  1.  Scanner coming on-sheet at the back side of the machine.  2.  Measuring the edge at one side and scanning the sheet till it detects the edge at the other side.  3.  Scanner changing direction and moving from front to back of the machine, passing over the edge for a second time.  4.  Continuation of scanning and passing over the edge at the back end of the machine and finally going off-sheet.  Figure (3.24) shows the decimated sequence when a rectangular window is used. At the top of the figure, all the 10 scans are overlaid. The figure with overlaid scans represents the true measured C D profile, since the sheet was held stationary during the experiment. The characteristics of the filter are shown in the middle row and the C D profile obtained by exponential filtering is shown at the bottom of the figure. The bounds shown on the bottom row are calculated based on the variability of each point in M D position. This variation indicates the repeatability of the measured points. Figures (3.25) and (3.26) show the decimation of the data using a Hamming and a Hanning window. The second row shows the filter characteristics in time and frequency domains, and the bottom row shows the C D profile with the 2a bounds. Profiles obtained by different filtering techniques are compared with the Measurex C D profile in Figure (3.27). The solid line represents the Measurex profile, and the dash line is the profile obtained with different windows. The 2a bounds are calculated as before and are an indication of the repeatability of the measurement points. The results of this comparison is summarized in Table (3.3). The values in the table are the sum of squared error. The increased variance of the Hanning and Hamming windows is a result of greater frequency bandwidth. 56  Chapter 3: Testing and Experimental Methods Raw calibration data before processing 4.5 |  1  1  0.2  0.4  1  1  1  1  1.5 1 -  0.5 -  0I 0  i l l  1  1 _  0.6  1  i  i  0.8  1  1.2  Samples  i 1.4  >  i  i  1.6  1.8  I 2 |Q5  Figure 3.23: Calibration data collected by the data acquisition system before processing.  Profile calculation scheme  Harming window  Hamming window  Rectangular window  S (error from mean)"2  9.20  5.61  3.65  Table 3.3 Comparison of the Measurex profile with different profiles in the squared error sense. In order to compare the results of obtained profiles with the hand measurements of the basis weight, the profiles had to be further decimated. The length of the hand measured profile was 11 samples, and the profiles were subsampled to a length of 11 samples. The results are shown in Figure (3.28). The solid line is the hand measured profile, and the dash line is the decimated profile. The bounds are based on the variations of hand measured samples. The results of this comparison are summarized in Table (3.4), and the values in the table are the sum of the squared error from the mean value. When the signal is decimated to 11 samples, it does not represent the variations of the signal and should only be used for calculation of mean basis weight. The results in Table (3.4) show  57  Chapter 3: Testing and Experimental Methods  Profile calculation i . scheme S (error from mean)"2  Hamming window 1.41  Measurex  0.97  Hanning  Rectangular  window  window  0.90  0.78  Table 3.4 Comparison of hand measurement of basis weight with different profiles in the squared error sense, that the boxcar filter has the smallest variance because of smaller bandwidth.  3.6 Summary In order to establish the link between filtering techniques and C D profiles, we need to concentrate on Figures (3.27) and (3.28). A close match was expected between a C D profile based on a rectangular window and the Measurex profile, since this is the window used. The existence of error indicates that filtering by other elements in the scanner is taking place. An interesting observation in that although there is a good match between the two profiles in one scan direction (back-to-front) there is a mismatch that can be viewed as a shift in the other direction (front-to-back). This may be due to speed variations of the scanner at the front side of the machine. A comparison of Figures (3.24), (3.25) and (3.26) shows that the Hamming-and Hanning windows have better attenuation at higher frequencies than the rectangular window. This difference in attenuation characteristics will increase aliasing effects in the case of a rectangular window or the Measurex system at high frequencies. This effect cannot be observed from the results, primarily because of the sampling rate reduction. In the case of the Measurex system, the sampling rate reduces from 20 Hz to approximately 3 Hz (subsampling by 7) and the Nyquist frequency will be at 1.5 Hz. The aliasing that is caused by difference in filtering characteristics is mainly in the high frequency range and will not contribute significantly to the C D profile. There is also a second stage of filtering by the exponential filter. The comparison of the hand measurements and each of the profiles showed that the hand measurements are close in some portions of the profile, but that the resolution was not fine enough to display more than the trend. The closest profile to the hand measurements in the squared error sense was the rectangular window profile. A profile formed with such a coarse resolution would be best suited for determining the mean basis weight of the sheet.  58  Chapter 3: Testing and Experimental Methods  In summary, the following conclusions are drawn from the calibration: 1.  Based on the current resolution requirements for the C D control, the Measurex results conformed with those of the independent measurements.  2.  For improved calibration and verification of the Measurex C D profile, finer samples for the manual measurement of the basis weight are required.  3.  Better filtering techniques will improve the results for a high resolution profile.  4.  Although moisture variations are probably small, dry weight would provide a better measure for comparison of the results. In this chapter, results of data collection and data processing in the sensor in preparation for the  experimental part of the thesis are presented. Since the dynamic weight and moisture measurements are carried out by a paper machine scanner, it was first necessary to investigate the scanner data processing and calibration scheme. The data processing and calibration scheme have been explained. Backgrounds for processing the signals digitally and the multirate system theory, have been briefly mentioned and the importance of filter design in multirate systems was considered for the analysis and processing of the signals. This investigation was necessary to establish the extent of filtering present when industrial data is used in identification section of this thesis.  59  Chapter 3: Testing and Experimental Methods  Calibrated datafilteredwith a Rectangular window (All 10 scans )  10  15  20 25 30 Data box position (slice position)  Rectangular window used to filter data  0  50 Samples  100 >  35  40  45  Spectrum of Rectangular window  150  0.05 0.1 0.15 Normalized frequency (Nyq=l)  0.2  CD profile with exponentialfilterand Rectangularfilterfor averaging  15  20 25 30 Data box position (slice position)  35  40  Figure 3.24: CD profile calculated based on a Rectangular window averaging scheme.  60  45  Chapter 3: Testing and Experimental Methods  Calibrated data filtered with a Hamming window (All 10 scans )  10  15  20 25 30 Data box position (slice position)  Hamming widow used to filter the data  20  40  60 Sample  80  >  35  40  45  Spectrum of Hamming widow  100 120  0.05 0.1 0.15 Normalized frequency (Nyq=l)  0.2  CD profile with exponential filter and Hamming filter for averaging 2 a bound  <c 2 8  it  1  s 0  -3  l  10  15  20 25 30 Data box position (slice position)  35  40  Figure 3.25: CD profile calculated based on a Hamming window averaging scheme.  61  45  Chapter 3:  Testing and Experimental Methods  Calibrated data filtered with a Harming window (All 10 scans ) 74 |  69 l  31  1  1  1  1  1  1  5  10  15  1  1  1  1  1  1  i  i  20 25 30 Data box position (slice position)  ;  I  35  i 40  CD profile with exponential filter and Hanning filter for averaging 1 1 1 1 r 2 a bound  2 a bound I  5  Figure 3.26:  10  15  I  I  I  20 25 30 Data box position (slice position)  I  L  35  40  CD profile calculated based on a Hanning window averaging scheme.  62  J 45  Chapter 3: Testing and Experimental Methods  Comparison between Measurex profile and a Rectangular window profile with 2 sigma bounds 1  2  1  1  1  1  \ / N x / " '-\\ /  '  V X - A ^ '  >--^^ > 2 cr bound  \ '  ''  v  2-  i  2 a bound  ,-  - / "  i  1  i  i  i  10  15  .  '  X  \  -"z^ ' N ^ — , '  /  V ' X - - '  Measurex profile  "  Rect. window profile  1  5  \ t/ \ — /7T\ \ \  i  i  i  20 25 30 Data box position (slice position)  35  40  45  Comparison between Measurex profile and a Hamming window profile with 2 sigma bounds r~ I i i i i [ r ,^ 2 a bound  10  15  20 25 30 Data box position (slice position)  35  40  45  Comparison between Measurex profile and a Hanning window profile with 2 sigma bounds T 1 1 1 1 r  J  5  10  15  |  |  |  20 25 30 Data box position (slice position)  L  35  40  45  Figure 3.27: Comparison of the Measurex profile with profiles obtained from different averaging methods.  63  Chapter 3:  Testing and Experimental Methods  Comparison of the hand measurement of the sheet and the Rectangular window profile 1  -  1  i  1  i  i  i  v 2 a bound  _____  -^^^^  I  1  1  ~ '',  Hand measurement Rect. window profile  2 a bound i  2  i  i  3  4  i  i  5  i  i  i  6 7 8 Points measured in CD position  9  i  10  11  Comparison of the hand measurement of the sheet and the Hamming window profile 1  1  03-1  i  1  i  i  i  _ 2 a bound  -  1  1  1  -  2 a bound i  2  3  i  4  i  5  i  Hand measurement Hamming window profile"  i  i  i  6 7 8 Points measured in CD position  i  i  9  10  11  Comparison of the hand measurement of the sheet and the Hanning window profile , 2 a bound  Hand measurement Hanning window profile  2 a bound _i i 5  6  10  7  Points measured in CD position  Comparison of the hand measurement of the sheet and the Measurex profile 1  i  I 1  i 2  1  i 3  1  i 4  1  i 5  i 6  1  1  i 7  1  i  i  8  9  Points measured in CD position  Figure 3.28: Comparison of different basis weight profiles with that of the hand measurements of the sample sheet.  64  J 10  I 11  Chapter 4: System Identification  Chapter 4 System Identification  4.1 Introduction This chapter discusses the parameter estimation of the C D response model. The value of the C D response model used here is its ability to link the physical phenomena to a closed form solution. An exact solution is possible only by finite element analysis. The most simple approach uses an empirical model that is fitted to the test data. The model used here combines physically meaningful parameters with a form of solution suitable for identification methods. The standard approaches to system identification for linear systems are not applicable in this case since the weight response is not represented by a linear time series.  The physics of the sheet formation at the wet end  suggests an analogy between wave propagation on the Fourdrinier and transmission and reception of communication and radar/sonar tracking signals. This analogy leads to some techniques available in parameter estimation of signals in radar tracking and application to the parameter estimation of the C D response model on the paper machine.  First, the system identification problem is  addressed and the parameters that need to be estimated are presented. The background theory for the parameter estimation is explained, and the identification approach is posed as a continuous wavelet decomposition problem. Some background on the wavelet transform is explained, and the dry weight model previously developed in chapter 2 is shown to be an admissible wavelet basis function. The estimation procedure is demonstrated with simulated bump tests for light and heavy grades. The estimation method based on wavelet decomposition is verified with identification of model parameters using industrial data. The chapter is concluded with a summary.  4.2 System Identification Approach 4.2.1 The Problem The dispersive wave model, developed as a solution of the Laplace equation, represents the surface of the slurry on the Fourdrinier table. It was shown that at the dry line, the surface represents 65  Chapter 4: System Identification  the shape of the C D response to slice adjustments and the C D response could then be formulated as equation (4.53) :  f}(x,t) - A  e' ^ kd  dry  cos(kx =F U * t )  - h < x < +&i  dTy  The unknown parameters in the above equation are : A , dry  (4.53)  k , k and UJ. The independent variable t^ d  xy  is equal to the time it takes for the sheet to reach the dry line from the beginning of the wire. One of the assumptions in the developed C D response model was that the shape of the initial disturbance at the surface of the slurry is the same as that of the deflected slice lip. This shape is given in equation (4.54), which is derived by setting t=0 in equation (4.53).  •q{x, 0) = A e~  kd  0  The spatial frequency k and damping factor k  d  W cos (kx).  (4.54)  are variables and change as the surface wave  disperses to the sides. This variation is caused by the nature of the drainage on the Fourdrinier. Since the depth of the liquid at the top of the wire changes at each point in machine direction, the speed and frequency of the dispersive wave are affected. It was shown in chapter 2 that the dispersion equation (equation 2.27) defines the relation between the time and spatial frequencies of a dispersive wave. It was also pointed out that the phase velocity or the speed of the wave front is proportional to the square root of the depth. As a surface disturbance approaches the dry line, the depth decreases and the wave slows down. This effect translates to variable spatial frequency k and damping factor k , d  therefore, the value of k and k at the slice are not the same as those at the dry line. In chapter 2 it was d  shown that k is a complex function of depth. The wave amplitude, A , and the damping factor, 0  k, d  also change as the disturbance approaches the dry line; therefore, exact initial values are not critical in the identification process, and the focus will be on the identification of the values of the parameters at the dry line. Correlation between the initial parameters at the slice lip and the parameters at dry line is the subject of the C D control algorithm discussed in chapter 5. In the estimation algorithm it is assumed that the dispersion of the bump response is mainly due to variations in the spatial frequency k and not the phase shift u>t. This assumption applies to typical bump responses on light weight grades such as newsprint and a range of grades on fine paper machines. In these cases an  66  Chapter 4: System Identification  equation such as equation (4.54) is used (phase shift=0). The bump response on heavy grades such as liner board machines can be characterized as an M-shape response. The two maxima in the response are mainly due to the phase shift introduced by the slow speed of the machine which allows more time for the surface wave to disperse (flatten) and for the two maxima to separate (development of phase shift). In these situations an equation similar to (4.53) is used. In summary, for light grade responses with single extremum, parameters to be estimated are: Amplitude A , dry  damping factor k  d  and spatial frequency k. For heavy grade responses with two extrema, in addition to the mentioned parameters, a fourth parameter u is also estimated to account for the phase shift.  4.2.2 Analogy  The parameter estimation scheme to be used was first developed for radar tracking, and, therefore, the analogy between radar tracking and surface wave dispersion is briefly explained.  In target  localization and tracking, a signal is sent towards the target and the reflected signal is received. Comparison of the original and the reflected signal should provide information about the speed and the shape of the target. If the target is moving, the reflection will not only change both the magnitude of the signal but also the phase and the frequency. The change in frequency is normally referred to as the Doppler shift in the case of pure frequency or narrow band signal. When the signal has a wide band, the shift in frequency is best treated as dilation in time [91]. In the communication and radar field the target, and the channel and the medium through which the signal is being transmitted are normally considered together. The model of wave generated as a disturbance at the slice and its dispersion on the Fourdrinier table developed in chapter 2 is analogous to a stationary but dispersive transmission channel in communication theory.  Although there are differences between the two  problems, in both cases, the measured signal is the delayed and dilated version of the original. The two cases are different in that the measurement of the dry weight profile on a paper machine is carried out at the dry end of the machine as opposed to the radar problem where the signal is normally measured at the transmission station. The standard technique for estimating the parameters of the reflected signal in radar signal processing is by cross-correlation. The correlation method is closely related to continuous wavelet theory, which will be explored in the next section.  67  Chapter 4: System Identification  4.3 Wavelet Theory and Correlation Analysis  4.3.1  Introduction  Fourier Analysis consists of the decomposition of the signal of interest in terms of a set of basis function. In any analysis scheme (Fourier, Laguerre, STFT, Wavelet), the shape and properties of the basis functions are characteristic of the signal and better understood giving insight into the properties of the original signal. The basis functions are normally chosen with properties that signify the characteristics of the signal under investigation. In Fourier Analysis the basis functions are pure frequency signals (Sin(x) and Cos(x)) which extend over (-co, + oo) time/space range. Once the basis functions are determined, the next task namely, approximation is considered.  The goal of  approximation theory is to use the least number of basis functions to get a good (in a sense normally defined by an optimum criterion) representation of the signal. Fourier Analysis is an ideal choice when there are long term periodic features in the signal; however, when there are short term periodic components, or a signal with a time varying frequency such as a chirp signal, Fourier Analysis does not localize these frequencies in time. The transform spreads the energy among several frequencies through the construction of a sinusoid with a rectangular windowing function. This effect is illustrated in Figure (4.29). In order to localize these variable frequencies we need to resort to the notion of frequency analysis with time. One approach is to use windowed Fourier Transform or Short-Time Fourier transform (STFT) in which the signal is multiplied with a finite duration window and the Fourier Transform is calculated. This process is repeated by multiples of the window length until the signal is completely covered by the window. The Short Time Fourier Transform maps a onedimensional signal into a two-dimensional space (time and frequency) which results in a square tiling of the time-frequency (t-f) domain. The area of each unit in t-f plane is lower bounded since the time-frequency localization cannot be arbitrarily small. A n increase in time localization (a short window length in time) or resolution will result in a decrease in frequency resolution [32]. In Figure (4.29), two sinusoids are present with 1 Hz and 3 Hz frequencies in adjacent time intervals and a discontinuity at t=l sec. The plot at the top of the figure shows the signal, and the second row shows the Fourier transform and STFT of the signal when mapped into frequency axis. 68  Chapter 4: System Identification  Two Frequencies at adjacent intervals ~r~  0.2  0.4  0.6  0.8  1 1.2 Time in Sees STFT of Combined Signal when Projected into Frequency Axis  Fourier Transform of Combined Signal  5  10 Frequency  15  20  Frequency  Hz  Hz  Scale of cobrs from MIN to MAX Figure 4.29: Distortion in the Fourier analysis for frequencies in adjacent intervals and improvement by Short Time Fourier Transform (STFT). Improvement of Wavelet Transform over STFT by detecting the discontinuity at sample 128 ( t=l sec). 69  Chapter 4: System Identification  In the last row, continuous wavelet transform (CWT) of the signal with Haar bases is shown. The discontinuity is evident in the C W T plot, while it has been distributed over different frequencies with STFT (2nd row). This discontinuity cannot be detected accurately with the use of STFT since the highest frequency is limited by the finite duration window length. To detect the discontinuity in the signal a small window is necessary since the window width establishes the time resolution. This approach, therefore, requires a priori knowledge of the range of frequencies sought in the signal. To overcome the problem of fixed analysis window, some scientists in mid 1980's adjusted the window width to be a function of frequency and developed the Wavelet Transform. Good reviews of wavelet methods may be found in [86], [63], [85], [45] and [77]. Here we present the general ideas in wavelet theory and formulate the C D profile identification problem in the context of wavelet theory.  4.3.2 Definitions and Terminology In Fourier analysis the basis functions are periodic and extend over the entire real time axis. Wavelets are used for representing signals with basis functions that are oscillatory and have amplitudes that quickly decay to zero. Wavelet basis functions are constructed by scaling and translation of a single function known as the mother wavelet.  Scaling corresponds to stretching or compressing  the mother wavelet to generate new basis functions.  Decomposing a signal in terms of wavelet  components is referred to as Wavelet Decomposition or Wavelet Transform.  The decomposition  is carried out by correlating the signal with scaled and translated versions of the mother wavelet. Representation of a signal at small time scales preserves details of the signal and at large scales represents a smoothed view of the signal. A physical analogy of the multiresolution property of wavelet transform as presented in [99] is star gazing with a telescope. The sky in this example is the signal, and the lens can be considered as the mother wavelet. Observing the entire sky, the lens is out of focus or zoomed out. This case represents large scale and gives a global picture of intensity variations of the sky. Zooming in on a particular location of the sky, more detail is observed; however, to see the entire sky at the same zoom level, it is necessary to move (translation process) the telescope across to cover the entire sky. As we zoom in, more details are observed, but the telescope needs to be moved in smaller steps and more often than before. Smaller scale requires more translation to  70  Chapter 4: System Identification  cover the entire signal. Wavelet approximation of a signal consists of the details at different levels, and the degree of accuracy of the approximation is directly related to the scale. The scale, as is apparent in this example, is closely related to the detail level or resolution. Scale can be viewed as a measure of window length of time/space resolution which is inversely proportional to frequency resolution. Scaling and translation when applied to an independent variable (time/space) is called an affine operation. By performing the affine operation on a mother wavelet a set of wavelet bases are constructed, provided the mother wavelet is admissible. The formulation of affine operation is [86]: V>,&(*) = l l " a  ^(~~)  1 / 2  a  •<= Affine Operation  (4.55)  Where a and b are the scaling and the translation parameters and a ^ 0 . In order for the mother wavelet to have desirable properties, typically oscillatory behavior and fast decay, it has to satisfy a set of conditions called the admissibility conditions. If the mother wavelet is ip(t) with Fourier Transform  a sufficient admissibility condition is denned as: co Cp = J  -  2  —  —  duo < oo  (4.56)  — CO  Note that the lower limit of the integration is — co instead of zero. This allows for cases where the mother wavelet is complex and non symmetric about zero. This condition is sufficient and not necessary so that there may be admissible functions that do not satisfy equation (4.56) [70]. If the signal is an energy signal which is normally the case for real signals, then the Fourier transform 4  of the wavelet function will be continuous. For a continuous i^(u>), equation (4.56) can only hold if the following is satisfied. CO  f(0) = 0 #  J  ip(t)dt = 0  (4.57)  — CO  The inner product (correlation) of the signal and the dilated and translated version of the mother wavelet are called the wavelet coefficients. CO  W^f(a,  b) = (f(t), ^ (t)) b  = lap / 1  2  j — CO  4  See appendix for definition of energy signals.  71  f(t)  $  dt  (4.58)  Chapter 4: System Identification  Equation (4.58) is the formulation of wavelet transform.  The notation (•) refers to inner/sealer  product operation. The wavelet coefficients are unique for a given mother wavelet. This means that, contrary to Fourier transform, the wavelet transform must always be associated with a mother wavelet and expressions such as " wavelet transform of a signal" must be identified with a specific mother wavelet. The inverse transform to reconstruct f(t) from the wavelet coefficients is a double integral since wavelet transform is a mapping from time/space domain to a scale-translation (shift) domain. CO  f(t) = C^  f  CO  j.W^f(a,b)i> , (t) ^ 1  a b  (4.59)  — OO — C O  The parameter C^  1  is the inverse of the admissibility function. The choice of mother wavelet depends  on the application. The choice should be based on the characteristics of the wavelet function and the features of interest in the signal. A common choice of wavelet in detecting sharp edges and discontinuities is the Haar wavelet because of its step-like shape. In the application of wavelets to system identification of the dry weight response shape on paper machines, we resort to the bump response as the basis function and show that the equation derived in chapter 2 for the surface wave satisfies the admissibility conditions. In the case of Fourier Analysis the basis functions are well localized in frequency domain. In the continuous case the change in frequency is continuous and any frequency localization can be achieved as accurately as desired. On the other hand, in time domain, Fourier bases extend from (-co, + oo), hence no time localization exists. A perfect localization in time domain can be achieved by the impulse function with poor localization in frequency domain (frequency band covers the whole spectrum). One advantage of the wavelet basis is to achieve a good time and frequency localization at the same time. By increasing the time domain window length we allow low frequencies to be captured by the window and effectively reduce the frequency range that could be distinguished with the window . In other words, as the time interval (window length) increases, the frequency interval decreases and vice versa. This suggests that there is a limit to the localization of time and frequency at the same time. Gabor [32] has shown that if a Gaussian function is used as the window function, the uncertainty in time and frequency will be lower bounded by \ ( Af At > i). The choice of other functions makes the product of the uncertainties in time and frequency larger than \ but reasonably 72  Chapter 4: System Identification  close to this limit. The minimum bounds on time and frequency in time-frequency domain constructs a rectangle referred to as the information cell. Wavelet Classification: There are different classes of wavelets, each suitable for a specific application.  The common features among them are the existence of a mother wavelet, and the  dilation and translation operation (affine operation). If the affine operation domain is considered as the classification criterion, the wavelet transform can be divided into continuous and discrete transforms. In the continuous case a and b, the scaling and translation parameters, and the time parameter t are changed continuously. In the discrete wavelet transform all the parameters are discretized as  « = <c => {fl>m,n{t)}  (4.60)  b = nb a™ 0  and the dilation equation similar to equation (4.55) is given by  ^ , n ( < ) = a~ l  m 2  ^{a~ t  - nb )  m  m,n£  0  Z  (4.61)  The discrete wavelet transform relation can now be defined as [DWT(f(t))}  mn  = W^f(m,n)  =  (/(*), Vw(0>  CO  - "~  m/2  J  fit)  i){a~  m  t - nb ) dt  ( 4 - 6 2 )  0  The discrete wavelet transform can further be divided into: non-orthogonal, semi-orthogonal and orthogonal wavelets. If the wavelets are redundant, they are not linearly independent. This redundancy can be advantageous in some applications since the overlap of information in wavelet coefficients reduces the accuracy requirement on the storage or transmission of them over a communication channel.  Semi-orthogonal bases are orthogonal between scales, but the orthogonality condition is not  satisfied within a fixed scale. Eliminating the redundancy results in an orthogonal basis where the wavelets are linearly independent. Orthogonal wavelets are important in the application of wavelet theory in signal processing and image analysis. Closely related to the orthogonality of the bases is the concept of frames. Frames are sets of discrete wavelets that satisfy the following condition.  Mmw < E i < / ( * ) . ^ . n ( * ) ) i 2  73  2  < B\\f(t)\\  2  (4.63)  Chapter 4: System Identification  highpass  Details Scale 1  Input  l o w p a s s  Details 2  highpass —I rT\  Scale S i  lowpass  Apprx. Scale 2  Apprx. Scale 1  Figure 4.30: Multiresolution wavelet decomposition scheme. A and B are the frame bounds and B/A is their ratio . Both A and B are positive numbers. The ratio 5  of frame bounds is a measure of redundancy in the wavelet bases. If the ratio is close to one the frame is called a tight frame and reconstruction of the original signal from wavelet coefficients can be achieved with good approximation. A large frame bound ratio represents a loose frame that offers robust reconstruction at the presence of noise [70]. The multiresolution scheme developed by Mallat [60] systematically constructs orthogonal basis. The generation of wavelet basis in multiresolution scheme involves construction of a lowpass filter called scaling function. The lowpass and highpass filters derived from scaling function and its highpass equivalent, the wavelet function, are the main operators for the generation of approximations and details of the signal at different resolutions. A schematic of the multiresolution scheme is shown in Figure (4.30). A more general approach to orthogonal decomposition is the Wavelet Packet Decomposition [17]. In multiresolution analysis the lowpass branch is decomposed and the details are retained. In wavelet packets the detail branch as well as the approximation branch may go through the decomposition process and every branch can be divided into details and approximation at a lower level. decomposition generates an orthonormal bases as a binary tree.  This  Coifman and Wickerhauser [17]  developed a fast algorithm to search for the best basis. The optimization criterion is based on the least information cost. A schematic of the wavelet packet algorithm is shown in Figure (4.31).  4.3.3 Continuous Wavelet Transform (CWT) The main analysis and synthesis equations for the continuous wavelet transform (CWT) have See appendix for the definition of norm  74  Chapter 4: System Identification  highpass highpass  Details  ) )  H  1  H  33>  L  ft  >  Input  * low-pass highpass H  lowpass  Apprx. L1  L  ft ft  Details H1H2  Apprx. H1L2  Details L1H2  Apprx. L1L2  lowpass  Figure 4.31: Wavelet packet algorithm, been presented in equations (4.58) and (4.59). The mathematics of the C W T is very similar to that of short-time Fourier transform (STFT). The main difference is in the transform domain where the STFT is always a complex function regardless of the type of input signal; but in the case of C W T the coefficients are real for real valued signal and the mother wavelet [23].  Continuous wavelets  are non-orthogonal. This means that unlike discrete wavelet transform there is no scaling function associated with the mother wavelet and perfect reconstruction is not possible. The main features that distinguish continuous wavelet transform from its discrete counterpart are the scaling and translation parameters. In the continuous case, a, b and t can assume any real value, while in the discrete domain the scale normally changes by factors of 2 and the shift/translation is a function of scale as shown in equations (4.60) and (4.61). In image analysis and multiresolution applications, the discretization by a factor of 2 is appropriate but in some applications a small fraction of change in scale may have a significant impact on the approximation. The numerical implementation of C W T on computers requires discretization of transform domain parameters. The scale parameter a can vary between the starting value and any desired level based on the accuracy requirement of the problem. The translation operation also moves the bases smoothly over the signal and is completely independent of the scale parameter. The coefficients of C W T show the similarity of the portion of the input signal and the basis function at the specific scale and shift parameter values. If the degree of similarity is high the coefficient is large, while in the case of dissimilar features the coefficient is small. This is the underlying idea of the identification of the C D response model.  75  Chapter 4: System Identification  4.3.4  Correlation Analysis  Correlation analysis is an important tool in identification and time-series modelling of the systems. In this section Wide-band Correlation Processing of the signals is discussed and shown to be similar to continuous wavelet transform. Based on the analogy drawn between radar detection of a moving target and a dispersive channel, the continuous wavelet transform is shown to be a useful tool for the identification of the parameters in the C D response model. Wide-band correlation processing is used in the identification of signals with fractional bandwidths, non-stationary signals, signals with large bandwidth, and signals reflected from a fast moving object [96]. The analogy presented here is between a dispersive channel and a moving target with definitions and terminology often referring to radar processing with a moving target by sending a known signal that is compared with the reflected signal. This comparison is carried out by correlating the input signal with its reflected counterpart. Knowing the signal transmission speed one can locate the object. If the target is moving with a constant speed, the frequency of the reflected signal is shifted by an amount known as the Doppler shift. The Doppler shift explains the frequency shift of the reflected signal if the incoming signal is narrowband. More generally, the term wide-band refers to a signal that has a wide range of frequencies that are not shifted equally when reflected off a moving target. A narrow-band signal contains frequencies that are very close in the frequency band so that when frequencies are shifted due to reflection off a moving target, all the shifts can be considered with respect to a single reference frequency (fundamental carrier frequency). In the case of a wide-band signal, since the original constituent frequencies are not close, the frequency shift after reflection off the target cannot be approximated as offset from a single reference frequency. More precisely defined, a signal is considered wide-band if the spectrum of the signal cannot be separated into two components, one being the centre or carrier frequency and the second an offset from the carrier frequency. In cases where the reflected signal undergoes dilation rather than frequency shifting the correlation process is carried out between the received signal and the dilated and delayed versions of the transmitted signal. This technique is referred to as matched filter processing and is equivalent to wide-band correlation processing which is the same process as in the continuous 76  Chapter 4: System Identification  wavelet transform. For matched filter processing and the continuous wavelet transform to correspond, the signal used as the mother wavelet must satisfy the admissibility condition.  4.3.5 Proposed System Identification Approach The identification process is based on the continuous wavelet transform. The Fourdrinier table at the wet end of the paper machine is considered as the transmission channel for the disturbance generated by slice lip deflection. A dispersive transmission channel acts in a similar manner to the Fourdrinier table and scales and translates the transmitted signal before it settles at the dry line to be finally measured by the scanner. It is clear from the dispersive nature of the Fourdrinier table that in order to estimate the parameters of the surface wave at the dry line one has to correlate the bump response with dilated and translated versions of the transmitted signal. This process can be viewed as a matched filtering scheme or a continuous wavelet transform of the bump response with respect to the mother wavelet obtained from surface wave equation. In order for the equation (4.53) to be a mother wavelet it has to satisfy the following admissibility condition. For the bump response to be an admissible function, equation (4.53) should satisfy condition (4.56). Condition (4.56) requires the Fourier transform of the mother wavelet. The Fourier transform of equation (4.53) is taking (A  = 1):  0  — CO  e~  ax  cos (kx) e~  3WX  dx  0 < x < +oo  o  +  0  - co < x < 0 — CO  Using standard Fourier transform tables: a + ju! (a + jwf  0 < x < +co  + k  2  + a- ju> (a — jw)  2  — oo < x < 0  + k  2  77  (4.64)  Chapter 4: System Identification  The admissibility function  can now be calculated by inserting in equation (4.56) which gives:  0  (4.66)  — OO  Admissibility condition requires that relations in equation (4.66) are upper bounded. By carrying out the integration and calculating the limits it can be shown that C.^ has a finite value . One of the 6  consequences of the admissibility condition is that the mother wavelet has a zero mean. This result is in agreement with the nature of the C D variations in general and bump response as a special case of C D variations. The mean value of the C D variations of the sheet (moisture, dry weight, basis weight, etc) is always attributed to machine direction and C D variations have zero mean. The bump response used as the mother wavelet, as a special case of C D variation due to a single actuator move will, therefore, have zero mean. Having shown that the C D response model is an admissible function, we now present the application of continuous wavelet transform (CWT) with the C D response model as the mother wavelet to identify model parameters.  It should be stressed that although the response model is  admissible it may not be a good candidate for other applications, especially for discrete wavelet transform (DWT). Regularity is one of the main features required in D W T applications. Regularity is directly related to the differentiability of the mother wavelet at x /1 = 0. The C D response model can be made to have a maximum at x-0  by introduction of a phase shift in equation (4.53), but  the second derivative of response is not zero. The regularity of the mother wavelet is related to the degree of accuracy of the bases to approximate the details at high frequencies (small scales). If the mother wavelet is not differentiable many times at zero, a point wise convergence will not occur at high frequency approximations. This issue does not surface in the application of C W T to dry weight response identification since a bump response is the highest bandwidth that is controllable with C D actuators and smaller scales are treated as uncontrollable variations and not approximated. This will  6  Details of the integration procedure and proof are presented in the Appendix.  78  Chapter 4: System Identification  be discussed when the controllability of the C D variations and bounds on the continuous scale for identification are discussed in chapter 5. A schematic of the identification procedure is shown in Figure (4.32) which shows the correlation of the bump response with the mother wavelet and its dilated versions. The selection criterion is maximum correlation. If there is more than one bump present in the data, then all the extrema of wavelet coefficients will be considered. In the diagram, the first section is the continuous wavelet transform of the bump response with the generic shape of the bump response taken as the mother wavelet (translation of C W T bases is embedded in the correlation process). The dilation property of the Fourdrinier table can be interpreted in two different ways. One view is, as the disturbance propagates towards the dry line, the number of samples stays constant and the sampling distance (sampling period) expands or dilates. This view treats the problem as a multirate system where the existing tools can be utilized to change the sampling rate. The other view is to keep the sampling distance constant and dilate the shape or change the parameters of the model.  In the parameter  estimation we resort to the second interpretation. The relation between the mother wavelet parameters and the estimated parameters is described in equation (4.67), where MW wavelet and wavelet basis respectively.  and WB stand for mother  Subscripts / and j refer to scale and shift in the wavelet  domain. The shift index j is between 0 and the number of samples of the basis function. The scale index / shows the number of divisions in scale parameter that varies between 1 and desired accuracy of the identification scheme. The trade off between the accuracy of identification with C W T and the resolution of scale index i is discussed later in the chapter. When the response is that of a light grade with a single extremum, the shift index j is zero. In the identification of a bi-modal response, a non zero translation index indicates a phase shift in the mother wavelet. In the bi-modal case, the extra parameter u can be estimated through the estimation of u>t\ according to the shift index j in wavelet domain. If the location of the dry line and machine speed are known MW  =WB {x) QJQ  WBijfa)  where  = e  _ o | a ; |  = e~  alxi1  i = scale  cos (kx)  cos (k ) Xl  - &_ < x < +b  1  - b < x < +b 2  j = shift/translation  79  and UJ can be determined.  2  (4.67)  Chapter 4: System Identification  At the dry line parameter x is scaled to x\ and the following linear relations hold between parameters at the slice lip ( k, a) and those at the dry line ( k], a/). taxi <  = a\x  ,  and xi =  ,  scale  a scale  •a  x  k scale For the bi-modal case, the relation between parameters at the slice lip and the dry line are given: MW  =WB (x,0)  = e"  Ofi  WB^jix-utx) where  =e  i — scale  cos (kx)  a|3:|  - a | x i 1  - h < x < +b  1  c o s ( f c i i T w i i ) - b < x < +b 2  &  2  (4.69)  j — shift/translation  The following linear relations hold between parameters at the dry line (k/, ai) and those at the slice lip (k, a). ax\ = a\x kx\ = k\x  <^ after u\,\ phase shift  = (j * sampling distance) ^  Dry line Distance j  Machine Speed a (ai = scale * scale  ^  X  l  ~ scale  4.3.6 Alignment Alignment in the context of the C D variation analysis refers to the adjustment of the position of the actuator response with respect to its expected location. When an actuator is moved, the expected response at the dry line should be such that the centre of the actuator and its response coincide. If the centres do not coincide, there is misalignment between the actuator and the C D profile. The normal method of determining alignment is through a bump test. 80  Chapter 4: System Identification  The misalignment can be caused by the movement of the sheet and also by shrinkage as the sheet dries in the dryer section of the paper machine. Measurement of distance of the actuator centres from the centre of the response locations is a measure of the misalignment of the actuators with the profile. The C W T can be employed to find out whether the actuators are aligned. The misalignment can be estimated at the same time as response parameter estimation is carried out. The extra effort required is to map the centres of the bump locations from actuator resolution to the C D response resolution which can be easily carried out with the multirate filtering scheme introduced in chapter 3. The first step is to change the sampling rate of the bump locations from actuator resolution to databox resolution. The C W T is performed as in the model parameter estimation scheme and the distance between the nearest bump location and the centre of the identified bump indicates the amount of misalignment for that particular actuator. A schematic of the alignment algorithm is shown in Figure (4.33).  CD bump shape g(x)  f(k1 x) CD bump L shape g(x)  f (k2 x)  <x> CD bump , shape g.x)  f (k3 x)  Selection Criterion  <x>-  Estimated Parameters  (Max CWT coefficient)  f (kn-1 x)  CD bump shape g(x)  CD bump shape g(x)  f (kn x)  Dilated Replicas of bump shape  Correlation / Translation^. Selection of best Process \ in CWT / matching parameters  Identification with Continuous Wavelet Transform  Figure 4.32: A schematic of the identification procedure with Continuous Wavelet Transform.  81  Chapter 4:  11 Actuator resolution  V  Multirate scheme sample rate change  1  Databox resolution Measure distance between centre of response and nearest actuator centre  CD response with bumps  Continuous Wavelet Transform  System Identification  Find Max I wavelet coef |  Misalignment Measurement  Selected coef. and shift indices  Figure 4.33: Determination of misalignment using CWT. Damping Parameter  Mother Wavelet  Mother Wavelet  Spatial Frequency k amplitude  Equation (4.54)  0.1  0.2  1.0  Table 4.5 Mother wavelet parameters used in Figure (4.36) 4.4 Parameter Estimation 4.4.1  Estimation with Simulated Data  In order to illustrate the identification approach, a set of simulated bump tests was generated with the help of the model developed in chapter 2. To make the simulation more realistic, a white noise sequence with zero mean and a variance of 0.1 was added to the response.  The simulated  responses were carried out in two different cases, a single mode response which characterizes a response on a light grade sheet and a bi-modal response that best describes bump responses on liner board machines. Figure (4.36) shows the single-mode bump representative of light grade sheet. The response was generated by running the simulation of the model for a slice opening of 1.5 cm and a headbox consistency of 0.4%.  A white noise sequence of zero mean and a variance of 0.1 was  added. The plot at the top shows the simulated bump response. The mother wavelet was chosen as in equation (4.54) with initial conditions shown in Table (4.5).  The middle plot shows the wavelet  coefficient in the wavelet transform domain. The x-axis indicates the zero shift/translation in C D position (samples) and the y-axis shows the scale parameter. The scale was chosen to be linear between 1 which represents the mother wavelet resolution and 20 with 100 divisions. The maximum coefficient occurred at sample 50 (zero translation) and scale level 3. The estimated parameters are  82  Chapter 4: System Identification  Estimated  Amplitude (Wavelet  Parameters at Dry line  Damping Factor (kdi) Spatial Frequency (k\) 0.07  0.14  Coef, A i ) 0.194  Table 4.6 Estimated parameters for the example in Figure (4.36). calculated based on relations derived in equation (4.67).  The last plot shows the estimated response  which is simply plotting the basis function at scale level 3 and the proper gain obtained from C W T of the response shape. The average two-sigma (2*standard deviation) of the estimation error was found to be 0.043. A plot of normalized autocorrelation of the estimation error with 95% confidence interval is shown in Figure (4.39). This plot is an indication of the goodness of fit. The plot indicates a high correlation at zero shift and small correlation for other shifts. This shows that the estimation error has similar characteristics to a white noise sequence. A summary of the estimated parameters for the model is shown in Table (4.6). Figure (4.37) shows the results of estimation algorithm applied to a simulated bi-modal response. This example was generated with the model in chapter 2 when the slice opening was set to 2.5 cm and a headbox consistency of 0.9% was used. A white noise sequence of zero mean and a variance of 0.1 was again added to the C D response. A careful analysis is required in the application of C W T algorithm to bi-modal responses. A bi-modal response is actually two single mode waves that are moving apart with a certain phase velocity. If the C W T is blindly applied to a bi-modal response the highest correlation is normally obtained with a single mode response at zero shift. In fact C W T tries to approximate the overall shape of the response with a wide single-mode basis. The correct application of C W T in this case is to consider only half of the response and treat it as the input to C W T algorithm. This approach isolates each mode and separates them into two individual modes. When the new input is processed with CWT, maximum correlation peaks at the correct translation index and scale. The first plot in Figure (4.37) shows the estimated shape with the solid line and the actual response with the dotted lines. The estimated response is in good agreement with the simulated response. The middle plot indicates the estimation error at different C D positions. The last plot in the bottom row shows the normalized autocorrelation of estimation error with 95% confidence interval. The graph indicates that the autocorrelation is within the 95% confidence band. The average  83  Chapter 4: System Identification  Damping Parameter  Mother Wavelet  Mother Wavelet  Spatial Frequency k amplitude  k  d  Equation (4.54)  0.08  0.22  1.0  Table 4.7 Mother wavelet parameters used in Figure (4.37) Translation Index Amplitude (Wavelet (phase shift in  Damping Factor (kai)  Spatial Frequency (k\)  0.05  0.14  Coef, Ai)  databox unit) 2  1.481  Table 4.8 Estimated parameters for the example in Figure (4.37). two-sigma of the estimation error was found to be 0.190. The estimated parameters are presented in Table (4.8) with the mother wavelet parameters shown in Table (4.7).  4.4.2 Estimation with Industrial Data In this section, the identification procedure with the continuous wavelet transform is applied to two sets of industrial data. The first set consists of three bumps across the lip on a machine manufacturing heavy grade paper. The plots in Figure (4.34) were produced by averaging the data over 10 scans before the bumps and subtracted from an average of 10 scans after the bumps had been performed. This method isolates the effect of the bumps from steady state variations. In the identification exercise the last bump in Figure (4.34) was isolated and considered as the input to the C W T estimation algorithm. The bumps were performed on a 81 lb/ream (3300 sqft) grade. The signal to noise ratio is small although the exact amount is not known. For the identification of the response shape the same procedure as in the simulated case is used. The initial parameters (mother wavelet parameters) are shown in the first row of Table (4.9). The estimation results are presented in Table (4.10) and Figure (4.38). The plot at the top of Figure (4.38) shows the bump and the actuator set point (actuator set point is not to scale). The two-dimensional image in the middle is the mapped wavelet coefficient in wavelet transform domain. The best fit happens at the scale level 25 and zero shift (sample 72). The two-sigma value of the estimation error was found to be 0.665. The twosigma value compared to that of the simulated case is an order of magnitude larger which indicates 84  Chapter 4: System Identification  the high noise content of the input signal. The normalized autocorrelation of the estimation error with 95% confidence interval is shown in Figure (4.39). The autocorrelation shows a dominating spatial frequency of 1/4 (Cycle/databox width). This frequency can be visually observed in Figure (4.38). The second set of industrial data was taken from [48]. Figure (4.35) shows the bi-modal response which is typical of liner board machines. The same treatment as in the case of simulated bi-modal response is given in the case of industrial data. The initial parameters used in the mother wavelet are given in the second row of Table (4.9). The identification results are presented in Figure (4.41). The top plot shows the overlay of the estimated response and the input to the identification algorithm. The middle plot indicates the estimation error. The maximum estimation error occurs at the negative lobes of the response. Although the overall shape of the response has been estimated to a reasonable level the model mismatch is obvious in the negative lobe locations. This sort of error is expected since the model has a single parameter for damping which only varies with each resolution (scale) while the input signal represents a sinusoid with a variable damping factor in C D position. The last plot shows the normalized autocorrelation of the estimation error with 95% confidence interval. Cross direction positions at which model mismatch is more significant fall outside the confidence limit. The overall two-sigma of the estimation error was found to be 1.13. The estimated parameters are listed in Table (4.11).  4.4.3  Model Validation In general, there are two approaches to model validation. These methods are; the use of plots  and visual comparison of the estimated model response with the measured data and statistical tests on the estimation error [84]. In the first approach a visual inspection of the overlays of the model and the measured output indicates the goodness of the fit. The deviation of the model and measured output come from the modelling errors and the disturbance. Although there are several statistical tests available for the second approach one standard technique has been used for this analysis. The main focus on the statistical approach is on the estimation error or residual. The discussed validation methods are based on the following assumptions [84]: 1.  e(t) is zero mean white noise sequence. 85  Chapter 4: System Identification  Bump test on 3 actuators across the lip  0.5  M  200 Databox 40  Differential actuator set point showing bump locations  20  'a.3 8 o  0  3  < -20  -40  10  15  20 CD Position  25  30  35  Figure 4.34: Results of a bump test at three different zones across the lip (Courtesy of Measurex Devron Inc).  Normalized bump response on a liner board machine taken from Karlsson.et.al  80  100 CD Position  120  140  160  180  Figure 4.35: Bump test on a liner board machine from reference [48] showing a typical M-shape response.  86  Chapter 4:  System Identification  Mother Mother  Damping  Figure #  Spatial Frequency k Wavelet  Wavelet  Parameter amplitude  Figure (4.38)  Equation (4.54)  0.09  0.2  1.0  Figure (4.41)  Equation (4.54)  0.15  0.4  1.0  Table 4.9 Mother wavelet parameters used in Figures (4.38 and 4.41). Estimated  Spatial Frequency  Parameters at Dry  Damping Factor (kdi) 0.016  Line  Amplitude (Wavelet  (Ai)  Coef, i )  0.035  1.96  x  Table 4.10 Estimated parameters for the example in Figure (4.38). Translation Index (# of databox) 7  Amplitude (Wavelet Damping Factor (kdi) Spatial Frequency (k\) 0.07  Coef A i ) 5.36  0.18  Table 4.11 Estimated parameters for the example in Figure (4.41). 2.  e(t) has a symmetric distribution.  3.  e(t) is independent of the past inputs.  4.  e(t) is independent of all inputs if the process is in open loop.  The two main statistical tests are autocorrelation and cross-correlation tests.  Based on the first  assumption, if eft) is a white noise sequence then the autocorrelation is zero everywhere except at r = 0.  When we are estimating the autocorrelation function from measured samples, the  autocorrelation function estimates will not be zero unless the number of samples tends to infinity. The autocorrelation estimate of a white noise sequence is expected to be large at the first sample and small at other samples (shifts). The bounds which determine what values of the autocorrelation should be considered small are found using the confidence interval. A common choice of confidence interval range is between 90% to 99%. This means that all the values of the autocorrelation that fall within the confidence band around the mean are to be considered small enough to satisfy white noise criterion and the risk associated with this decision is only 5% (confidence interval=95%). 87  Chapter 4: System Identification  The second method among the statistical techniques of model validation is the cross-correlation technique. This method is based on the first, third and fourth assumptions in the list. The underlying idea is that if the input to the system is correlated with residuals, useful information still exists in the residuals and the estimation is not accurate. In the case of good fit to the data, the residuals are white and the cross-correlation of the input with the residuals will be zero. This explains the case for positive shifts (r > 0). When considering the negative shifts interesting results can be deduced. If the input is white, cross-correlation of the input and the residuals for negative shifts will be zero regardless of the model accuracy. A nonzero cross-correlation for negative shifts points towards the existence of feedback during the identification process. The other use of this test is to check for the presence of feedback while the identification experiment is being carried out. Clearly this method relies on the information of the input. In the parameter estimation of the C D response the general form of the input was used as the mother wavelet and the specifics of the input was not considered. The correlation method explained earlier is a better choice for model validation of the C D response model and was performed on the simulated cases and the industrial data set.  Figures (4.39) and (4.40)  show the normalized autocorrelation of the residuals with 95% confidence interval. Figure (4.39) resembles characteristics of a white noise although few points cross the confidence band. Figure (4.40) appears to indicate a less accurate model, but further comparison of the estimated model and the original measurement show that the frequency that appears in the autocorrelation of the residuals is also present in the measurement. This points towards the presence of a disturbance rather than an inaccurate model. If a notch filter is used to filter out the cyclic disturbance, the autocorrelation of the residuals after the first sample remain in the confidence band. The comparison of the two identification tests on the bi-modal responses for the simulated and industrial data show that in the case of the industrial data the modelling is less accurate. The number of autocorrelation points that lie outside of the confidence band are larger than that in the simulated case.  4.4.4 Scale Bound Considerations for CWT  There are two issues regarding the discretization of the continuous wavelet domain parameters, the limits of spatial frequency and damping coefficient of the response model and the other the 88  Chapter 4: System Identification  discretization level in the scale parameter. The translation parameter is normally decided on the basis of sampling time/distance of the input signal since the highest accuracy in the translation parameter is achieved by shifting the wavelet bases for every sample. The limits of the scale parameter are closely related to the controllability of the system. If the scale parameter is less than one, the mother wavelet shrinks and as the scale (>1) increases, the mother wavelet dilates. In the context of the application in C D response parameter identification although the whole range for the scale parameter is theoretically valid, the physics of the problem can limit the scale. This limitation is beneficial since it reduces the computation cost. The limit is introduced by the maximum bandwidth that can be controlled by the C D actuator. Since the focus of this research has been on the C D response model of the weight profile on the Fourdinier machine, the discussion concentrates mainly on response to the C D weight actuators; however, the concept could easily be applied to other C D response models. In the case of C D weight profile, the main actuator is the slice lip, and deformation of the slice lip generates variation in the C D profile. The highest frequency or, in the wavelet terminology, the smallest scale that is attainable with the slice lip is limited to the response width generated by movement of a single actuator. From a controller point of view, this is the highest controllable frequency which limits the controller bandwidth. In the identification experiment, the scale that is close to the response width should be considered the lower bound for the scale parameter. This simply means that although the wavelet analysis could distinguish smaller scales in C D variations, since they are not controllable due to limitation of the slice spacing, they should be treated as noise. The upper bound for the scale depends on the application. If the purpose of the identification experiment is to find model parameters for a bump shape, there is no need to extend the upper bound of the scale parameter to cover the D C level. In the case of C D control, this upper bound should be larger compared to response identification case, since the objective of the control is to identify a wide range of frequencies and specially the range that covers the low end of the frequency spectrum. The quantization level of the scale parameter or the resolution, is closely related to the accuracy of the estimation algorithm. If the measurement is assumed to be noise free, and the divisions in the scale parameter are represented by Aa, then in the worst case, the estimated scale can be in error by ^ r . There is always a trade off between the quantization level of the scale and the computation cost. 89  Chapter 4: System Identification  4.4.5 Choice of Mother Wavelet Parameters The critical mother wavelet parameters in equation  (4.53) are the damping coefficient and the  spatial frequency. As a general approach, this problem can be formulated as an optimization problem on a two-dimensional space. This is a black box approach to parameter estimation, and does not make use of a priori information about the initial state of the system. Making use of the available information about initial parameters is especially attractive since identification with C W T of model parameters showed little sensitivity to the choice of initial parameters, if they are within range. The initial parameters of the bump shape at the slice lip are good choices of mother wavelet parameters.  4.5 Summary In this chapter, the identification of the model parameters was considered. Since the parameter estimation cannot easily be formulated as a time series, a novel identification algorithm borrowed from the radar processing field was introduced. The underlying algorithm used is the continuous wavelet transform. A brief introduction to the theory of wavelet transform was first presented, and the model developed in chapter 2 was shown to be a candidate for the mother wavelet function. A continuous wavelet transform with respect to this new basis was applied to identify single mode and bi-modal responses with simulated data. The application of the estimation scheme to two different industrial data sets gave good results.  Autocorrelation techniques were employed to validate the  estimation and modelling scheme. Finally, some issues regarding the choice of initial parameters for the mother wavelet and limits on the scale parameters to reduce the computation cost were discussed.  90  Chapter 4: System Identification  Simulated bump response with output noise  50 CD Position  60  Continuous Wavelet Transform of the simulated bump •  90 xo u 70 •a  a 60  1 J  § 40 o U  30 20 10 10  20  30  40  50 Shifts  60  70  80  90  100  Response identification using Continuous Wavelet Transform 0.25  1  0.2  T  1  r  i  H  0 1 0.05  i  \ \  // //  Estimated  \ V \\  //  /' /j  -  tl  £  i  - — Actual  0.15  1  i  i  Li  *  0  /  \  /\  i —  , —  -0.05 -0.1 -0.15 0  10  i  i  20  30  i 40  1  50 CD Position  1  1  1  1  60  70  80  90  100  Figure 4.36: Identification of a simulated response on a typical light grade with Continuous Wavelet Transform.  91  Chapter 4:  Identification of a simulated heavy grade (m_shape) response with CWT  40 Databox Estimation error for the identified m_shape response  10  20  30  40 50 CD Position (databox)  Normalized autocorrelation of estimation error with 95% confidence interval  Figure 4.37: Identification of a simulated response on a liner board machine with Continuous Wavelet Transform.  92  System Identification  Chapter 4: System Identification  Bump Response on an 81 lb/ream grade  60  80 CD Position  Continuous Wavelet Transform of a bump response from a 81 lb/ream grade  10  -10 60  80 Shifts  Response shape identification of industrial data with Continuous Wavelet Transform  60  80 CD Position  Figure 4.38: Identification of a bump test response on a 81 lb/3300 sqft grade with Continuous Wavelet Transform.  93  Chapter 4: System Identification  Normalized autocorrelation estimation error with 95% confidence interval  ' 0  10  20  30  40  50 Shifts  60  70  80  90  100  Figure 4.39: Normalized autocorrelation of estimation error with 95% confidence interval for the example shown in Figure (4.36).  T  1  Normalized autocorrelation of estimation error 1 1 i  r  c  J  I  I  L  20  40  60  80  100  120  Shifts  Figure 4.40: Normalized autocorrelation of the estimation error and 95% confidence interval for the example shown in Figure (4.38).  94  140  Chapter 4: System Identification  O  i  Identification of a bump response on a board machine with CWT i i i i i i i  2  -2 -4  / \/ \  -  0 L.  /  c ^ ^ x  — "\  V  /  Estimated -  \  /^x^v  \  -  -  0  20  40  1-51  1  1  i 20  i 40  _i  Actual  A A  4-  i  51  0  60  80  100  120  140  160  i 180  200  i 180  I 200  Estimation error for the identified m shape normalized bump response 1 1 1 1 1 r  I  60  i 80  i 100 CD Position  i 120  i 140  i 160  Normalized autocorrelation of estimation error with 95% confidence interval  Shifts  Figure 4.41: Identification of a bump test response on a liner board machine (Industrial data taken from [48]) with Continuous Wavelet Transform.  95  Chapter 5: Cross Directional Profile Control  Chapter 5 Cross Directional Profile Control  5.1 Introduction This chapter addresses the C D control problem using results obtained in chapters 2 and 4 on modelling and identification with CWT. The underlying theme of the control algorithm is to combine the wavelet decomposition of the C D profile, the shape of the slice at the headbox and the resulting simple correlation between these shapes (profiles at slice and scanner). The C D control problem is considered as a decoupling algorithm that is only concerned with the shape of the C D profile. A l l dynamic effects, including the time delay between actuator movement and the resulting variation at the scanner, are attributed to machine direction.  The traditional approach to profile  control is to use an interaction matrix to map the C D variations from scanner to actuator resolution. The mapped/decoupled profile is then subtracted from the desired profile and gives individual set points for each actuator. The decoupled set points are fed to a series of single-input single-output (SISO) controllers associated with each actuator. C D control is thus a multivariable problem which reduces to a series of single-input single-output systems after decoupling. In this chapter a new C D control/decoupling algorithm based on wavelet theory is presented, which allows the mapping problem to be viewed either as a change in the scale parameter or as a compression scheme. First, the possible approaches to the C D control problem using wavelet decomposition theory are considered, and then the approach based on developments of chapter 2 and 4 is presented. A n overview of the chosen control algorithm is given and the interpretation of profile mapping as a change in the scale parameter is described. The M D dynamics are modelled as first-order lag plus time delay, and the standard Dahlin algorithm is explained. Finally, the performance of the proposed C D control algorithm is evaluated by application of the algorithm to the bump test results used in chapter 4 for identification of the response shape. Although a true test of the algorithm is not possible in a simulated environment, the simulation can provide guidelines for potential applications in the field. This chapter is concluded with a summary. 96  Chapter 5: Cross Directional Profile Control  5.2 C D Profile Control with Wavelets In this section, two different approaches to C D control with wavelets are explained. The first method uses orthogonal polynomials to represent variations in the measured C D and actuator profiles. Since the number of actuators is normally smaller than the measurement points, the static input/output relation can be described with a non-square band-diagonal matrix of gains normally referred to as the interaction matrix. The decomposition of the input/output profiles with a series of orthogonal polynomials has been addressed by Grimble [37] for sheetmaking processes in steel industry, and by Kristinsson [54] and Rigopoulos et al [76] in papermaking applications. The underlying idea in this approach is to decompose the input profile in terms of an orthogonal basis (Gram polynomials in [54]) and relate it to the decomposed C D profile through the interaction matrix. The decomposition can be thought of as a series representation of the input/output profiles in terms of an orthogonal basis with the corresponding weighting coefficients.  The C D control problem then reduces to the adjustment  of input coefficients in order to achieve the desired set points expanded in terms of the orthogonal basis and the weighting coefficients. C D control can also be formulated as an optimization problem where the criterion index is a quadratic function of the deviations from the desired profile. Applying the same strategy using wavelets as the basis functions results in a new C D control algorithm. The orthogonality requirement suggests the discrete wavelets as the choice of basis functions since the continuous wavelets are not orthogonal. The choice of mother wavelet in this case should allow perfect reconstruction, and basis functions such as Daubechies or Haar wavelets are good candidates. Equation (5.71) shows the C D control problem formulation with discrete wavelets in its general form. In this equation, the C D profile is first represented in terms of its wavelet expansion and the estimation of the wavelet coefficients is formulated with a least-squares relationship. The first row in equation (5.71) gives the wavelet expansion of the C D profile. The least-squares estimation of wavelet coefficients is given in the second row. The decomposition of the slice lip profile with wavelet basis functions is shown in the third row with the static input/output relation is given by interaction matrix in the fourth row. If we substitute for C D profile in terms of the wavelet representation of the slice lip in the second row, an input/output relation between wavelet coefficients of the slice lip and the C D profile is derived in the last row of equation (5.71). The next step is to use a control strategy  97  Chapter 5: Cross Directional Profile Control  to change input coefficients to achieve the desired set points.  The second approach to the C D control problem makes use of continuous wavelet transform. This method is the preferred control scheme in this thesis. The main advantage of discrete wavelet transform (DWT) over continuous wavelet transform (CWT) is the perfect reconstruction property of DWT. This feature is not important in the representation of the C D profile due to limited control bandwidth of the actuators. The controller bandwidth only allows a lowpass version of the C D profile for control. The C W T may be more efficient in this application since the D W T representation will result in significant wavelet coefficients across different scales as a result of the subsampling operation. The subsampling is not present in the CWT, and variations can be represented with fewer coefficients for the appropriate basis. The second scheme is explained in more detail in the following sections.  y(x) =  ipm,n( )  Pi> — {$m,nty ,n) m  u(x) = ip  mi  n  D W decomposition of profile  Pip  x  * V'm.n y( )  wavelet coefficient estimation for profile.  x  (x)  D W decomposition of slice lip.  y(x) — G u{x)  static input/output relation. (5.71)  =  {ipm,n' r' ) l  1  mn  V'm.n G u(x)  wavelet coefficient estimation for profile substituted for y(x).  $4,  = {^m,ntp , ) m n  ^ ^m,n  G  ^  ,nX ) x  m  i  M  relation between input/output  ^  wavelet coefficients.  In the above equation the variables are defined as:  y(x)  and u(x) = C D & actuator profile and  = Wavelet coefficients  ip  m<n  ip  mn  =  Wavelet bases  = Transpose of wavelet bases matrix  flip = Estimated wavelet coefficient G = Interaction Matrix. 98  Chapter 5: Cross Directional Profile Control  5.3  Overview of the Control Scheme A schematic of the overall profile control algorithm is shown in Figure (5.42). This algorithm  relies on identification procedures obtained in chapter 4. The C D profile control depends heavily on the form of the response, which is also an indication of the controllability of the system. The standard method of determining the response width is by performing a bump test as is shown in the first stage of the C D control procedure. It is not necessary to perform a bump test in every control cycle, but the knowledge of the response width is a requirement of the control systems and a bump test needs to be performed whenever a grade change or a major variation in one of the process variables, especially M D parameters such as slice lip opening or speed, is observed. The parameters of the response shape are identified with the scheme introduced in chapter 4. The identified parameters are the settings for the new mother wavelet when the wavelet decomposition of the total C D profile is carried out. This simply means that the correlation process is carried out using the C D profile and the response shape and its dilated and shifted versions. In the identification scheme in chapter 4, the maximum correlation value determined the parameter estimation. In the control scheme the purpose is to identify all the variations in the C D profile that are significant. This requires elimination of some of the detected variations in the C D profile using a thresholding scheme. Thresholding is an established method for data compression with discrete wavelets [25]. In fact, in the continuous case we resort to a similar approach to the de-noising scheme [25] used in the D W T to reduce the detected variations in C D profile and eventually to minimize the actuator movements. Thresholding is a nonlinear filtering strategy which sets those coefficients that are less than the threshold value to zero. The next step is to reconstruct the filtered version of the C D profile after thresholding using the wave model developed in previous chapters. The wavelet transform relates the parameters of the response model at the slice and the dry line/scanner (the dry weight does not change within the dryer section of the paper machine) by a simple scaling operation. The same argument presented in the identification part will be given here. As discussed previously, the dilation can be viewed in two different ways, one as changing the sampling length (independent variable in the C D is space rather than time) and the other as changing the model parameters while keeping the sample length constant. In the latter view, the mapping problem from scanner resolution to actuator resolution (which can also  99  Chapter 5:  MW'= Mother Wavelet  Start Bump i Test I  Cross Directional Profile Control  Set Point Smoothing/ Beam Model  Perform Bump Test Response ID with CWT  Input  New MW* Based on Response Shape  | System Identification  CWT of Profile with New MW* Parameters  CD Profile  Input CD Control = Decoupling Algorithm j Actuator ' Moves ; Output  Actuator ji Set Points; i Profile to set point Conversion  Dahlin Controller MD Control  |i  Ifl  .  .  .  Profile  JL  Figure 5.42: Schematic of CD/MD control of dry weight profile with Continuous Wavelet Transform. be viewed as a sampling rate reduction in a multirate system context) is replaced with a change in the scale parameter while the number of samples stays constant. The compressed/decoupled profile can now be fed into a time delay controller, often chosen to be a Dahlin controller. In order to correlate the parameters at the slice with those of the response shape, a smoothing algorithm is required at the slice to generate the shape of the bump. One approach for the smoothing algorithm is to use the deflection curve of the slice lip ([3], [64] and [53]).  5.4 Control Algorithm In this section, blocks shown in Figure (5.42) that have not been discussed previously are explained in more detail. There are ten different steps/blocks in the complete profile control scheme, some of which need to be calculated only once, which provide information that will be required in every control cycle. Blocks 1, 2 and 3 involve the identification of the response shape parameters as discussed in chapter 4, and the details will not be repeated in this section. These blocks are required only once at the beginning of the control cycle where the mother wavelet equation ( as found in identification stage) is calculated with the new set of parameters. The mother wavelet shape required  100  Chapter 5: Cross Directional Profile Control  for control is the response shape at the dry line. Blocks 4—10 will be given in more detail in the following sections.  5.4.1  Slice Lip Deflection Model  The deflection of the slice lip was used as the initial condition for the surface wave. The known parameters of the slice are normally the actuator settings. This requires a modelling of the slice lip as a beam deflected under loading by slice screws. A good treatment of the subject has been reported in [3], [64], [54] and [53]. The slice lip is modelled as a deflected beam based on the results presented in the above references. Here we determine the elastic curve of the slice beam based on the method of integration. The curvature of a plane curve at a point Q(x,u) is — ~ Pi  d x 2  '<M 1 2  -  (5.73)  3/2  where u(x), in the above equation, represents the deflection, and pi is the radius of curvature. If the second order variations are neglected, the curvature varies with ^ r , and equation (5.73) may be related to the governing differential equation of the elastic curve while making use of the relation between the bending moment and the curvature of the natural surface in pure bending [5].  In the above equation M(x) is the bending moment and EI is the flexural rigidity of the beam (E = Young's Modulus of elasticity and / = Moment of Inertia about neutral axis). If the variation of M(x) along the beam is known, with proper boundary conditions, u, the deflection curve can be found by integrating equation (5.74) twice. The slice is represented as a beam with free supports at the ends as in Figure (5.43).  The reactions R  a  and Rt, can be calculated in terms of P by summing  the moments at both ends of the beam.  (5.75) *  =  P  L  The bending moment along positive x-axis for a coordinate system placed at a is E I ^  = - P ^ - x  101  + P(x-(L-b))  (5.76)  Chapter 5: Cross Directional Profile Control  A u  ! i i  U-  Jr^L a  tp  3—:  Ip  a  Rbt  b  Figure 5.43: A schematic of the slice lip as a beam with loading, where (.) is defined as the singularity function and is equal to zero if the argument is less than zero and should be treated as parenthesis when the argument is positive. Assuming a constant flexural rigidity for the beam and integrating equation (5.76) twice, the relation for the deflection of the slice is given by:  U  ^  =  ~ ^TE1 P  X3  +  P  {  X  ~ "^ (JL  B  3  +  C  L  X+  (- >  C2  5  77  By using the boundary conditions at both ends, the equation for the force in each actuator c\ and c can be determined. If the actuators are flexible rods they can be modelled as stiff springs. This 2  assumption is true when thermal rods are used instead of motor driven actuators.  The relative  movement of the actuator can be described as the difference between the deflection of the rod if it were rigid and a secondary movement due to spring action (secondary movement = actuator set points). The resultant deflection can then be modeled as:  u(xi) = ui(xi)  Pi - —  (5.78)  u\(xi) represents the secondary movement of the beam at location / due to spring action and ks is the spring constant. The boundary conditions can now be applied by calculating the composite deflection at both ends and c\, c  2  are given as:  . . P(L-b) c = ui(a) + —— 2  (ft)- (q)'  ttl  ttl  UL-b),  2  102  2  ]  1  (  5  J  9  )  Chapter 5: Cross Directional Profile Control  By eliminating P between equations (5.78) and (5.77) , after substituting for the constants, one can find the deflection equation in terms of actuator set points ui(x{).  ^ { ( L  - b)(L*  -x  2  - b )x + L{x 2  - (L  Equation (5.80) shows this relation [54].  - b)) } 3  +  - b) +  ^  L(X ) 1  i + 3757(W - b)) + ^-((L-b) •tufa) - «i(a) -  i  W  u L  b)  J  + ui(a) +  + 2^6) L  (  5  8  0  )  Xi  This relation shows the deflection due to a single actuator movement.  The same concept can be  extended to the whole beam by making use of the principle of superposition. A careful look at equation (5.80) shows that the relation between set points and the slice shape can be viewed as a smoothing operation on  M_(.T.).  This becomes evident when the equation for the complete slice  is written, since the overall equation will be a weighted sum of u\(xi)  at different locations. The  process can thus be viewed as an interpolation of the actuator set points with a suitable window. The bandwidth of the window is a function of the shape of the bump on the slice. A flexible slice lip requires a narrow-width window while a stiff lip should be interpolated with a wide window. In the implementation of the C D control algorithm one needs to calculate the smooth profile from set points and vice versa. The inverse process (going from the smooth profile to actuator set points) can be interpreted as a sampling process or, in the discrete case, a sampling reduction scheme where the new sampling distance is the actuator width. This novel interpretation will reduce the computation cost of the beam model. This summarizes the main ideas of Block 4 in Figure (5.42).  5.4.2 CD Profile Decomposition with CWT The algorithm in Block 5 is the same as that used in chapter 4. The mother wavelet used for decomposition of the C D profile is the response shape. By limiting the scale parameter to that of the response shape, one is performing a lowpass filtering on the C D profile. The lowpass filtering which is embedded in the decomposition procedure is the implicit implementation of the controllability of C D actuators as a feature in CWT, recognizing that C D actuators can not control any variations narrower than the response shape and those uncontrollable variations should not be passed on to the controller. Clearly, as the response width becomes wider the controllability of the controller is  103  Chapter 5: Cross Directional Profile Control  reduced. The bandwidth of the controller is interpreted as the inverse of the response width. Another point to note is that actuators can only control those variations that are aligned with the centre of actuators, so that the phase of the response (the location of the variations) is also a limiting factor in controllability. In a Fourier domain, this is one case where both the magnitude and the phase of the response are critical in controller bandwidth. In the context of wavelet theory, information about the phase (location of the signal) is represented in the translation or shift parameter. In order to limit the variations to those in phase with actuators, coefficients that have shift parameters corresponding to the actuators should be chosen. This process can be carried out by multiplication of the wavelet coefficients along the shift axis with a comb filter. The same result will be obtained by multiplication of the reconstructed C D profile with a comb filter. By performing the C W T on C D profile and limiting the scale, a lowpass and controllable version of the profile is generated.  5.4.3 Thresholding  Wavelet coefficients of the decomposed C D profile indicate the degree of similarity of a segment of the profile to the basis function at a specific scale and translation. In the continuous case, for each scale there are at least the same number of coefficients as there are samples in the signal. Since the exact nature of the response is not known and since there are always modelling errors, a necessary element in a C D control algorithm is to filter out significant variations. A control scheme that proposes many actuator movements in order to respond to small variations is not considered an efficient control methodology. Increased variations due to controller action caused by either a poorly tuned controller or by model mismatch is not uncommon in process industries. The goal of the thresholding scheme used here is to retain those wavelet coefficients that contribute the most to the energy of the profile signal and eliminate the rest. There are two main schemes for thresholding [86]. Hard thresholding sets the coefficients with magnitudes within a band to zero while leaving those outside the band unchanged. In soft thresholding, the band is the same as that in hard thresholding but coefficients outside the band are decreased in magnitude by the bandwidth value. The graphical representation of these two schemes is shown in Figure (5.44) [86]. The thresholding algorithm is given by equation (5.81). There are different criteria for choosing the threshold level. The four more  104  Chapter 5: Cross Directional Profile Control  important selection scheme are: Stein's Unbiased Risk Estimate, fixed threshold, Minimax criteria and a combination of these schemes [86]. In the C D control algorithm the threshold level can be chosen based on the steady state noise level of the C D profile.  {  x(t),  \x\ >. 6  0,  |x| < 6  rsign(x(t)) Mag|  s o / t  \x\ > 6  = i  U, sign(x)  (\x(t) - 6\),  (5.81) a: | > 0  \x\<6  = ar| < 0  Hard Thresholding  Soft Thresholding  Figure 5.44: Schematic of the Hard and Soft thresholding algorithms.  5.4.4 Reconstruction of CD Profile for Control  After a second stage of filtering by thresholding, correlation is used to eliminate redundancy in the wavelet coefficients. A trimmed version of the C D profile can then be constructed that is not only controllable, but also minimizes the actuator movements (Block 7). The reconstruction procedure is simply that of multiplying the wavelet coefficients by the wavelet basis at the corresponding scale and shift parameters. The variations in this profile are controllable but the profile first requires further processing to change resolution from scanner to actuator level as explained in the next section. 105  Chapter 5: Cross Directional Profile Control  5.4.5 Change of Scale (Profile Compression)  By making use of the wave model in chapter 2 and the ideas in chapter 4, profile mapping reduces to a simple linear relation of scaling in the context of wavelet theory. This interpretation resolves all the issues associated with inversion of the interaction matrix and its possible ill-conditioning due to its band-diagonal characteristics [57]. The physics of the problem and the interpretation of the Fourdrinier table as a dispersive channel allow one to consider the propagation of a disturbance at the slice on the Fourdrinier table as a reduction in frequency of the generated wave. The decrease in the spatial frequency is equivalent to an increase in scale parameter in the wavelet. C D control requires the inverse process, which is reduction of the resolution from scanner to the actuator level. This translates into a decrease in scale or an increase in the spatial frequency of the response. The fundamental assumption in this scheme is that the change in scale is independent of the initial value of the spatial frequency (the bump response is assumed to be the same across the machine). This means that if there are two different disturbances at the slice with spatial frequencies of k] and that change to ks and ^ at the dry line respectively, the following holds:  r =r «3  ( 5  -  8 2 )  A, 4  This assumption is not limiting for the algorithm, and the scheme is valid even if the scale parameter varies at different C D locations. The only requirement in the latter case is that bump tests should be performed at those different C D location at which the response shape is expected to vary. By making use of the bump test results the initial spatial frequency at the slice can be related to the response shape at the dry line by simple division. This ratio is used when the lowpass profile at the scanner resolution is changed to actuator resolution. In this case, instead of having a single scale ratio between spatial frequency at slice and the dry line, several scale ratios exist for different regions across the machine. In transforming the profile from scanner resolution to actuator resolution the number of samples stays constant but the width of profile is reduced by the scale factor. The term compression refers to this change of scale and should not be confused with the terminology in data communication (data compression). The C D profile at the actuator resolution is subtracted from the  106  Chapter 5: Cross Directional Profile Control  desired profile to generate the error profile. The error profile is then sampled to give the actuator set points. The decoupled set points are transferred to the M D control block.  5.4.6 Dahlin Control  The C D control task is completed after the decoupling of the C D profile is carried out. The decoupled set points are then transferred to the M D controller dedicated for each actuator (Block 10). One of the characteristics of the M D dynamics is the time delay due to transportation lag for the sheet to reach the dry end of the paper machine. The time delay in current paper machines, although heavily dependent on the machine speed, is approximately between 40 seconds to one minute. The common approach to modelling the M D dynamics is as a first-order lag plus time delay. In many installations of computer control systems on paper machines today, the dominant control algorithm is the Dahlin controller, although a few multivariable controllers for M D have been reported recently [94]. The Dahlin controller is a model following controller that assumes the closed loop response of the desired system is a first-order plus a time delay. The standard formulation of the Dahlin controller is shown in equation (5.83)[79]. The sum of the coefficients of the denominator of Dahlin controller is equal to zero indicating that the denominator of GDC is always divisible by (l - z~ ), x  and so includes an  integrator. Although an integrating element in the controller eliminates steady state error, and the controller guarantees the model following feature, potential windup of the controller must be avoided. Although implementation of thresholding algorithm avoids most of the small errors being passed on to the controller, still an antiwindup algorithm is required. In addition to the implementation of the anitwindup algorithm explained in [54], other methods such as those suggested in [79] can be applied.  Desired 1st — order response  1 - D z-i p  - (1 - DJz-"-  (1 - D ) a  -  g  Jl-Daz- ) 1  P  Kg{l -D )  K (l  1  1 - Dpz-  1  - (1 -  107  D )zv  D )z~ -i N  a  (5.83)  Chapter 5: Cross Directional Profile Control  where: D  = e  D  = filter pole for the system to be controlled  p  a  ''  discrete filter pole  Kg = system gain HG{z)  = Z O H equivalent of system  D(z) = controller transfer function GDC = D a h l i n controller A = Desired closed loop time constant. 5.5  Simulation of C D Control Algorithm with Industrial Data  The main stages of the control algorithm are now applied to the bump test results shown in plot (a) and (b) in Figure (5.45). Although it is not uncommon to use a bump profile for C D control evaluation in the field, this is not a typical profile to be seen on a daily basis. The purpose of this section is to demonstrate the proposed control algorithm in a realistic situation. The C D profile is the gain or the steady state value of the M D control scheme.  Since the focus of this thesis is on  C D control, the implementation of the Dahlin controller is not included. In Figure (5.45) all the control actions are shown at the final stage after the M D transients have died out. Figure (5.45) shows the main steps involved in the proposed C D control algorithm based on continuous wavelet transform (CWT). Plots (a) and (b) show the C D profile to be controlled. Plot (c) is the average of the first 8 scans before the bump test was performed. This plot is an indication of C D variability under normal conditions before the system is put in manual mode to perform the bump test. This graph serves two purposes. First, it is an indication of system baseline noise level under control, or the controller bandwidth with specific settings for that particular day. Secondly, it is used to determine the threshold level when compression of wavelet coefficients are considered. Blocks 1—3 in Figure (5.42) require the identified response shape parameters for performing the CWT. Since this information is already available as identification results in chapter (4), they are not shown in Figure (5.45). Applying the response shape parameters to the mother wavelet equation, one can perform the C W T of the profile. The wavelet coefficients are shown in plot (d). The scale range was between 1  108  Chapter 5: Cross Directional Profile Control  and 40 with 100 divisions. Block 5 of the C D control algorithm is the thresholding or compression stage. For compression, hard thresholding is the most common strategy to use because the threshold level can be determined globally. For this purpose a threshold level of 5 was considered. This value is based on results of correlation of baseline noise with the mother wavelet. After thresholding it is necessary to extract the features of the thresholded wavelet coefficients. By locating the maximum correlation in each region, the coefficients that contribute to the signal the most are found. The results are shown in plot (e). By reducing the threshold level, the number of regions at which features are extracted increases, and the controller attempts to remove smaller variations. Plot (f) indicates the reconstructed profile based on extracted wavelet coefficients. The distance between large variations was linearly interpolated to reach a smooth reconstructed profile. The lowpass reconstructed profile captures the main features of the measured profile. This reconstructed profile is at the scanner scale. In order to convert it to actuator set points the reconstructed profile needs to be rescaled to the slice lip resolution. Since the scale parameter for this data set was calculated in chapter 4, the same value is used now to rescale the C D profile to the actuator resolution. Plot (g) shows the reconstructed profile before and after rescaling. In this plot only the scale operator has acted upon the profile and the gain is the same for both profiles. This was carried out to give a better visual feel for rescaling. The rescaled profile after gain adjustment is shown in plot (h). The rescaled profile indicates the slice lip shape and not the actuator settings.  As mentioned earlier, the process of moving from  actuator set points to slice shape can be modeled as a smoothing operation. In this case, we need the reverse procedure which is a sampling task. Another way of interpreting this operation is by a staircase approximation of slice lip profile with Haar wavelets. The actuator set points will be the approximation of the shape at the scale determined by the actuator width. Figure (5.46) shows the effect of implementation of C D control on the profile. The residuals show that the variations are at the same level or smaller variations as compared to plot (c) in Figure (5.45). There are still some low frequency variations with magnitude close to the noise level. This can be adjusted by changing the threshold level as explained before. The second graph in Figure (5.46) shows the control action required to reduce the variations assuming a flat profile is desired.  109  Chapter 5: Cross Directional Profile Control  5.6 Summary In this chapter, results from previous chapters were combined to develop a novel C D control algorithm based on the continuous wavelet transform. It was pointed out that dry weight control or the control of any other quality variable of the sheet consists of essentially two different schemes. C D control is mainly a decoupling algorithm and, after decoupling, set points are passed on to a series of single-input single-output time delay controllers such as Dahlin algorithm. Since the information about the width of the response is critical to the design of a good C D controller, the identification algorithm previously introduced is used to find the parameters for the mother wavelet required for C D control algorithm. By limiting the scale parameter an inherent lowpass filtering is performed on the C D profile. This lowpass filtering is essential in order to avoid a saw-tooth actuator profile.  Thresholding was used as a compression technique for the wavelet coefficients, and also  helped to limit the actuator movements by blocking some of the low energy variations. The main feature of the C W T based approach to C D control is the elimination of the mapping problem which requires the inversion of a large interaction matrix. The mapping in this case reduces to a simple scaling operation on C D profile between the scanner and actuator resolution. The proposed control scheme was demonstrated by applying it to a set of industrial data and the results showed a good performance of the control scheme.  110  Chapter 5: Cross Directional Profile Control  CD profile to be controlled  100  200 300 CD position (a) Actuator set points  5  10  15 20 25 CD position (b)  Localizing main CWT coefficients after thresholding 201  200 300 400 Shifts (e) C D profile reconstruction at measurement scale 100  400  30 35  100  Steady state profile with noise band  100  200 300 400 CD position (c) CD profile CWT coefficients  Figure 5.45:  200 300 Shifts (d)  400  Comparison of CD profile at actuator & slice scale  100  100  200 300 CD position (f)  200 300 CD position (g)  400  Actuator set points & scaled slice lip profile  100  400  200 300 CD position (h)  Different stages of CD control algorithm shown through plots (a) to (h).  Ill  400  Chapter 5: Cross Directional Profile Control  Figure 5.46: Residuals after the implementation of CD control and control action.  112  Chapter 6: Conclusions and Future Work  t  Chapter 6 Conclusions and Future Work  6.1 Conclusions  Results of research on paper machine wet end modelling for dry weight, identification of the model parameters, and control based on the identified model were presented. The main purpose of the modelling effort was to develop a parametric model that is a good representation of the physical process in as much details as allows a closed form solution. The model combines the surface wave theory with the drainage principles to represent the physics of the fibre distribution profile. Since the C D control plays a major role in the production of a good quality paper, it became the focus of this study. The following points are the main conclusions drawn from this research:  •  A survey of modelling strategies for the wet end of the Fourdrinier was presented. It was stated that so far there have been two main approaches to modelling of the wet end, first, an approach using fundamental fluid mechanics to obtain numerical solutions, and second, approximation of the process response using time series. The approach presented in this thesis for the first time takes advantage of both techniques and shows that the physical principles can lead to a better and more compact model as compared to standard modelling schemes used for C D control.  •  The model considers the disturbance at the surface of the slurry as a small amplitude wave. The combination of the surface wave with the drainage theory resulted in a damped dispersive wave that extends the one-dimensional drainage theory to two dimensions. The surface wave theory applies only up to a certain point. The limiting factor in this case is the fluid depth on top of the wire. Since the developed model is a two dimensional dispersive wave it can not represent rush/drag phenomenon.  •  The model can be used to simulate different design scenarios by changing drainage characteristics of the Fourdrinier table and changing the slice lip opening and headbox consistency. Variable consistency headboxes can be modelled. There are four parameters in the model, each with a  113  Chapter 6: Conclusions and Future Work  physical representation. Variation of these parameters can give insight to the operation of paper machines. •  Quality variables of the sheet are measured at the scanner, and a close look at the data processing scheme at the scanner was carried out. Different filtering schemes were tested and it was shown that FIR filtering is the best choice due to its ease of implementation and linear phase property. The choice of filter type was not critical in the final result when the data points were decimated by a large factor.  •  The identification of model parameters was carried out with a novel approach using the continuous wavelet transform. The continuous wavelet transform proved to be efficient in identifying model parameters. The parameters for two different types of response shapes, a single mode and a bi-modal response, were estimated with simulation using industrial data.  •  The continuous wavelet transform requires a generic basis function from which the other bases are generated. The model developed previously was shown to be a candidate mother wavelet.  •  Based on the continuous wavelet transform, a novel C D control algorithm was proposed. The control algorithm requires the response shape as a starting point.  The identified parameters  were used in the control scheme. Hard thresholding was used to reduce the number of wavelet coefficients and minimize actuator movements. The mapping issue in C D control was resolved with a simple linear scaling relation as a result of the continuous wavelet transform. The scale limits and resolution and the thresholding band are the ciritical parameters in the proposed control algorithm.  A summary of the main contributions of this thesis are:  1.  Although there are excellent results on paper machine modelling available, the model developed in this thesis is the only known parametric physical model of the C D response on paper machines that is developed for the C D control purpose.  2.  The main application of the wavelet transform has been in signal and image processing and its use in the context of system identification in general, and in the paper machine area in particular, is unique to this thesis.  114  Chapter 6: Conclusions and Future Work  3.  The proposed C D control algorithm based on wavelet theory is a novel approach with the benefit of eliminating the mapping problem in C D control, and its associated issues such as ill-conditioning.  6.2 Future Work Some issues have been left open ended. Some interesting follow ups from the work presented here could be in the following areas. Dilution Headbox Modelling: In the modelling area, there are a significant number of installation of consistency profiling headboxes that demand some modelling effort. The developed model is capable of simulating headboxes with different consistencies across the machine. One of the issues that is worth investigating is the mixing of adjacent jets with different consistencies on the Fourdinier table. Correlation of Fibre Orientation Angle and Dry Weight Profile: A problem that is closely related to mixing of jets of different consistencies is the relation between fibre orientation angle and the C D weight profile. One of the main ideas behind the development of consistency headboxes was to reduce the effect of slice lip deflection on formation. A n extension to the developed model in this thesis can be used to look at the correlation of fibre orientation angle and the weight profile. Rush/Drag Modelling in Terms of Variable Physical Parameters: One of the limitations of the developed model is the absence of an explicit modelling of the rush/drag phenomenon. The rush/drag process introduces an element of complexity in the model since the current two-dimensional model needs to be extended to three dimensions.  Investigation of effect of rush/drag on the drainage  characteristics of wire is a worthwhile extension to the developed response model. C D Control with D W T : Extensions of the proposed control strategy can be considered. The general formulation of C D control in a D W T framework was presented in Chapter 5.  The potential use  of discrete wavelet transform (DWT) for C D control is a natural extension of this work. Due to possible difficulties with limited translation parameter in DWT, application of wavelet packets is a worthwhile investigation.  115  Appendix A  A . l Solution of Laplace's Equation by Separation of Variables The following is a summary of the methodology presented in [18], [24] and [8].  Laplace's  equation in two dimensions is denned as follows:  d  w  +  =°  x and y in the above refer to independent variables.  -  (A 85)  The standard method of solving Laplace's  equation is by separation of variables technique. In this approach, we assume that the solution can be written as a spareable solution of three different functions (a function of x, y and time) <t> = X(x)Y(y)T(t) The function T(t), represents the dependence of solution on time variable.  (A.86) This dependence is  normally expressed in terms of lateral boundary conditions. The standard lateral boundary condition is a periodic function of time. The following defines the lateral boundary condition. (f>(x, t) = 4>(x + L,t) (A.87) <f>(x,t) = <f>(x,t + T) Since 4> has to be periodic in time by the lateral boundary condition, T(t) can be defined as: T(t) - sin (tot)  (A.88)  Substituting in equation (A.85) form (A.88) and (A.86) and dividing by 4>, we reach the following result.  dx  X  2  Y  dy2  Since the left and the right-hand sides of the second row in the above equation are independent of one another, both terms must be a constant. dX 2  1  dr  X  ax  A  2  Y  dy2  = X2  116  ^  A  9  (  )  ^  Value of separation constant  Differential Equation  Solutions  A  dX — - +A X = 0 ax dY  Real  2  x  2  2  X(a-) = A cos  (A.T) +  y(j,) = Ce  + i? -  B sin (Az)  z  A > 0 2  Xy  2  A y  e  = Ax + B  d?X dx  2  A= 0  Y(y) = Cy + D  _ o  c J 2 y  dy  2  dX 2  A <0 2  , A = j|A|  1  dx l  P  1  |  X(x)  2  ~  2  dy  +  |  x  Y{y) = C cos | A \y + D sin |A| y  2  2  + Be~^  x  = °  + IAI Y - 0 - °  Y  = Ae^  A  | 1  Table A.12 Possible solutions of Laplace's equation depending on the value of separation constant from Reference [24]. The expressions above can now be written as: dX 2  + \X 2  dx dY 2  =0 (A.91)  2  - A y = o 2  dy2  Both equations are second order differential equations with the general solution given in the following: X{x)  = Ii\ e~  lXx  +K  e  jXx  2  (A.92) Y(y) = C y x  ie  + C e~  Xy  2  Depending on the value of A, the separation constant, different solutions can exist. Table (A. 12) shows the general form of the possible solutions for equation (A.92) [24].  When A > 0 2  &  Real, the  general form of the potential function (the solution of Laplace's equation) is given as:  (Ce  Xy  + De~ ^ Xy  (A sin (\x) + B cos (Xx))  117  sin(wi)  (A. 100)  Since the superposition of the solutions is also a solution of Laplace's equation, after some algebra, a different form of the potential function results. 4> = ( A i cosh(Xy) + Bi sinh(Xy))  cos(Xx-u)t)  (A.101)  By applying the boundary conditions at the free surface and the bed, A can be defined in terms of the depth h, wave number k, and damping factor kj.  A.2 Signal Processing Scheme of the Measurex Scanner The elements of the Measurex scanner were shown in the system diagram of Figure (3.18) in Chapter 3. The signal processing and data analysis associated with the scanner is carried out in the control computer and the V F C module (voltage to frequency converter). The major steps involved in processing the incoming data in the Measurex system are : Analog to digital conversion in the V F C . •  Digital filtering prior to data compression.  •  Data compression or many-to-one mapping.  •  Calibration.  •  C D profile calculation.  Digital filtering and data compression are done by an averaging process in a single step; however, it is convenient to separate the digital filtering from the data compression.  A.2.1 A/D Conversion in VFC The analog voltage has to be converted into digital form for processing in the computer. In the Measurex scanner this is done by a voltage-to-frequency converter (VFC). The V F C is particularly suitable for such applications since it has a high noise immunity. This type of A / D conversion is known as an integrating A / D technique and operates by an indirect conversion method. The function of the V F C is to generate a series of pulses, the frequency of which is proportional to the input voltage or current. In the scanner, the current is first converted to voltage and then fed to the V F C . The input signal is averaged over a certain time (integrated) and is used as a single measurement 118  V;in ,  O  R1  WW  Plxintegrator S  precision pulse timer  —O  X  pulse output  Figure A.47: Schematic diagram of a voltage to frequency converter.  point. The averaging process will introduce some lowpass filtering of the noisy analog signal and explains the noise immunity of the V F C . A schematic of a typical V F C is shown in Figure (A.47). The A / D conversion technique is known as charge-balancing where a single capacitor is used to keep track of the ratio of the input voltage to a reference voltage. The output of the V F C is proportional to the average input voltage over a fixed time and can be formulated as the following  f = -  where /  fe  1  V  R  T  V f  #i  in  re  2  r are the output frequency and the pulse width respectively.  (A. 102)  In summary, the A / D  conversion by V F C will introduce some amplitude attenuation due to lowpass filtering.  A.2.2 Data Compression In each scan the scanner collects approximately 600 data points. These points are then averaged over every 7 points to obtain 83 databoxes and displayed on the operator's console. The averaging process, in the context of a multirate system, can be viewed as a decimation or subsampling process. Averaging is a lowpass filtering with a sliding rectangular window and a width equal to the sequence compression ratio (7 to 1 in this case). Figure (A.48) shows the schematic of an averaging process. The equivalent of the averaging process in the multirate systems will be a decimator preceded by a rectangular lowpass filter as shown in Figure (A.49). If sequence A is represented by A(n), where  119  sequence a  • 4 A i . wW ' / I \ \ \ \  < 4 l  Sequence A represented with 3 points by averaging , 1 1 Rectangular window |  sequence A  Sequence A with 21 points  Figure A.48: Schematic of the averaging process viewed as a mapping of many-to-one process.  Decimation block  Prefiltering with a rectangular window  Subsampling by 7  Figure A.49: Straight averaging viewed as a decimation block. n — l....m,  the averaged sequence will be represented by 1  a(k)  m  l=m  (A.103) .=0  The filtering can be improved if a weighted sum is used for averaging and if over lapping methods are used. After this stage of processing, the signal is ready for the calibration process. Since the basis weight signal was considered in the experiment, only the calibration of the basis weight signal is explained in the next section.  A.2.3 Basis Weight Calibration The measured voltage signal has to be converted to basis weight units for display on the operator's console. The calibration is carried out using a nonlinear curve. The x-axis of the calibration curve is the basis weight, and the y-axis is a variable called the calibration  or transmission ratio. A typical  calibration curve is shown in Figure (A.50). This curve is sometimes referred to as  Transmission  curve and can be approximated by a negative exponential function. The calibration ratio is defined as:  R  air sheet  measurement measurement  120  (volts) (volts)  (A. 104)  A  / \  v  A  ~-jjc  Clean curve  Dirty curve  \  -  it  /  >;  —Expected range of R  /  ^  ^  ^  ^  Range of basis weight Basis weight  Figure A.50: A typical transmission curve used in Measurex calibration scheme. For each grade in production, the highest and the lowest values of R are calculated and using the transmission curve, the upper and lower bounds of the basis weight are determined. The region identified by the calculated bounds can be used to linearize the curve if need be. The calibration scheme in the Measurex system uses a polynomial fit to the transmission curve, and for each grade a different set of parameters is used. This scheme is the first stage of calibration. The Measurex algorithm accounts for four different disturbances that could introduce some error in the basis weight measurements. •  The corrections taken into accounted are:  Dirt factor. Temperature factor.  •  Position alignment factor.  •  Z direction factor.  As the dirt builds up around the source, receiver and the air gap, the intensity of the beam collected by the receiver body is reduced and introduces error in the calibration. The error introduced by the dirt factor is compensated by another curve, the "dirty calibration curve" that is used together with the original calibration curve. This curve is also represented by a polynomial, the coefficients of which are accessible through the Measurex computer. Temperature is another source of error. The density of air changes inversely with temperature and results in a varying transmission ratio. The location of the databoxes (slice positions) may also be a source of error due to scanner speed variations. The  121  Z direction factor accounts for the variation in distance between the source and the receiver as the scanner traverses the sheet. These correction factors are available for each grade of paper and have been used in our calibration of the basis weight signal.  A.2.4 Moisture Calibration The operation of moisture sensors is based on infrared radiation. The source emits infrared radiation that is reflected through the sheet several times before entering a beam splitter and then being guided towards two infra red filters. One filter passes only a wavelength of 1.8 micron, and the other is sensitive to 1.9 micron. A wavelength of 1.8 microns is used as a control reference, and 1.9 micron is the wavelength that is heavily absorbed by water present in the sheet. A weak current corresponding to the amount of detected radiation for each wavelength (1.8 and 1.9 micron) is fed to the V F C , and equivalent frequencies are sent to the computer. The ratio of these frequencies (1.8 f/1.9 f) is used in the calibration process.  The calibration curve is a linear curve, the slope  and the offset of which need to be determined. The calibration parameters are grade dependent and are specified for a range of moisture levels. The calibration procedure is carried out by taking grab samples of the sheet which are immediately sealed and weighed in the lab. The samples are then placed in an oven for a day to remove moisture and weighed again immediately after removal. The difference in weight before and after the drying step indicates the weight of the water. Figure (A.51) shows a typical moisture calibration curve. The parameters of the calibration curve are adjusted to match the moisture readings of the scanner with those of the laboratory tests.  A.3 Definition of Norm Definitions presented in this section are taken from [80] and [58]. Before the norm operator is defined, a summary of different spaces that a norm is associated with is defined. The spaces that are common in the control theory have a certain structure that is presented in the following. 1.  A linear space is a set where the addition of elements and multiplication of an element by a scaler is properly defined, and the results of those operations always belong to the set. 122  ii  *  Ratio R=(1.8 f/1.9 f)  Figure A.51: A typical calibration curve used in Measurex scanners. 2.  A space where the distance of the elements has been properly defined as: d(yi,y )  = 0  d{yi,y )  = d(y , yi)  d(yi,y )  < d{y yz)  2  2  2  Vi = V2 & > 0 if  2  (A.105)  2  u  yi^y  + rf(y ,2/i) 3  is called a metric space. 3.  A set which is a linear space and consists of open sets including the empty set with proper definition is a topological vector space.  4.  A topological vector space where the norm is defined and is also complete is called a Banach space.  5.  A Banach space with an inner product defined and the norm is derived from the inner product is called a Hilbert space.  The norm will be more precisely defined in the next paragraph.  Loosely defined, a norm is an  element of a space which is a measure of magnitude of the element and is the generalization of the absolute-value concept. A useful concept that is an indication of the length of a vector, a matrix, a signal or a system is a function called norm denoted by the symbol ||-||. The most common norm is the Euclidean norm which is the distance between two points in space. A norm should have the following properties: 1.  ||e|| > 0. 123  2.  ||e|| = 0 & e = 0.  3.  \\a • e\\ = \a\ • \\e\\.  4.  ||ei + e || 2  < ||ei|| + ||e ||. 2  In this section only the vector-norm and (/,Z)-norms are defined.  If we consider a vector b the  general /?-norm for this vector is defined as the following:  S>r)  I/P (A. 106)  When p=l, 2 and oo, the /?-norm is referred to as sum-norm (1—norm), Euclidean-norm (2—norm) and infinity-norm (oo-norm) respectively.  The infinity-norm is the largest magnitude in the vector  (to be precise the max operator should be replaced with sup operator). A Hilbert space that the inner product is defined for infinite sequences is referred to as l -space. v  The / -norm is defined for p  signals and is given by:  (  oo  J  \  1/P  dr\  ^(T^  (A.107)  If the signal is defined as an infinite sequence, the norm is represented with lower case /. In the case of continuous signals ( functions instead of infinite sequences ), the norm is represented by upper case L. The 2-norm in this case refers to the energy of the signal, and, if it is bounded, the signal is called an energy signal: iSKJ  won = 2  (A. 108)  \  The wavelet frame bounds in chapter 4 were defined as  L2-norms.  A.4 Calculation of Admissibility Function The admissibility condition for a function tp to be a mother wavelet was defined in chapter 4 as:  du> <  124  00  (A. 109)  In order for the model developed in chapter 2 to be a candidate for mother wavelet it has to satisfy equation (A. 109). The first step is to calculate the Fourier transform of the model function. The following shows the calculation of Fourier transform. -f oo  ^ ( w ) = J e~ W cos(kx)  e~ dx  a  lwx  — OO + 00  = J e~  cos (kx) e-  ax  dx  JWX  0  0 < x < + 00  + 0  J e  cos (kx) e~  ax  dx  JWX  - 00 < x < 0  ( A . l 10)  -co  + JUJ  a  (a + JOJ)  0< x  +  2  < +00  k  2  + a ~ jv (a - jcj) The admissibility function  — co < x < 0  + k  2  2  can now be calculated. Insertion of the above in equation (A. 109) gives: Q  00  (a +a; )  I  (a +fc -a; ) + 4 a ^  _  ^  2  2  2  2  2  2  u  J  Lip =  2  2  / -( a  j  2  '  +fc +u 2  du  —00  The calculation of limits in equation (A.l 11) when u> -* ^ 0 0 results in indefinite forms  (00—00  and  2|) and the integrals are solved by resorting to L'Hospital's rule. co  v  l i m  2  dw  2  2  2  ,  2  du>  0 J 0  = A  W )  2  2  I -  C,b = -r~  -++00  lim C,p a;->-oo  ( q  [ ( a + A: -u; ) +4a a,  d  OJ 0 0  J/  (a +a. ) (a +a. ) (  a  2  2  2  2  2  +  W  r  —LO  rfc  W  —00 co  v  f  l i m Q, = J  H w->-co  r  0  n J  /  d  2  2  2  2  2  a  —  div  1  dw  2  I (a + A: -w ) +4 u;  /  ^  (A. 112)  (a +a; ) 2  /  /  \  -  L  -  2  ,  ,  I  —  du  < 00  OJ  i  ( a  2  + ^ )  M  —to  -CO  125  ^  ^ J  <  0  O  The last two rows of equation (A. 112) prove the admissibility of the developed model as a mother wavelet function.  126  Bibliography  [I]  K. J. Astrom. Computer Control of a Paper Machine - an Application of Linear Stochastic Control Theory. IBM Journal, July Issue, 1967.  [2]  LP. Bachand and S.M Cadieux. Pressure Profiles of Modern Drainage Elements — Effect on Water Removal. In Proceedings of the Technical Association of the Pulp and Paper  Industry,  1982 Annual Meeting. [3]  R. Balakrishnan and D. McFarlin. The Application of a Model of the Slice Lip as an Elastic Beam for Weight Profile Control. Tappi Process Control Symposium, pages 105-109, 1985.  [4]  A. E . Beecher and R. A . Bareiss. Theory and Practice of Automatic Control of Basis Weight Profiles. Tappi Journal, 53(5), pages 847-852, May 1970.  [5]  F.P. Beer and E.R. Johnston. Mechanics of Material.  McGraw-Hill, Inc. 2nd Ed., New York,  N Y , U S A , 1992. [6]  L . G . Bergh and J.F. MacGregor. Spatial Control of Sheet and Film Processing. The Journal of Chemical Engineering,  [7]  Canadian  Vol. 65, pages 148—155, 1987.  W. L . Bialkowski. Newsprint Uniformity and Process Control — The Untapped Competitive Edge. In CPPA Conference, Montreal, Canada, Feb. 1990.  [8]  W . E . Boyce and R.C. DiPrima. Elementary Problem.  [9]  Differential  Equations  and Boundary  Value  John Wiley & Sons, Inc. New York, N Y , U S A , 1986.  T. J. Boyle. Practical Algorithm for Cross-directional Control. Tappi, Vol. 61, pages 77—SO, 1978.  [10]  R. M .R. Branion et al. A Current State of Drainage Monograph.  [II]  Theory up to JANUARY  1978, A  Omni Continental LTD. New Westminister, B . C . Canada. 1978.  G. Burkhard and P. E . Wrist. The Evaluation of Paper Machine Stock System by Basis Weight Analysis. Tappi, Vol.37, 1954.  127  [12]  G. Burkhard and P.E. Wrist. Investigation of High Speed Paper Machine Drainage Phenomona. Pulp and Paper Magazine  [13]  of Canada, pages 100-118, March 1956.  S.M Cadieux and J.P. Bachand. Comparative Pressure Profile Measurements on Modern Drainage Elements. Pulp & Paper Canada 84: 12, pages T258-T262, 1983.  [14]  S. C. Chen. Cross Machine Profile Control for Heavy Weight Paper. EUCEPA,  Stockholm,  pages 77-94, May 8-11, 1990. [15]  S. C . Chen and R. G. Wilhelm. Optimal Control of Cross-machine Direction Web Profile with Constraints on the Control Effort. ACC, Vol. 3, pages 1409-1414, 1986.  [16]  N . P. Cheremsisinoff. Pocket  Handbook for  Solid-Liquid  Separations.  Gulf Publishing  Company, Houston, Texas, U.S.A., 1984. [17]  R. Coifman, Y . Meyer, and M . V . Wickerhauser. Wavelet Analysis and Signal Processing. Ruskai et al, Jones and Bartlett, Boston., pages 153—178, 1992.  [18]  G. D. Crapper. Introduction to Water Waves. Ellis Horwood Limited, Chichester, West Sussex, P019 1EB, England., 1984.  [19]  R. E . Crochiere and L.R. Rabiner. Multirate  Digital Signal Processing.  Prentice-Hall Signal  Processing Series, Prentice-Hall Inc. Englewood Cliffs, New Jersey 07632, 1983. [20]  W . H . Cuffey. Some Factors Involved in Basis Weight Uniformity. TAPPI Journal  40 (6),  pages 190A-196A, June 1957. [21]  I. G . Currie. Fundamental Mechanics  of Fluids. McGraw Hill, U.S.A., 1974.  [22]  K . Cutshall. Nature of Paper Variations. Tappi Journal, pages 81—90, June 1990.  [23]  Xue-dong Dai, R. L . Motard, and B. Joseph. Introduction to Wavelet Transform and TimeFrequency Analysis. In Wavelet Applications Kluwer Academic Publishers,  [24]  in Chemical Engineering., Ed. Motard & Joseph,  Ma, U.S.A., pages 1—32, 1993.  G. R. Dean and A . R. Dalrymple. Water Wave Mechanics for Engineers and Scientist. PrenticeHall, Inc., Englewood Cliffs, New Jersey 07632, USA., 1984.  [25]  D . L . Donoho. De-Noising by Soft-Thresholding. Preprint, Department of Statistics, Stanford University,,  1992.  [26]  G. A . Dumont. Lecture Notes on Paper Machine Estimation and Control. Department of Electrical Engineering, University of British Columbia, Canada., 1987.  [27]  G. A . Dumont, K Natarajan, and M . S. Davies. An Algorithm for Estimation of Cross and Machine Direction Moisture Profile. Technical report, Post-Graduate Research Lab Report PGRL 388, Pulp and Paper Research Institute of Canada, 1987.  [28]  S. R. Duncan. Estimating the Response of Actuators in a Cross-Directional Control System. Control Systems 96, Halifax, Nova Scotia, Canada, pages 19-22, April 30-May2, 1996.  [29]  L . Eriksson, J. Hill, and I. Lundqvist. Process Control in Groundwood Production and Control of Quality Profiles on Paper Machines. Das Papier 32, No. 10A, pages 157—164, Nov. 1978.  [30]  S Fielding and L . Baldwin. Establishing Drainage Curves on the Forming Table. Paper Technology, pages 28—30, Oct. 1992.  [31]  J. Ghofraniha, M.S. Davies, and G.A. Dumont. CD Response Modelling for Control of a Paper Machine.  [32]  The 4th IEEE Conference on Control Applications, Albany, N . Y . , Sept. 1995.  D. Gobor. Theory of Communication Part I. The Institution of Electrical Part III, pp 429-441,  [33]  Engineers,  Vol. 93.  1946.  G. C. Goodwin and K. Sin. Adaptive Filtering  Prediction  and Control. . Prentice Hall, Inc.  Englewood Cliffs, New Jersy, U S A , 1984. [34]  G.C Goodwin, S.J. Lee, A . Carlton, and G. Wallace. Application Coating Mass Estimation. IEEE Conference A C C .  [35]  of Kalman Filtering  to Zinc  1994.  A . Graser. Cross Profile Control - Predetermination of the Control Results Regarding Slice Lip Bending. Pulp & Paper Canada 96:1, pages 35—39, 1995.  [36]  S.I Green and R.J. Kerekes. Numerical Analysis of Pressure Pulses Induced by Blades in Gap Formers. In Proceedings of Tappi Engineering  [37]  Conference, pages 185—192, 1996.  M . J. Grimble and J. Fotakis. The Design of Shape Control Systems for Sendzimir Mills. IEEE Transactions on Automatic Control, 27:, pages 656-666, 1982.  [38]  R. Hall. Coordinated Cross Machine Control. Paperija 241-248, 1991.  Puu, Paper and Timber 73: 3, pages  [39]  S. T. Han and W. L . Ingmanson. A Simplified Theory of Filtration. Tappi Journal, No. 4, April  Vol. 50,  1967.  [40]  S.T. Han. Drainage in a Vacuum Former. Tappi, Vol. 45, No. 4, April 1962.  [41]  W.P. Heath. Self-Tuning Control for Two-Dimensional Processes. PhD thesis, Control Systems Centre, UMIST, 1992.  [42]  Robert. L . Heider. How Industrial Digital Controllers Can Cause Process Cycling. American Control Conference, Arlington.  [43]  Chicago, pages 375—378, 1982.  D. Hunter. Papermaking, The History and Technique of an Ancient Craft. Dover Publication, New York, N . Y . , 1978.  [44]  W. L . Ingmanson. Filtration of High-Consistency Fiber Suspensions. Tappi Vol. 47, No. 12, December  [45]  1964.  B. Jawerth and W. Sweldens. An Overview of Wavelet Based Multiresolution Analysis. Technical report, Dept. of Mathematics, University of South Carolina., 1993.  [46]  I. J. Jonsson. Estimation and Identification of Moisture Content in Paper. Master's thesis, Dept. of Electrical Engineering, U B C , Canada, 1991.  [47]  H . Karlsson and L . Haglund. Optimal Cross-Direction Basis Weight and Moisture Profile Control on Paper Machines. In The Third International Pulp and Paper Process Symposium,  [48]  Vancouver, BC, Canada, May 2-4,  Control  1983.  H . Karlsson, I. Lundqvist, and T. Ostsman. Principles and Potentials of Cross Direction Basis Weight Control. In Proceedings of EUCEPA  Symposium, Stockholm, Sweden, May 11—14,  1982. [49]  G. Kastanakis. Prediction and Verifying the Best Quality Control Under Particular Conditions. Tappi Journal, 73(12), pages 203-207, Dec. 1990.  [50]  G. Kastanakis and A . Lizr. Interaction between the Machine Direction and Cross Direction Control Processes in Papermaking and Plastics Machines. Tappi Journal, pages 77—83, Feb. 1991.  130  [51]  M . A. Keyes. An Integrated Digital System for Distributed Parameter Paper Machine Control. ISA Trans. Vol. II, No. 2, 1972.  [52]  M . A. Keyes. Computer Control Census. Tappi Journal, 58 (6), pages 71—74, June 1975.  [53]  A.P Kjaer, M . Waller, and P.E. Wellstead. Headbox Modelling for Cross-Directional Basis Weight Control. ISA Transactions, 33, pages 245-254, 1994.  [54]  K. Kristinsson. Cross Directional Polynomials.  Control of Basis Weight on Paper Machines  Using Gram  PhD thesis, Dept. of Electrical Engineering, University of British Columbia,  1994. [55]  K. E . Kwok and Ping Li. M D - C D Interaction Modelling and Control of Paper Machines. The Canadian Journal of Chemical Engineering,  Vol. 75, February., pages 143—151, 1997.  [56]  Sir Horace Lamb. Hydrodynamics.  Dover Publications Inc., New York, N.Y., 1932.  [57]  D. L . Laughlin and M . Morari. Robust Performance of Cross-directional Basis Weight Control in Paper Manufacturing. ACC, Vol. 3, pages 2122-2128, 1989.  [58]  J.R. Leigh. Functional Analysis and Linear Control Theory. Academic Press Inc. (London) Ltd. 24/28 Oval Road, London N W L , 1980.  [59]  C. Lindeborg. A process Model of Moisture Variations. Pulp and Paper Canada, Vol. 87, No. 4, pages T142-147, 1986.  [60]  S. Mallat. A Theory for Multiresolution Signal Decomposition: the Wavelet Representation. IEEE Trans. Pattern Analysis and Machine Intelligence, 11, pages 674-693, 1989.  [61]  M . J. McCann. Modeling and Control of Distributed Parameter Systems. Technical report, SCR Report 124-C-67-55, Case Institute of Technology,  [62]  1967.  H . Meyer. A Filtration Theory for Compressible Fibrous Beds Formed from Dilute Suspensions. Tappi, Vol. 45, No. 4, April  1962.  [63]  J. Motard and B. Joseph. Wavelet Applications Publishers, Ma, U.S.A., 1994.  in Chemical Engineering.  [64]  K. Natarajan et al. On the Manipulation of the Slice Lip for the Control of Basis Weight. Technical report, PGRL 409, PAPRICAN, July 1988.  131  Kluwer Academic  [65]  R. W. Nelson. Approximate Theories of Filtration and Retention. Tappi Journal, Vol. 47, No. 12, pages 752-764, 1964.  [66]  Z. Nesic. Paper Machine Data Analysis Using Wavelets. Master's thesis, Dept. of Electrical Engineering, Univeristy of British Columbia, Canada, Feb. 1996.  [67]  D. Newcombe. Getting an Edge on Edge Waves. PIMA Magazine, April,  [68]  J. N . Newman. Marine  Hydrodynamics.  1994.  The MIT Press. Cambridge, Massachusette, U S A ,  1977. [69]  Bo Norman. Overview of the Physics of Forming. Cambridge Symposium,  [70]  S. Palavajjhala, R. L . Motard, and B. Joseph. Computational Aspects of Wavelets and Wavelet Transforms. In Wavelet Applications Academic Publishers,  [71]  1989.  in Chemical Engineering., Ed. Motard & Joseph,  Kluwer  Ma, U.S.A., pages 33—83, 1993.  J.R. Parker and B . A . Epton. Wet end Barring. CPPA Annual Conference, Montreal,  Canada,  Feb. 1976. [72]  R. S. Peterson. Improving Basis Weight Uniformity With Deckle Wave Control. Tappi Journal, July,  [73]  1992.  E . C. Pires et al. Computational model for Water Drainage in Fourdrinier Paper Machines. Tappi Journal, April 1988.  [74]  J.G. Proakis and D.G. Manolakis. Digital Signal Processing. Macmillan Publishing Company, 2nd Edition., 1992.  [75]  R. O. Ried and K. Kajiura. On the Damping of Gravity Waves over a Permeable Seabed. Trans. American Geophysists Union, Vol 38, 1957.  [76]  A . Rigopoulos, Y . Arkun, and F. Kayihan. Control relevant Disturbance Modeling Machine  of Paper  Full Profile Properties using Adaptive PC A. Control Systems 96, Halifax, Nova  Scotia, Canada, pp: 35-39, April 30-May2, 1996. [77]  O Rioul and M . Vetterli. Wavelets and Signal Processing. IEEE Signal Processing pp 14-38, Oct. 1991.  132  Magazine,  [78]  N . N . Sayegh and T. O. Gonzalez. Compressibilty of Fibre Mats during Drainage. Paprican Report, PPR 941, June,  [79]  1992.  D. E . Seborg, T. F. Edgar, and D.A. Mellichamp. Process Dynamics and Control. John Wiley & Sons, Inc. New York, N Y , U S A , 1989.  [80]  S. Skogestad and I Postlethwaite. Multivariable  Feedback Control, Analysis and Design. John  Wiley & Sons, Ltd. Baffines Lane, Chichester, West Sussex P019 1UD, England., 1996. [81]  G. A. Smook. Handbook Vancouver, B C , Canada,  [82]  Pulp  and Paper  Technologists.  Angus Wilde  Publication,  1989.  G. A . Smook. Handbook of Pulp and Paper Terminology. Angus Wilde Publication, Vancouver, B C , Canada,  [83]  for  1990.  R. M . Snyder, S. C. Chen, and Jr. Wilhelm. Adaptive Profile Control for Sheetmaking Process. 6th IFAC Conf. on Industrial pages 77-83,  Automation  in Paper, Rubber, Plastic and  Polymerization,  1986.  [84]  T. Soderstrom and P. Stoica. System Identification. Prentice Hall International (UK) Ltd., 1989.  [85]  G . Strang. Wavelet Transforms versus Fourier Transforms. Bulletin of the American matical Society, 28, 2, pp 288-305,  [86]  Mathe-  1993.  G Strang and T Nguyen. Wavelets and Filter Banks. Wellesley-Cambridge Press, Wellesley, Ma, U.S.A., 1996.  [87]  G. I. Taylor. Drainage of a Table Roll. Pulp and Paper Magazine 3,  [88]  G. I. Taylor. Drainage at a Table Roll and Foil. Pulp and Paper Magazine  of Canada, Vol.  1958.  M . L . Tyler and Morari M . Estimation of Cross-Directional Properties: Scanning vs. Stationary Sensors. AIChE,  [90]  Vol. 57, No.  1956.  59, 172, [89]  of Canada,  Vol. 41, No. 4, April  1995.  P.P. Vaidyabathan. Multirate Systems and Filter Banks. Prentice Hall, Englewood Cliffs, New Jersey 07632, 1993.  [91]  H . L . Van Trees. Detection, Estimation and Modulation  Theory. Wiley, New York., 1968.  [92]  R. Vyse, J King, M . Heaven, and S. Pantaleo. Consistency Profiling,  A New Technology for  CD Basis Weight Control. Pulp and Paper Canada, 97:9, pp:62—66, 1996. [93]  D. Wahren. Recent Highlights in Paper Technology. Tappi Journal, 69(3), pages 36-45, 1986.  [94]  J. Ward MacArthur. RMPCT: A New Robust Approach to Multivariable  Predective Control for  the Process Industries., pages 53—60. Control Systems 96, Halifax, Nova Scotia, Canada. April 30-May2, [95]  1996.  R. Wasler et al. Performance Analysis of Hydrofoils with Blades of Various Widths. Pulp and Paper Magazine of Canada,  [96]  71, T183, 1970.  L . Weiss. Wavelet and Wideband Correlation Processing. IEEE Signal Processing  Magazine,  pp 13-32, Jan, 1994. [97]  P. Wellstead and W.P. Heath. Two-Dimensional Control Systems: Application  to CD and MD  Control Problem. Control Systems 92 Conference, Whistler, B.C., Canada, pp: 39-43. 1992. [98]  W. Westermeyer. Modeling Flow in the Headbox Slice and on the Fourdrineir  Wire with Regard  to Cross Direction Basis Wieght Profile of the Paper Web. Das Papier 41 (11), pp 591-601, (English Translation by Institute of Paper Chemistry, Appleton, Wisconsin, U.S.A., Nov. 1987. [99]  R. K. Young. Wavelets Theory and Its Applications. Kluwer Academic Publishers, 2nd printing, Ma, U.S.A.,  [100]  1993.  R . H . Zhao and R.J. Kerekes. Pressure Distribution Between Forming Fabrics in Blade Gap Formers: Thin Blades. Paperi ja Puu, 78, pages 36—38, 1996.  134  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items