REAL TIME PAPER MACHINE DATA WAVELET ANALYSIS A N D CROSS D I R E C T I O N A L A C T U A T O R RESPONSE IDENTIFICATION By Ka Keung Kan B.A.Sc. (Electrical Engineering), The University of British Columbia, 2000 ., A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF A P P L I E D SCIENCE in T H E FACULTY OF G R A D U A T E STUDIES D E P A R T M E N T OF E L E C T R I C A L AND C O M P U T E R ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA April 2004 © Ka Keung Kan, 2004 Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. l2/04</ Zcro<f f Name of Author (please print) Degree: Date (dd/mm/yyyy) Year: H^CCL Department of ]^£cbrXcsJL The University of British Columbia Vancouver, BC Canada e QUAA GOXA^DIA^TPM (rwPAU4J?y(* Abstract The development of an online automated diagnostics tool for analyzing paper machine data is useful to mill operators when monitoring the performance of a cross directional (CD) control system. The use of Object linking and embedding for Process Control (OPC) interface as a medium for client-server based data acquisition provides an efficient, robust and accurate data transfer between a CD control system and a monitoring workstation. Therefore, a real time diagnostics toolbox has been developed to help paper mill companies reduce costs and save time on troubleshooting. The Real Time Wavelet Toolbox (known as Waveplot and developed over several years at UBC, [62]) provides color maps of raw, wavelet and residue profiles for both basis weight and moisture. Moreover, machine directional (MD) and cross directional (CD) control performance indices as well as the standard deviations for controllable and non-controllable CD components are calculated and displayed. Other important plots such as power spectral density and two-sigma plots are also included in the toolbox for further analysis. Waveplot also helps mill operators identify the high-resolution CD actuator response shape and CD mapping through bump tests. The identification is achieved by using a continuous wavelet transform with the chosen CD actuator response shape model as the mother wavelet. The control performance index reflects recent views on CD performance index in terms of CD controllability given the knowledge of the identified CD actuator response shape. Furthermore, the identified actuator response centers could help mill operators investigate any problems due to CD mapping misalignment and shrinkage, particularly at the paper edges. The displays and identification results are validated and demonstrated by collecting paper profiles through the OPC server from a multi-grade, operating paper machine. ii Table of Contents Abstract ii List of Tables vi List of Figures vii Acknowledgements xS Chapter 1 Introduction 1 1.1 Background 1 1.2 The Canfor Specialty Paper Machine 2 1.3 Motivation and Objectives 6 1.4 Previous Use of Wavelets on Paper Machine Data Analysis ... 7 1.5 Thesis Outline 8 Chapter 2 Wavelet Theory 9 2.1 From Fourier to Wavelet Transform 9 2.2 Wavelet Bases 12 2.2.1 Admissibility 12 2.2.2 Orthogonality 13 2.2.3 Compact Support 13 2.2.4 Symmetry 13 2.2.5 Smoothness 14 2.2.6 Number of Vanishing Moments 14 2.3 The Discrete Wavelet Transform 14 2.3.1 Multi-resolution Analysis 15 2.3.2 Filter Banks 16 2.4 De-noising 2.4.1 18 Thresholding Methods iii 20 2.4.2 Threshold Value Selection 20 Chapter 3 Cross Directional Actuator Response Identification 23 3.1 Introduction 23 3.2 CD Actuator Response Identification Using Wavelets 25 3.2.1 Problem Formulation and Assumptions 25 3.2.2 CD Actuator Response Model 27 3.2.3 Correlation Analysis With Wavelets 31 3.3 CD Actuator Response Identification Procedure 32 3.3.1 Extraction of Bump Test Profiles from Industrial Data . 33 3.3.2 Peak Search 3.3.3 Correlation Analysis using Continuous Wavelet Transform 39 36 3.4 Parameter Estimation Results and Discussions 44 3.4.1 Identification Results 44 3.4.2 Response Model Comparisons 45 3.4.3 Response Width 50 3.4.4 Significance of Identified Response Shape 51 3.5 Summary 56 Chapter 4 Real Time Paper Machine Data Wavelet Analysis 58 4.1 Introduction 58 4.2 Overview of Object Linking and Embedding for Process Control (OPC) 59 4.3 Waveplot for Real Time Paper Machine Data Analysis Tool ... 60 4.3.1 Real Time Data Access Methodology 60 4.3.2 Overview of Main Features of Waveplot 66 4.3.2.1 Color Maps Display 67 4.3.2.2 MD and CD Performance Indices 4.3.2.3 Standard Deviations iv . 72 74 4.3.3 4.3.4 Special Features 76 4.3.3.1 Time Stamps 76 4.3.3.2 Menu Bar 76 4.3.3.3 Crosshair 79 4.3.3.4 Turn-up Marks and Machine Status 79 4.3.3.5 Real Time Display on Workstations 82 4.3.3.6 Historical Storage 82 A Real Time Diagnostics Example 83 4.4 Summary r- Chapter 5 Conclusions and Future Work 89 90 5.1 Conclusions 90 5.2 Future Work 92 Appendix A 93 A.l Proof of Admissibility of Mother Wavelets . 93 A.l.l ip(x) = e~ ^cos{-irx) 93 A.l.2 ip{x) = e~ cos(-nx) 94 a ax2 A.2 Computation of MD Performance Index 95 A.3 Computation of Power Spectral Densities of CD Profiles 99 Bibliography 101 v List of Tables 3.1 Actuator Profile Peak Locations and Corresponding Wavelet CD Profile Peak Locations 37 3.2 Squared Estimation Error for Identified CD Response Shapes with CWT. . 50 3.3 Parameter w for Identified CD Response Shapes with CWT 51 3.4 Sum of Squared Error of Estimation Error with Different Response Lengths. 52 3.5 Important Data for Frequency and Phase Plots for Averaged Identified Actuator Responses 55 4.1 Required Paper Machine Data and the Number of Tags vi 64 List of Figures 1.1 A Simplified Diagram of a Paper Machine 2 1.2 Sensor Path 3 1.3 The Canfor Specialty Paper Machine 4 1.4 The Canfor Specialty Paper Machine (Close Up) 4 1.5 MH Headbox-C CIV Chamber 5 2.1 Time-frequency Plane: FT, STFT and WT 11 2.2 Nested Vector Space Spanned by Wavelet and Scaling Functions 16 2.3 Two-stage Two-Band Wavelet Decomposition Tree 17 2.4 Two-stage Two-Band Wavelet Synthesis Tree 18 2.5 Schematic Representation of a Noise Reduction System Using Thresholding. . 19 2.6 Hard Thresholding 21 2.7 Soft Thresholding 21 3.1 Arrangement of a Bump Test 24 3.2 CD Actuator Response Model #1: r(x) = e~ \ \cos{-Kx) 29 3.3 CD Actuator Response Model #2: r{x) = e~ "'COS(TTX) 30 3.4 Extraction of Bump Test Profiles 33 3.5 Differential Wavelet CD Profiles of a Bump Test (1-4) 34 3.6 Differential Wavelet CD Profiles of a Bump Test (5-8) 34 3.7 Differential Wavelet CD Profiles of a Bump Test (9-12) 35 3.8 Differential Wavelet CD Profiles of a Bump Test (13-16) 35 a x ax 3.9 Peaks of Differential Actuator Setpoint and its Corresponding Wavelet CD Profile of a Bump Test for Identification 38 3.10 Actuator-databox Alignment Line 39 3.11 Normalized Wavelet Bump Response #6 41 3.12 Color Map for Wavelet Coefficients of Response #6 with Scale 1-100 and Response Model #1 41 vii 3.13 Color Map for Wavelet Coefficients of Response #6 with Scale 20-40 and Response Model #1 42 3.14 Color Map for Wavelet Coefficients of Response #6 with Scale 1-100 and Response Model #2 42 3.15 Color Map for Wavelet Coefficients of Response #6 with Scale 20-40 and Response Model #2 43 3.16 Identified CD Response Shapes 43 3.17 Identification of CD Responses with Response Model #1 using Continuous Wavelet Transform. Dotted Line:True Response; Solid Line:Identified Response 45 3.18 Identification of CD Responses with Response Model #2 using Continuous Wavelet Transform. Dotted Line:True Response; Solid Line:Identified Response 46 3.19 Average Identified Response for Response Model #1. Dotted Lines:Identified Responses for 9 Bumped Actuators; Solid Line:Average Identified Response. 47 3.20 Average Identified Response for Response Model #2. Dotted Lines:Identified Responses for 9 Bumped Actuators; Solid Line:Average Identified Response. 48 3.21 Identification Results Using Continuous Wavelet Transform. Dotted Line:True Response; Solid Line:Estimated Response 49 3.22 Frequency Response for the Identified Actuator Responses Using Model #1. Blue Line: Average Response; Red Lines: All Bumped Actuator Responses. 53 3.23 Frequency and Phase Response for the Identified Actuator Responses Using Model #2. Blue Line: Average Response; Red Lines: All Bumped Actuator Responses 54 3.24 Power Spectral Densities for Basis Weight 56 4.1 OPC Client/Server Relationship 61 4.2 Real Time Data Access Diagram 62 4.3 Diagram for Data Transfer Between OPC Server and Matlab Client 63 viii 4.4 Flowchart for the Data Retrieval Process 65 4.5 Waveplot GUI Window 68 4.6 Raw Image for Basis Weight Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) 69 4.7 Raw Image for Moisture Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) 70 4.8 Wavelet Image for Basis Weight Collected on November 12, 2003 at 23.38pm. (Scan Vs Databox) 70 4.9 Wavelet Image for Moisture Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) 71 4.10 Residue Image for Basis Weight Collected on December 2, 2003 at 23:38pm. (Scan Vs Databox) 71 4.11 Residue Image for Moisture Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) 72 4.12 Wavelet Decompostion Tree and the Wavelength for each Node 74 4.13 Wavelet Color Map and MD and CD Performance Indices 75 4.14 Standard Deviations 77 4.15 Time Stamps in Waveplot GUI Shown at the Edit Boxes 78 4.16 Basis Weight Actuator Setpoints 80 4.17 Crosshair 81 4.18 Mill Workstation 82 4.19 Wavelet Color Map for Basis Weight Collected on November 10, 2003 at 08:10am. (Scan Vs Databox) 84 4.20 Wavelet Color Map for Basis Weight Collected on December 2, 2003 at 17:23pm. (Scan Vs Databox) 84 4.21 Wavelet Color Map for Moisture Collected on November 10, 2003 at 08:10am. (Scan Vs Databox) 85 ix 4.22 Wavelet Color Map for Moisture Collected on December 2, 2003 at 17:23pm. (Scan Vs Databox) 85 4.23 MD Performance Indices for Basis Weight Collected on November 10 and December 2, 2003 86 4.24 MD Performance Indices for Moisture Collected on November 10 and December 2,2003 87 4.25 CD Performance Indices for Basis Weight Collected on November 10 and December 2, 2003 87 4.26 CD Performance Indices for Moisture Collected on November 10 and December 2, 2003 88 x Acknowledgements I would like to take this opportunity to thank my supervisor Dr. Guy Dumont for his guidance and support. I would like to thank John Ball from Canfor pulp and paper mill for all his efforts in making the installation of Waveplot successful. I truly respect him as a knowledgeable process control engineer professional with endless patience in helping me understand the complex paper machine process for making my stay in Canfor enjoyable. Special thanks to Zoran Nesic for helping me at the beginning of my masters project. I would also like to thank the administration staff at the UBC Pulp and Paper Centre, especially Ms. Lisa Hudson, for helping me make travel arrangements for my visits to Prince George. I appreciate Canfor IT Support Help Desk staff especially Stace Clark for helping me on setting up workstations and remote connections. My appreciation also goes to the professors and students in the Process Control Group at UBC for giving out seminars and Dr. Michael Davies for organizing the group meetings. My heartfelt thanks to my best friends Dierdre Tse and Meoni Poon for their precious friendships in terms of emotional support, sharing good and bad times as well as encouragement that have enlightened me to carry on till the last minute. Thanks to Stephen Shapiro for his encouragement to help me through the hardship of graduate school. Last but not the least, I especially thank my parents, my sister and my brother for their patience during all those long years believing in me and allowing me the time to get things right. The achievement for finishing my thesis would not be possible without their endless support and love. xi Chapter 1 Introduction 1.1 Background A paper machine is used to transform a mixture of water and fibres into a sheet of paper efficiently. Figure (1.1) illustrates the main parts of a paper machine. Thick stock consisting of approximately 0.5% of fibres and 99.5% of water is pumped into the headbox of the paper machine. The purpose of the headbox is to distribute the stock as uniformly as possible across the machine. There are also actuators in the headbox distributed across the machine to maintain uniformity of the stock. The slurry then travels along the wire mesh belt at the speed of 50km/h. The wet pulp is dewatered on its way to the mid-section of the machine. The pulp with approximately 20% of solid is then heated by steam boxes and pressed by the counter-rotating rollers. A sheet of paper leaving the dryer section after a series of dewatering and pressing processes contains 5 to 9% of water content by weight. Finally, the sheet of paper enters the calendar stack where the thickness of the paper is controlled and the final product is wound at the reel. For a complete description for this complex process, refer to [35, 20, 21]. There are three main paper properties in paper production namely basis weight in g/m , 2 moisture expressed as percentage of water content of the sheet and the thickness of the sheet of paper expressed commonly in fxm. These properties are measured by a scanner sensor which traverses back and forth across the paper machine. As the direction of travel of the scanner is perpendicular to the direction of the papermaking process, the trajectory of the measurements produces a zigzag pattern as shown in Figure (1.2). Consequently, the measurements contain a combination of Machine Direction (MD) and Cross Direction 1 Chapter 1. Introduction 2 Figure 1.1: A Simplified Diagram of a Paper Machine. (CD) components. This 2-dimensional control problem has been the subject of rigorous research since the past decade because the decoupling of two orthogonal components and the control in two directions simultaneously are not easy tasks. 1.2 The Canfor Specialty Paper Machine The paper machine that this thesis deals with is located at the Canfor pulp and paper mill in Prince George. This sack kraft paper machine was built in 1965. In the mid 90's, the company experienced a declining market for their traditional quality unbleached product. In order to remain competitive, the mill decided to rebuild the machine to increase production and to diversify the product line. Nowadays, the paper machine can produce more profitable paper products with different grades of unbleached, bleached and colored paper [19]. The Canfor Specialty Paper Machine now has high consistency refining, electric drives, 5 new drying cans, new fan pump, new screens and approach piping as well as an upgraded Chapter 1. Introduction 3 Machine Direction Q o k*i k k-1 Figure 1.2: Sensor Path. ABB Smart Platform scanner [19]. The basis weight and moisture scanners have 600 databoxes at 9.85mm width resolution. The paper machine is also equipped with a new Mitsubishi Heavy Industries CIV dilution headbox with Beloit Concept IV design which was installed in the spring of 2001. Figure (1.3) and Figure (1.4) show the pictures of the headbox that has been installed in the Canfor Specialty Paper Machine in Prince George. Figure (1.5) shows the headbox with consistency profiling. The headbox is a hydraulic type unit which allows stock flow from the fan pump comes through the piping system and into the head rounder (not shown in the figure) which is located behind the tube bank. The multiple tube bank which is located in front of the tube bank provides uniform CD distribution of the stock into the headbox chamber [31]. The uniformity of the stock flow is heavily dependent on the turbulence of the stock flow. The turbulence must be reduced to an acceptable level to achieve uniform dispersion of the stock in the tube bank. Control of turbulence intensity of the stock flow is achieved by using several mixing devices in the headbox. The flow passes through the tube bank and then into the headbox chamber where further macro-turbulence occurs. The flow sheets Chapter 1. Introduction Figure 1.4: The Canfor Specialty Paper Machine (Close Up). 1 Chapter 1. Introduction Figure 1.5: MH Headbox-C CIV Chamber. 5 Chapter 1. Introduction 6 which are made of lexan (polycarbonate) indicated in Figure (1.5) are therefore installed as converging channels within the headbox chamber to reduce unwanted macro-turbulence so that a uniform dispersion of the stock can still be maintained [31]. The flow sheets have to be designed with correct thickness and length to achieve desirable turbulence of the stock flow otherwise excessive stock mixing in the cross direction could happen. The headbox is also equipped with Consistency Profiling system shown in Figure (1.5). This system provides narrow-band profile control by injection of dilution water near the tube bank. This causes a lowering of consistency in the adjacent headbox tube and reduction of basis weight at the particular cross direction position. The Consistency Profiling system eliminates the need for mechanical slice-profiling adjustment [31]. The headbox consists of 78 actuated water injection tubes with spacing of 82mm for controlling basis weight. Furthermore, a steam box with 40 CD zones is available for moisture control. 1.3 Motivation and Objectives Mill companies have always looked for ways to save time and costs in troubleshooting. A control problem could take days for mill engineers to solve since the cause of the problem could be unknown in which they could have no idea on how to narrow down the problem. Furthermore, control problems are often not discovered immediately, leading to further loss of costs in labour on troubleshooting while additional costs have already been lost in producing paper with unacceptable quality. Therefore, there is a need to analyze paper machine properties in real time such that mill operators can be notified immediately if control problems occur. The main objective of this thesis is to develop an online automated diagnostics tool using wavelets for analyzing paper machine data. The use of wavelets not only gives a better representation of paper machine data, but also provides mill operators a great visual tool as quality images of two-dimensional data (CDxMD) to examine events such as sudden Chapter 1. Introduction 7 change in paper quality or streaks. The other objective would be to identify the CD actuator response for the particular type of headbox mentioned above. The identified CD actuator response shape is useful for CD variations separation so that control performance can be further improved. 1.4 Previous Use of Wavelets on Paper Machine Data Analysis The use of wavelets on analyzing paper machine data started almost a decade ago. Nesic [62, 63, 61] first proposed using discrete wavelet transforms for paper machine data handling. They showed that with proper choice of wavelets, the performance of profile estimation, visualization and storage are enhanced comparing, with standard approaches such as exponential filtering. Waveplot was alsofirstdeveloped to analyze historical paper machine data to display color maps of raw, wavelet and residue profiles of basis weight and moisture for diagnostics. By taking advantage of high computational speed of the algorithm, Waveplot has been implemented to allow real time analysis and this is also the main theme of the thesis. Jiao [34] furthered the use of wavelets in trim-loss optimization to reduce the cost of raw fibre and energy in the papermaking industry. Besides using wavelet analysis, she also proposed using wavelet packets to enhance the performance of the signal decomposition. Moreover, she applied wavelets and wavelet packets on CD variation separation. The separation into CD controllable and non-controllable variations is necessary for performance monitoring. Gorinevsky and Gheorghe [23] as well as Ghofraniha [33] used continuous wavelet transform for system identification. The CD actuator response shape can be identified by correlating the CD profiles collected from bump tests with an appropriate CD actuator response model as the mother wavelet. The identified CD actuator shape is very useful in CD variation Chapter 1. Introduction 8 separation and calculating CD performance indices which will be explained in the thesis later. Recently, Dumont, Ball, Emmond and Guillemette [19] proposed using wavelets for MD/CD separation instead of using exponential filtering of raw scanned data obtained from basis weight and moisture gauge scanners for CD profile estimation. The separation is achieved by using the Discrete Wavelet Transform (DWT). The DWT can be implemented as a filter bank to decompose a signal into approximations and details in different resolution levels. The CD profile can be estimated by finding a wavelet level that separates CD variations into controllable and non-controllable components. The authors mentioned above made significant contributions in applying wavelet theory to paper machine analysis. This thesis continues the use of wavelets to develop an online automated diagnostics tool for immediate CD performance monitoring with new concepts introduced to reflect the most recent development in this area. 1.5 Thesis Outline A brief introduction to wavelet theory is given in Chapter 2. Chapter 3 focuses on the CD actuator response identification. In Chapter 4, a detailed description on the real time diagnostics tool known as Waveplot is given with a brief overview on Object linking and embedding for Process Control (OPC). The thesis concludes by stating research results and future work in Chapter 5. Chapter 2 Wavelet Theory The wavelet theory was first developed by Morlet and Grossman [3] in the early 1980's. In mid 1980's, Meyer and Lemarie [42] suggested a new orthogonal wavelet expansion for a signal function. However, it was Daubechies [28] who solidified the wavelet theory by proposing a method to construct wavelets withfinitesupport and high degree of regularities in 1988. The theory was then adopted in different areas within the applied mathematics communities namely signal processing, statistics and numerical analysis. Many researchers have published their work on wavelets. Daubechies [28, 27, 30, 29], Meyer [42, 59], Cohen [2] and Wickerhauser [40] have significantly contributed to the development of wavelet theories. In this chapter, an introduction to the wavelet theory including discrete wavelet transform, wavelet bases and filter banks as well as applications such as data de-noising are presented. 2.1 From Fourier to Wavelet Transform The Fourier Transform theorem was first developed by Fourier in 1807 which is still widely used in many disciplines of science because it has the ability to analyze the frequency components of a linear time-invariant signal. The sinusoidal waves e ZUJt act as eigenvectors for decomposing a signal into different frequencies. A typical Fourier integral for a stationary signal f(t) is given by: '+0O (2.1) 9 Chapter 2. Wavelet Theory 10 Equation (2.1) measures the amplitude or intensity that the signal fit) oscillates at a certain frequency u>. However, the main drawback for the Fourier Transform is that it is not suitable to analyze signals with transient characteristics. Therefore, in 1946, Gabor [8] proposed the windowed Fourier Transform or Short-Time Fourier Transform (STFT) by decomposing a '-non-stationary signal over a time-frequency plane which is defined as following: (2.2) where g(t — r) is the window function centered at a time r. The STFT takes one step forward to accommodate time-varying signals by evaluating their frequency spectrum at certain time t. The energy of f(t) is therefore localized over a time interval. However, the window size for the STFT is fixed resulting in only one resolution in both time and frequency. Wavelet basis functions are therefore created to represent a signal with different resolutions and time localization capability. Similar to the STFT, Wavelet Transform can measure time-frequency variations of spectral components, but in different resolutions. The basis functions of Wavelet Transform for multi-scaled signal representation are created by scaling and translating a mother wavelet basis function ip(t) which'is called mother wavelet. The wavelet basis function is defined as: (2.3) where a is the scaling factor with a ^ 0, and b is the translating factor. The scaling factor a can be viewed as an indication of the frequency resolution. In fact, the time resolution is inversely proportional to the frequency resolution. When the scale gets smaller, more translations are required to cover the entire signal. On the other hand, larger scale implies less translations. Chapter 2. Wavelet Theory 11 The time and frequency resolutions are defined as: n~tw(t)\dt At = /£V|* (w)|du; (2.4) 2 Au (2.5) where ^(u) is the Fourier Transform of the mother wavelet tp(t). The limitation for how refined can the resolution be is determined by the uncertainty principle or Heisenberg inequality [41]: . . 1 2 AtAu > _ (2.6) Figure (2.1) shows time-frequency planes for Fourier, STFT and Wavelet Transform to represent a signal in frequency domain. From the figure, At of the Fourier time-frequency plane is infinite which implies there is no time localizaton. For STFT, At is finite but it is fixed which implies there is only one fixed resolution. The time-frequency plane for wavelet transform shows varying Ato and At meaning that it can represent a signal with more than one resolution. Figure 2.1: Time-frequency Plane: FT, STFT and WT. Chapter 2. Wavelet Theory 2.2 12 Wavelet Bases The wavelet basis functions defined in Equation (2.3) have important properties which become foundations in developing subsequent theories such as multi-resolution analysis and filter banks, as well as wavelet packets. For details of the wavelet basis functions properties, see [4]. 2.2.1 Admissibility A mother wavelet ip(t) has to be oscillatory in nature so that it can capture frequency components of a signal. It should also be a function with fast decay for better time localization. These two properties can be ensured by checking the admissibility condition of ip(t) [24]: However, this condition is sufficient but not necessary [55]. In order for Equation (2.7) to hold, the following has to be satisfied: +oo ip(t)dt = 0 (2.8) The wavelet coefficients can be computed by correlating the signal /(£) with the dilated and translated version of mother wavelet function tp(t) as shown in the following: +oo / f{tW , {t)dt=<f{t)^ {t)>, j k k (2.9) •oo where < . > is the inner product operator and ipj^it) are the wavelet basis functions defined in Equation (2.3). Chapter 2. Wavelet Theory 2.2.2 13 Orthogonality A wavelet basis function defined in Equation (2.3) is orthogonal if the set ip(t — k)\k € Z forms an orthonormal basis: < ^,k(t),il} (t) >= 8(j -jo)5{k- k ) jolko 0 (2.10) Filter banks are constructed based on the orthogonal property such that perfect reconstruction is guaranteed. The orthogonal property is also important in terms of numerical computations and error minimizations. The fast wavelet transform, which decomposes successively the approximation of f(t) in subspace Vj into a coarser approximation Vj+i plus the wavelet coefficients of f(t) in subspace Wj+i or recovers f(t) in subspace Vj from V j+i and Wj+i subspaces, is a unitary transformation with the help of the orthogonal property. Therefore, the transform will not increase the error in the initial data. For details of the fast wavelet transform, refer to Section 2.3. 2.2.3 Compact Support The compact support of scaling functions (see Equation (2.13)) and wavelet functions guarantee finite impulse response filters resulting in finite fast wavelet transform. This property also ensures the perfect reconstruction of the filter banks. 2.2.4 Symmetry Symmetry of the wavelet filters is essential to have linear phase otherwise phase distortion could result. Chapter 2. Wavelet Theory 2.2.5 14 Smoothness A function is said to have M degrees of smoothness if the M derivative of the function th is continuous at all points. This property is significant for de-noising and compression applications. The reconstructed signal after de-noising and compression processes will be much smoother if the wavelet function has a high degree of smoothness. For more details, refer to [46]. 2.2.6 N u m b e r of V a n i s h i n g M o m e n t s The number of vanishing moments k for a sequence {gk} is defined as the number of moments of wavelet that are zeros: (2.11) fc The higher the number of vanishing moments, the greater is the ability in detecting singularities and giving smoother wavelet [1, 7, 54]. 2.3 The Discrete Wavelet Transform The wavelet basis function in Equation (2.3) is considered to be continuous in time variable t which implies the variables a and b to be continuous as well. The wavelet transform can be implemented in a discrete manner where variables of time, translation and dilation in the wavelet basis function are sampled with appropriate step sizes. Equation (2.12) shows the discrete wavelet basis function. i(j , (t) = 2 ^i){2 t-k), j j j k (2.12) where j is the scaling factor with j ^ 0, and k is the translating factor; j, k € Z. For analyzing the paper machine data, the discrete wavelet transform (DWT) and inverse Chapter 2. Wavelet Theory 15 discrete wavelet transform (IDWT) are used. The DWT and IDWT are closely related to the multi-resolution analysis since they are implemented with the use of Mallat's pyramid algorithm [53] or filter banks which are originated from the concepts of multi-resolution formulation. 2.3.1 Multi-resolution Analysis The concept of multi-resolution is based on the analysis of a signal in different resolutions. The paper machine property data collected from the scanner can be approximated in different levels of resolutions to preserve important features so that a smoother signal can be obtained. Any function f(t) in the square integrable space £ {1Z) can be expressed as a 2 sum of scaling functions which represent coarser information and wavelet functions which represent finer resolution. The scaling function with scaling factor j and translating factor k is: (2.13) A set of scaling functions can be generated to span over k in the subspace V} [7]: Vj = Spon {ip (pt)} k k (2.14) Figure (2.2) illustrates the scaling function vector spaces Vj and wavelet vector spaces Wj. In multi-resolution analysis, any square integrable function can be decomposed using scaling and wavelet functions which are depicted graphically in Figure (2.2) with the following properties [7, 34]: 1. a: V = V ®W x 0 0 b: V = V ®W 2 1 c: W ±W ± r 0 l V 0 Chapter 2. Wavelet Theory 2. V, C V j+1 ; 3. f{t) € V, ^ jeZ /(2t) € V 4. f(t) eV 4=> 6. i + 1 f(t+l)£V 0 5- n 16 0 z ^ - { 0 } ; U ^ = £ (^) 2 ; e ; S V such that 0 — k)} z forms an orthonormal basis for V ne 0 Figure 2.2: Nested Vector Space Spanned by Wavelet and Scaling Functions. 2.3.2 Filter Banks Filter bank algorithms are based on the multi-resolution analysis to compute the orthogonal wavelet coefficients of a signal measured at afinerresolution. The algorithm also allows the reconstruction of wavelet coefficients at V,_i from coefficients in subspaces Vj and Wj. For wavelet decomposition, a signal s(n) with approximation coefficients Cj_i is split into approximation coefficients Cj and detail coefficients dj by passing low passfilterh (—n) and 0 high pass filter h\(—n) with a downsampling-by-2 operator implemented after each FIR filter. The signal can be further decomposed by splitting the approximation coefficients Chapter 2. Wavelet Theory 17 into Cj+i and d i and so on. Therefore, the algorithm can be viewed as an iterative j+ process. For reconstruction, the coefficients Cj_i can be recovered from Cj and dj. Both Cj and dj are first upsampled by 2 and then they are further processed by passing another set of low pass filter / (n) and high pass filter /i(n) respectively. c,-_i is finally computed by 0 summing two sets of coefficients coming from each filter output. Figure (2.3) and Figure (2.4) show the block diagram for two-stage two-band DWT and IDWT adopted from the Mallat's pyramid algorithm respectively. There are conditions in choosing the coefficients of the low pass and high passfiltersto achieve perfect reconstruction. The two main conditions for perfect reconstruction are: Alias Cancellation: F (z)H (z) + Fi(z)Hi(z) = 2z~ ; I odd (2.15) l 0 0 No Distortion: F (z)H (-z) 0 0 + F {z)H {-z) 1 1 (2.16) = 0, where H (z), H\(z), F (z) and Fi(z) are z-transforms of h (n), h\(n), /o(n) and /i(n) 0 0 0 W h,(-n) <B w. h (-n) V. M-n) M-n) <5> Figure 2.3: Two-stage Two-Band Wavelet Decomposition Tree. Chapter 2. Wavelet Theory 18 W i Figure 2.4: Two-stage Two-Band Wavelet Synthesis Tree, respectively. For details see [25]. 2.4 De-noising Noise reduction is one of the most important applications of wavelets. De-noising using wavelet transformation refers to the recovery or reconstruction of a signal from noisy data by shrinking the empirical wavelet coefficients towards zero followed by inverse wavelet transform [1]. Basically, the wavelet regression estimator [61] is implemented by the following steps: 1. Compute the discrete wavelet transform of a finite signal sequence {y} as: w = Wy, (2.17) where W represents the wavelet transform. 2. Modify the wavelet coefficients according to appropriate criterion which could be highly non-linear. The new set of wavelet coefficients is given by w*. Chapter 2. Wavelet Theory 19 3. Compute an estimate of y by applying inverse discrete wavelet transform to the modified wavelet coefficients w*: y = W W* (2.18) T Figure (2.5) illustrates the block diagram of a noise reduction system using wavelet transformation. ) 1 k Discrete Wavelet Transform w Wavelet Coefficients Thresholding w* Inverse Discrete Wavelet Transform •x(t) Figure 2.5: Schematic Representation of a Noise Reduction System Using Thresholding. Different methods for de-noising have been proposed and implemented by Donoho and Johnstone at Stanford University. For details, see [10, 11, 12, 13, 14]. They suggested that a signal be modeled as: Vi = fi + o-ei, i= l,...,n (2.19) where yi is the noised measurements, fa is the noise-free version of y a is the noise level it and €i is the zero-mean, unit-variance Gaussian white noise. In order to compute the de-noised estimate y, non-linear thresholding is applied to the wavelet coefficients. A real signal is usually represented by a small number of large wavelet coefficients. The advantage of thresholding is that since noise content is represented by small wavelet coefficients, the noise from the signal can be greatly reduced by cutting them off based on suitable threshold selection criterion. Chapter 2. Wavelet Theory 2.4.1 20 Thresholding Methods There are two thresholding methods namely hard and soft thresholding. Hard thresholding keeps all the coefficients whose values are above a threshold value A while it sets the remainder where the values are below A to zero. Mathematically, it is expressed as: x if x > A WhardiX) 0 otherwise (2.20) For soft thresholding, it sets the coefficients which are smaller than the threshold to zero and shrink the others by the threshold value. It is given by: x — A if x > A Wsoftix) = 0 if |x| < A (2.21) x + A if x < — A Figure (2.6) and Figure (2.7) illustrate the plots of Equation(2.20) and (2.21) respectively. Equation (2.20) and (2.21) suggest that these methods are highly non-linear. Due to the sharp discontinuities that hard thresholding imposes, the de-noised signal could contain irregularities so it does not look smooth. Soft thresholding could lead to a smoother noisefree signal, but it has larger bias. 2.4.2 Threshold Value Selection The choice of the threshold A is critical since it could affect the quality of the estimates. Donoho and Johnstone in [10] and [13] have suggested the Universal Threshold or VisuShrink method. The threshold calculation is given by: i = a«y/21og(n), (2.22) Chapter 2. Wavelet Theory 21 Threshoided Wavelet Coefficients Figure 2.6: Hard Thresholding. Threshoided Wavelet Coefficients Figure 2.7: Soft Thresholding. Chapter 2. Wavelet Theory 22 where n is the sample size of the signal which is zero-padded to the nearest power of 2 (2 ) and a is the noise-level estimate. There are three main different ways to compute the N noise-level estimate. The standard deviation method, abbreviated as "std", computes the standard deviation of the wavelet coefficients as the estimate. Another method, Median Absolute Deviation (MAD), calculates the estimate based on the following: b = Median{\w \)/0.6745, f (2.23) where Wf is a set of wavelet coefficients of the signal. This estimator produces smooth and robust results. Hence, the estimates are visually appealing so this method is also called Visushrink. The drawback for this method is that the same threshold is applied to all levels of wavelet transformation. This could cause the data to be over-smoothed. Therefore, another method is suggested and is called MultiMedian Absolute Deviation (MultiMAD). This method is an adaptive threshold selection procedure meaning that different thresholds are computed at different resolution levels of the wavelet transformation. Chapter 3 Cross Directional Actuator Response Identification 3.1 Introduction In cross directional (CD) control systems, an array of actuators is distributed across the paper web to spatially influence the paper properties such as basis weight and moisture. The quality of CD control can not be ensured without a full understanding of CD spatial actuator response. For this reason, CD identification plays a major role in CD control problems since discrepancies in identifying the CD actuator response shape could lead to producing paper with low quality or even damaging the paper machine. In addition, the alignment between the actuators and the measurements is critical for CD control systems since mapping errors could lead to system instability. The correct alignment relies on a sophisticated CD identification scheme to identify the geometrical parameters of a paper machine. CD identification becomes important to analyze the shrinkage behaviour at the paper edges and the wandering of the paper web. In fact, paper shrinkage is a non-linear phenomenon so that its accurate identification is needed. In the paper-making industry, the identification is carried out by conducting bump tests meaning that several actuators are manually stepped to observe the process response. Alternatively, actuators can also be probed with a sequence of setpoint changes with random duration known as pseudo-random binary sequence (PRBS) [49]. The probing sequences of multiple CD actuators must be statistically independent random sequences. In this thesis, bump tests have been used to observe the responses of bumped actuators. Figure (3.1) 23 Chapter 3. Cross Directional Actuator Response Identification 24 Figure 3.1: Arrangement of a Bump Test. shows graphically how a bump test is arranged [23]. In the figure, actuators located at the headbox are bumped and the response to the bumped actuators is observed based on the data collected from the scanner. Bump tests are usually conducted when the process is in open-loop. Most recent methods for identifying temporal and spatial responses use the interactive matrix approach to model the sheet-forming processes. The parameters are then identified by using traditional parametric identification techniques such as least squares and maximum likelihood as well as modern time-series analysis [50, 56, 15, 44]. Recently, Chen [51] proposed a new method to model the sheet-forming processes such that simultaneous identification for both temporal (MD) and spatial (CD) responses can be achieved. Although this method allows identification of the responses of all actuators across the paper machine, its recursive nature and the use of large number of parameters result in a huge need for computational resources. As the actuator technology nowadays allows enhanced Chapter 3. Cross Directional Actuator Response Identihcation 25 resolution by increasing the number of actuators and measurements across the paper web, these methods would become time consuming and computationally expensive. Moreover, the process noise corrupting the identification data decreases the integrity and accuracy of the identification. Therefore, in this thesis, an identification method is chosen based on using relatively simpler model with fewer parameters so that the algorithm is computationally efficient and easy to implement. Furthermore, with the help of wavelets, identified CD response shape can be obtained with acceptable computational loss and the noise on the identification results can also be greatly decreased. This chapter discusses the CD identification using wavelets. The identification problem is first addressed and then the approach using continuous wavelet transform is explained. The algorithm is based on methods proposed by Gorinevsky [23] and Ghofraniha [33] but it is applied on paper machine equipped with dilutionflowactuators. The identification results are then illustrated based on the industrial data. Results such as the identified actuator response shape could have a significant impact on CD controllability which is assessed by performance indices. The calculation of the performance indices will be discussed in Chapter 4. 3.2 3.2.1 CD Actuator Response Identification Using Wavelets Problem Formulation and Assumptions The CD measurements are collected using a scanner sampled at regular intervals across the sheet of paper. The data points collected from the scanner traversing from one edge of the paper web to the other end are described as a scan. In mathematical form, the CD profile of m measurements at scan t is described by a vector p(t) € ~R . The dynamics for m the process of the paper machine is given by [23]: p(t) = Gp(z- )u(t), 1 7 (3.1) Chapter 3. Cross Directional Actuator Response Identification 26 where u(t) € lZ is the actuator profile at scan t, z~ is the delay operator, g(z~ ) is a n x l single-input-single-output transfer function describing the time response behaviour of the paper machine, G is the CD interaction matrix with dimension m-by-n and 7 is the overall gain of the process. It is assumed that the process can be separated into dynamical (MD) and spatial (CD) responses where g(z~ ) represents the dynamical response which is modeled as first-order l lag with time delay and G represents the latter [48, 50, 56, 15]. Equation (3.1) can be interpreted as having two responses connected in a cascade form. It is also assumed that the dynamical and spatial models are parametric. The ultimate goal for identification is to identify the parameters defined in these two models so that the CD response can be computed using Equation (3.1). Consequently, the parametric models can be set up with 6T and 8cD being the model parameter vectors for dynamical and spatial responses respectively with: giz-^^giz-^Or), (3.2) G = G(6CD)- (3.3) In this thesis, the time response is assumed to have reached steady state which means the parameters in 9T are constants. This assumption is justifiable since MD process dynamics can always be assumed to be fast compared with the measurement sampling time [9]. As a result, the overall identification problem can be reduced to a CD identification problem. This assumption is widely accepted by the paper-making industry. For CD spatial responses, it is assumed that the process dynamics are linear which implies that the system response can be evaluated as the superposition of the actuator responses [9]. Moreover, it is also assumed that all actuators have the same spatial response. However, these two assumptions are only valid to a certain degree since non-linearity such as paper shrinkage Chapter 3. Cross Directional Actuator Response IdentiRcation 27 still exists in the papermaking process. 3.2.2 C D Actuator Response M o d e l The objective of CD identification is to identify the CD actuator response shape and to assess the CD mapping between the actuators and the measurements. The model for CD response isfirstformulated followed by describing the identification algorithms using continuous wavelet transform. The interaction matrix G describes the spatial response of all actuators. With the assumption that each actuator has the same response shape and the response shape is symmetrical, the matrix is in a form of banded diagonal matrix [37]: G 0 Po Pi ••• Pk 0 Pi Po Pi '•• '• '• Pi "•• '•• '• 0 Pk '• '•• '•• '• Pk 0 •. •. ••. • •• • 0 ••• 0 Pi Pi p k (3.4) Po ••• Pi pi po In a simpler form, the entries of matrix G can be written as: G J K = r](d k-Cf,ecD), s (3.5) .where n{x; OCD) is the CD response function, d is the databox width, Cj is the CD response s center in databox coordinate for the actuator j. The goal is to find such a CD response function and its unknown parameters. Over the years, process control researchers have put considerable effort on modeling the CD response and describing the model in mathematical terms. The modeling method is Chapter 3. Cross Directional Actuator Response Identification 28 generally based on the dispersive wave theory and Darcy's law. The combination of these two theories with appropriate boundary conditions results in having Laplace equations set up for the equation of motion of the slurry on the wire for Fourdrinier paper machines. Solving the Laplace equations gives the solution of the CD response model. For details of the derivations, equations setup and a list of boundaries conditions, see [33]. The CD response shape model is in a form of exponentially damped sinusoidal function which is given by: (3.6) r(x) = OCD e- cos{T\x) ap{x) = [w a 5} . T (3.7) (3.8) In Equation (3.6), g is a normalization constant. The parameter w defines the width of the response. The larger the w, the wider is the response. The parameter 5 is called the divergence parameter which defines the positions of two maximas and the distance between these two maximas. The need for identifying 5 is only applicable for actuators with bi-modal responses which occur with heavy grade papers [22]. Equation (3.7) states the solution of the Laplace equations. The attenuation parameter a defines the size of the negative lobes. If a is small, the response would have prominent negative lobes. On the other hand, the CD response would have little or no negative lobes if a is large. In this thesis, two types of CD responses are considered. In Equation (3.7), the function p(x) can be in one of the following forms: p(x) = |x|, or (3.9) Chapter 3. Cross Directional Actuator Response Identification 29 CD Response Model #1 — r(x) = exp(-a*abs(x/w))cos(pi*x/w) -10 -5 0 Databox —> 5 10 Figure 3.2: CD Actuator Response Model #1: r(x) = e~ ^cos(-Kx). a p{x)=x . (3.10) 2 Figure (3.2) and Figure (3.3) illustrate the response shape with large a and p{x) being \x\ and x respectively. As clearly seen in Figure (3.2), there is a discontinuity present 2 at the maximum for the case when p(x) — |x| while the top area of the response shape shown in Figure (3.3) for the case when p(x) = x is smooth and differentiable. In this 2 thesis, CD identification using two different response models will be illustrated followed by comparisons between these two response functions. The CD response model discussed above can be further simplified so that the computation for the identification is efficient. In this thesis, only the case when the response shape has only one maxima is considered. It implies that the parameter 8 is set to zero. Moreover, the parameter a can befixedto a large value (for example, 7) based on the assumption that no negative lobes present in the response. Therefore, the parameter vector "in Equation Chapter 3. Cross Directional Actuator Response Identification 30 CD Response Model #2 — r(x) = exp(-a'(x/w) )cos(pi*x/w) Figure 3.3: CD Actuator Response Model #2: r(x) = e cos(irx) ax2 (3.8) can be modified as follows: 0CD = [w 7 Op (3.11) Consequently, Equation (3.6) will become: n{x]6 D) = C ge- cos{'Kx) ap{x) (3.12) As a result, the only parameter that needs to be identified is w and g. With the use of wavelets, the values of w and g can be found for each bumped actuator response so that the best fit between the CD response function and the actuator response can be achieved. The concept of using continuous wavelet transform as a medium for correlation analysis is explained in the next section. Chapter 3. Cross Directional Actuator Response Identification 3.2.3 31 Correlation Analysis W i t h Wavelets The algorithm for identification presented in this thesis is based on correlation analysis of the CD response. The correlation is achieved by using the continuous wavelet transform (CWT). In Chapter 2, a brief introduction to wavelet theory has been presented. Similar to Fourier Transform which uses e IUJi as basis functions, wavelet bases are constructed by performing affine operation on a predefined mother wavelet as described in Equation (2.3). The wavelet coefficients are computed using the inner product of the signal and the dilated and translated version of the mother wavelet: VMa, b) =< f(t),i> , (t) >= \a\- ' J °° f(t)^ 1 2 dt + a b (3.13) It is assumed that the mother wavelet function is real and symmetric about zero. For continuous wavelet transform (CWT), the variables t , a and b are continuous so as the integral in Equation (3.13). The CWT in Equation (3.13) can also be viewed as the correlation between the signal f(t) and the mother wavelet basis functions. The mother wavelet functions have to satisfy the admissibility condition described in Equation (2.7) and (2.8). The proof for admissibilty for the two proposed CD response shapes as mother wavelet functions in Equation (3.7) is demonstrated in Appendix A.l. To reconstruct a signal f(t) from the wavelet coefficients computed in Equation (3.13), inverse continuous wavelet transform (ICWT) is used and is defined as a double integral form since the CWT has both time (shift) and frequency (scale) localizations of the signal. The ICWT is given by: fit) = r+oo / J-co flfi fib /.+oo / J—oo W^(a,b)^ (t)^, 0i6 (3.14) ^ where C^ is the inverse of the admissibility function defined in Equation (2.7). 1 From Equation (3.13), it is clearly seen that if the mother wavelet function tp(t) and the signal /(£) are real, the wavelet coefficients are also real [1]. The wavelets used in Chapter 3. Cross Directional Actuator Response Identification 32 the CWT are non-orthogonal. Unlike the case of Discrete Wavelet Transform described in Section (2.3), where scaling functions can be created to form orthogonality with the wavelet basis functions using Multi-resolution analysis and Filter Banks theories, the wavelets used in CWT have no scaling functions. Therefore, perfect reconstruction is not possible. Moreover, another difference between CWT and DWT is the nature of shift and scale parameters. In the continuous case, the shift and scale parameters can be any real values. However, the scale parameter from DWT is usually increased or decreased by a power of 2 and the shift parameter is usually an integer. The nature of CWT makes the correlation analysis possible. Correlation is based on the degree of similarity between two signals: the higher the similarity, the higher is the correlation. The concept of correlation can be applied to CD identification by using CWT. Bump test data can be extracted as an input signal y(t) with the chosen mother wavelet defined in Equation (3.7). The correlation analysis is carried out by correlating the input signal with the dilated and translated versions of the mother wavelet through evaluating the CWT. Large wavelet coefficients show that the degree of similarity between the input signal and the mother wavelet at certain scale and shift is high. This concept will be further illustrated in the next section with industrial bump test data as input signals. 3.3 CD Actuator Response Identification Procedure This section explains the identification procedure based on real industrial data. The industrial bump test data, including the wavelet CD profiles and basis weight dilution flow actuator profiles, are collected in real time using Waveplot. Details on the implementation of the real time diagnostic tool which has been successfully installed at Canfor pulp and paper mill in Prince George will be given in Chapter 4. The bump tests were performed on a 65g/m grade of paper. A set of raw bump test data were first collected from a gauge 2 scanner. The raw data were then further processed using wavelets and the processed data were displayed on Waveplot. Chapter 3. Cross Directional Actuator Response Identification 33 Bump Test Starts • Bump Test Ends take average take average X X X X X I I r 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 take average Scan # Figure 3.4: Extraction of Bump Test Profiles. 3.3.1 E x t r a c t i o n of B u m p Test Profiles f r o m I n d u s t r i a l D a t a The whole bump test procedure lasted for about 21 scans (10.5 minutes). The differential wavelet CD profiles were produced by first averaging the data over 5 scans before the bump test started. The averaged data before the bump test were then subtracted from an average of 5 scans after the bump test started on a sliding window. Figure (3.4) illustrates graphically how the extraction of CD profiles is achieved. The same method applies to the extraction of the actuator profiles. The reason for averaging the profiles is that the effect of the bumps from steady state variations can be isolated [33]. The subtraction of the averaged profile before the bump test from the averaged profile after the bump test is based on the assumption that superposition holds for the papermaking process. Figure (3.5), (3.6), (3.7) and (3.8) show the differential wavelet CD profiles extracted from the bump test. The extraction method results in a total of 16 different profiles. From Figure (3.5) to Figure (3.8), it is obvious that only the last 4 profiles are usable. This is due to the transport delay between the actuators at the headbox and the measuring device at the reel. For demonstration of the identification procedure, the last profile is used as the input signal. Chapter 3. Cross Directional Actuator Response Identification 34 Bump Test CD Profile #1 Bump Test CD Profile #2 Bump Test CD Profile #3 Bump Test CD Profile #4 0 100 200 300 400 500 600 CD Position —> Figure 3.5: Differential Wavelet CD Profiles of a Bump Test (1-4). Bump Test CD Profile #5 I L Bump Test CD Profile #6 Bump Test CD Profile #8 i 0 100 200 1 300 400 500 CD Position—> Figure 3.6: Differential Wavelet CD Profiles of a Bump Test (5-8). Chapter 3. Cross Directional Actuator Response Identification 35 Bump Test CD Profile #9 Bump Test CD Profile #10 Bump Test CD Profile #11 j i i Bump Test CD Profile #12 0 100 200 300 400 500 600 CD Position — > Figure 3.7: Differential Wavelet CD Profiles of a Bump Test (9-12). Bump Test CD Profile #13 Bump Test CD Profile #14 Bump Test CD Profile #15 i 1 Bump Test CD Profile #16 0 100 200 300 400 500 600 Figure 3.8: Differential Wavelet CD Profiles of a Bump Test (13-16). Chapter 3. Cross Directional Actuator Response Identification 3.3.2 36 Peak Search The next step of the identification procedure is to identify the actuator response shape for each bumped actuator individually. In this case, 9 bumped actuators have been used. Therefore, there are a total of 9 different actuator response shapes to identify. It is necessary to extract the 9 different responses from the wavelet CD profile. The method presented here is to use a peak search algorithm to search for local maximas or minimas of the differential wavelet CD profile. Then the responses are extracted with the peaks being the response centers and with specified response widths. The response width of a bumped actuator is determined such that the squared estimation error between the wavelet CD profile and its estimated profile is minimum. This requires a few trial-and-error runs and paper-making experience to validate the result. The results for determining the optimal response width is discussed in Section (3.4). By trial-and-error method, the response width for a bumped actuator is found to be 25 databoxes. The algorithm forfindingthe peaks of the wavelet CD profile and the actuator profile is to find the locations such that the slope at these locations is zero. Mathematically, with differential wavelet CD profile and differential actuator profile described as y(t) and u(k) at scan t respectively, the local maximas or minimas are found such that their first derivatives are zero: y'it) = o, u'(fc) = 0. (3.15) (3.16) Figure (3.9) shows the last set of extracted bump test data of the differential actuator setpoint and its corresponding differential wavelet CD profile with the peaks indicated in crosses. The peak search algorithm assumes that the number of peaks equals the number of bumped actuators. However, as shown in Figure (3.9), the disturbances and noise from Chapter 3. Cross Directional Actuator Response Identification 37 the CD profile could make the algorithm pick a wrong peak. Based on the assumption that the amplitudes of the actuator response are much higher than that of the disturbances, the 9 highest maximas or minimas are chosen as the peaks of the response. Table (3.1) shows the locations of the peaks of the actuator profile in actuator resolution and their corresponding locations of the peaks of the CD profile in databox resolution. Actuator Profile Peak Locations Wavelet CD Profile Peak Locations 7 42 16 107 24 171 32 235 40 299 48 362 55 418 62 474 70 534 Table 3.1: Actuator Profile Peak Locations and Corresponding Wavelet CD Profile Peak Locations. Note that when the differential actuator setpoint is positive, the corresponding peak of the wavelet CD profile is negative. For dilution flow actuators equipped with consistency profiling system, injecting dilution water (positive setpoint) into the header near the headbox tubes results in a reduced basis weight (negative peak) at affected CD positions [31]. The algorithm for searching the peaks for both wavelet CD profiles and the actuator profiles is necessary to check whether the center of the actuator coincides with its response center. If it is not the case, there is a misalignment between the wavelet CD profile and actuator profile. Misalignment is a common problem in paper-making and is caused either by shrinkage at the edges of the paper as the sheet dries in the dry section or by sheet wander. To check for alignment, the peaks of the differential actuator profile are plotted against the peaks of the differential wavelet CD profile. Figure (3.10) shows the alignment Chapter 3. Cross Directional Actuator Response Identification 38 Peaks of Actuator Profile -1 5 I 0 I 100 I 200 I 300 C D Position I 400 : I 500 I 600 > Figure 3.9: Peaks of Differential Actuator Setpoint and its Corresponding Wavelet CD Profile of a Bump Test for Identification. line with crosses indicating the actuator-databox peak pairs. If all the bumped actuators are aligned with the peaks of the CD profile, the shape of the line should be straight. In this case, the alignment in the middle section is satisfactory as shown in the figure. However, there is a misalignment at the edges of the paper machine judging from the fact that there are deviations from the straight line. The alignment plot shown in Figure (3.10) provides mill operators a quick reference on mapping diagnostics. The degree of deviations from the alignment line at the edges also gives useful information visually on how strong the edge effects would be. Chapter 3. Cross Directional Actuator Response Identification 39 Bumped Actuator Positions Vs Wavelet CD Peaks 600 550 500 450 1 i/i 400 350 | 300 Q 250 is 200 150 100 50 10 20 30 40 Actuators —> 50 60 70 Figure 3.10: Actuator-databox Alignment Line. 3.3.3 Correlation Analysis using Continuous Wavelet Transform After searching for the peaks of the differential wavelet CD response, the CD profile is divided into 9 different responses for each of the 9 bumped actuators. As mentioned above, the width for each response is assumed to be 25 databoxes. The peaks found from the peak search algorithm are used as individual response centers. The 9 different response shapes are identified individually using correlation analysis with Continuous Wavelet Transform (CWT). In order to illustrate the identification approach, results from one of the responses are shown. Figure (3.11) shows the 6 extracted normalized bump response. Two sets of th wavelet coefficients are computed with the bump response being the input signal using two different mother wavelets as shown in Figure (3.2) and Figure (3.3). The scale was chosen between 1 and 100. Large wavelet coefficients imply high similarity between the scaled mother wavelet and the bump response. Figure (3.12) shows the wavelet coefficients Chapter 3. Cross Directional Actuator Response Identification 40 represented in a color map for Response Model #1 (p(x) = \x\) as mother wavelet. In the color map, the area where large coefficients occur is between scale 28 and 36. Figure (3.13) illustrates the zoomed color map which shows the wavelet coefficients between scale 20 and 40. A simple search algorithm is carried out to search for the exact scale and shift that gives the largest wavelet coefficient. It turns out that the highest similarity occurs at scale level 32 and zero shift. However, in order to have a finer resolution on the scale level, a new set of wavelet coefficients is computed with smaller step size for the scale level around level 32. That could give the scale level up to 2 decimal places. The high accuracy of the scale level is one of the advantages of using CWT since the scale level and the shift parameters are continuous. As a result, the scale level that gives the best fit between the mother wavelet and the bump response is at 32.45 with shift of zero. The same procedure applies to the computation of the wavelet coefficients using Response Model #2 (p(x) = x ) as mother wavelet. Figure (3.14) shows the color map of wavelet 2 coefficients with scale between 1 and 100. The largest coefficient occurs around scale between 18 and 22 with zero shift. The finer color map depicted in Figure (3.15) and the search algorithm indicate that the scale that gives the best fit is 20.43. The dilated versions of the two mother wavelets at the scales found above are the identified CD response shapes. Figure (3.16) shows the original input bump in dashed line, the identified CD response shape using Response Model #1 in dotted line and identified response using Response Model #2 in solid line. It is clear from the figure that using Response Model #2 provides a better fit compared with the one using Response Model #1 as mother wavelet. The results above combined with results of identification of other bumped responses will be shown and discussed further in the following section. Chapter 3. Cross Directional Actuator Response Identification 11 Normalized Wavelet Bump Response #6 -0.2 0 ' 5 ' 10 ' 15 ' 20 1 25 Databox — > Figure 3.11: Normalized Wavelet Bump Response #6. Color Map for the Wavelet Coefficients of Response #6 with Response Model #1 Figure 3.12: Color Map for Wavelet Coefficients of Response #6 with Scale 1-100 and Response Model #1. Chapter 3. Cross Directional Actuator Response Identification 42 Color Map for the Wavelet Coefficients of Response #6 with Response Model #1 Shifts > Figure 3.13: Color Map for Wavelet Coefficients of Response #6 with Scale 20-40 and Response Model #1. Color Map for the Wavelet Coefficients of Response #6 with Response Model #2 Shifts > Figure 3.14: Color Map for Wavelet Coefficients of Response #6 with Scale 1-100 and Response Model #2. Chapter 3. Cross Directional Actuator Response Identification 13 Color Map for the Wavelet Coefficients of Response #6 with Response Model #2 Shrfts > Figure 3.15: Color Map for Wavelet Coefficients of Response #6 with Scale 20-40 and Response Model #2. Identified CD Response Shapes 1.2 i 1 V / ft It ji it J1/ 08 •a V ////y //.>* 02 0 '0 5 10 15 20 Databox — > Figure 3.16: Identified CD Response Shapes. 25 Chapter 3. Cross Directional Actuator Response Identification 3.4 3.4.1 44 Parameter Estimation Results and Discussions Identification Results The identification procedure with the continuous wavelet transform described above is applied to other wavelet CD responses. Figure (3.17) and Figure (3.18) illustrate the identified CD response shapes for 9 bumped actuators using Response Model #1 and Response Model #2 as mother wavelet respectively. Clearly by visual inspection, the performance for using Response Model #2 as mother wavelet gives a better fit for the bumped actuator responses. The inspection can be further justified quantitatively by evaluating the sum of squared errors for the 9 responses. For identified response denoted as fi(x) and the input wavelet CD response r^x) for each bumped actuator i, the sum of squared error e» is given by: n ei= \ri(x)-fi(x)\ , 2 (3.17) x=—n where x represents CD position with response length of 2n+ 1. Table (3.2) shows the squared errors for 9 bumped actuator responses using two different response models. Except for the third response, the squared errors for using Response Model #2 are lower than that of using Response Model #1 as mother wavelet. In order to obtain thefinaloverall CD response shape, the average of 9 responses for each response model is computed. Figure (3.19) and Figure (3.20) depict the average response for both response models in solid lines and their individual CD responses in dotted lines. The variations of the identified response shapes are high due to the variability of each individual actuator and uncertainty of the identification such as noise and disturbances. The importance of computing the overall average response will be explained later. At this stage the response shape for each bumped actuator has been identified and in turn, parameters w and g which need to be identified from Equation (3.11) have been found. Table (3.3) shows the identified value of w for each response. The values of w for the Chapter 3. Cross Directional Actuator Response Identification 45 Figure 3.17: Identification of CD Responses with Response Model #1 using Continuous Wavelet Transform. Dotted Line:True Response; Solid Line:Identified Response. average responses at the last row of the table are found by performing the correlation analysis discussed above with the average responses being the input signals. Last but not the least, the identified individual responses are assembled to form the estimated CD wavelet profile for the bump test as shown in Figure (3.21). 3.4.2 Response M o d e l Comparisons To compare the identification results using response models illustrated in Figure (3.2) and Figure (3.3), the overlays between the response shape and the bumped response are visually Chapter 3. Cross Directional Actuator Response Identification 46 Figure 3.18: Identification of CD Responses with Response Model #2 using Continuous Wavelet Transform. Dotted Line:True Response; Solid Line:Identified Response. Chapter 3. Cross Directional Actuator Response Identification 47 Figure 3.19: Average Identified Response for Response Model #1. Dotted Linesddentified Responses for 9 Bumped Actuators; Solid LineAverage Identified Response. Chapter 3. Cross Directional Actuator Response Identification 48 Figure 3.20: Average Identified Response for Response Model #2. Dotted Linesddentified Responses for 9 Bumped Actuators; Solid LineAverage Identified Response. Chapter 3. Cross Directional Actuator Response Identification 49 Estimated CD Response for a Bump Test using Response Model #1 with CWT 0 100 200 300 Databox —> 400 500 600 Estimated CD Response for a Bump Test using Response Model #2 with CWT Databox —: Figure 3.21: Identification Results Using Continuous Wavelet Transform. Dotted Line:True Response; Solid Line: Estimated Response. Chapter 3. Cross Directional Actuator Response Identification Bumped Actuator 1 2 3 4 5 6 7 8 9 50 Squared Error for Response Model #1 Squared Error for Response Model #2 (p(x) = |x|) (p(x) = X ) 0.2829 0.1363 1.0654 0.9982 0.1440 0.2249 0.2879 0.1520 0.3272 0.1395 0.6706 0.0312 0.4955 0.1889 0.5182 0.2194 1.5835 0.5970 2 Table 3.2: Squared Estimation Error for Identified CD Response Shapes with CWT. compared. Based on the results from Figure (3.16), (3.17) and (3.18), Response Model #2 gives a betterfitand the results from Table (3.2) confirm the visual judgement. However, the latter model is merely based on extensive experience with paper machine processes [23]. This model has not yet been derived and validated based on physical principles even though the results above suggest the model gives a more accurate estimation. On the other hand, Ghofraniha [33] has derived the CD response model (Response Model #1) using physical principles namely dispersive wave theory and Darcy's law. Therefore, Response Model #1 is supported by physical theories. However, the derivation in the thesis is based on paper machines equipped with slice lip actuators instead of dilution flow actuators. As a result, although Response Model #1 has physical sense, it has not been proved that it can also apply on paper machines equipped with dilution flow actuators. Physical derivations on CD response model for paper machines with dilutionflowactuators are therefore required. 3.4.3 Response W i d t h The response width for the CD response model is chosen such that the sum of squared error is minimum. A few trial-and-error runs have been conducted. Table (3.4) shows the squared error for both response models. They both have a minimum on the squared error when the response width is 25 databoxes. Chapter 3. Cross Directional Actuator Response Identification 51 Bumped w for Response Model #1 w for Response Model #2 (p{x) = \x\) Actuator (p[x) = X) 1 51.06 25.29 2. 42.30 25.38 3 51.11 26.01 4 34.77 20.99 5 43.35 25.38 6 32.45 20.43 ' 7 25.19 15.90 8 51.11 25.91 9 42.55 25.43 Average 41.00 25.38 2 Table 3.3: Parameter w for Identified CD Response Shapes with CWT. 3.4.4 Significance of Identified Response Shape The identified actuator response shape is useful to determine the control bandwidth of the paper machine. A general method in determining the CD controllable components is to find out which spatial frequency separates the CD profile into controllable and non-controllable components. In paper-making industry, it is assumed that a CD controller could remove disturbances with a wavelength longer than two actuator spacings, 2X . In other words, a the control bandwidth is estimated to be up to the Nyquist frequency f which is defined c as: For this particular paper machine in study, the Nyquist frequency for the basis weight paper property is: / c = 2(fe) « 0.006mm- 1 (3.19) However, this rough assumption does not consider practical issues such as model uncertainty, robustness and actuator limit. According to the identification results discussed Chapter 3. Cross Directional Actuator Response Identification 52 Response Squared Error for Response Model #1 Squared Error for Response Model #2 Length (p(x) = (p(x) = X ) 15 16.7487 12.8830 17 15.6788 12.5004 19 14.3338 12.9421 21 14.3909 12.7323 23 14.3787 12.6406 25 14.1863 12.2424 27 14.5292 12.2705 29 14.3775 12.5143 31 14.6098 12.4994 33 14.3032 12.5363 35 14.3918 12.7440 2 Table 3.4: Sum of Squared Error of Estimation Error with Different Response Lengths. previously, the model uncertainty mostly comes from variations of the identified response shapes and uncertainty of identification. If the control bandwidth was set up to the Nyquist frequency, the large model uncertainty could drive the system to instability resulting in having a picketing actuator profile. Moreover, the considerable model uncertainty could cause the process noise or disturbances with frequencies around the Nyquist frequency not being filtered out properly so that they would enter into the controller. In fact, it is not always a good idea to control noise and disturbances as the controller could amplify the noise and disturbances which could lead to catastrophic outcomes. Therefore, it is wiser to choose the control bandwidth conservatively to accommodate large model uncertainty. The way of estimating the control bandwidth is to consider the Bode plots of the identified response shapes as well as the power spectral densities of the CD profile in spatial frequency domain. Figure (3.22) and Figure (3.23) illustrate the Bode plots of identified actuator responses for Response Model #1 and #2 in logarithmic scale. Table (3.5) highlights important parameters for the Bode plots of the average identified response shapes for both models that are useful for control performance analysis. The figures shown confirmed that it is inap- Chapter 3. Cross Directional Actuator Response Identification 53 Figure 3.22: Frequency Response for the Identified Actuator Responses Using Model #1. Blue Line: Average Response: Red Lines: All Bumped Actuator Responses. propriate to choose the Nyquist frequency as the control bandwidth which is the frequency for separating the CD profile into controllable and non-controllable components. According to Table (3.5), at the Nyquist frequency, the magnitudes for both models are around 30dB below the 3dB magnitudes. If the control bandwidth were set at that frequency, an amplifier of 30dB is needed to bring up the magnitudes which is quite large. Given such amplifier with high gain, the noise and disturbances accompanied with considerable model uncertainty could be largely amplified resulting in instability of the control system. The most conservative approach for choosing the right spatial control bandwidth is to select the 3dB frequency. The 3dB frequencies for the average responses using Response Model #1 and #2 are 0.0018 and 0.002 mm -1 respectively. However, this choice suffers large model uncertainty ranging from 17 to 20dB due to extensive variations of individual actuator responses. Control bandwidth or control cutoff frequency should be chosen where the model uncertainty is the minimum with the process gain being high enough. From Figure Chapter 3. Cross Directional Actuator Response Identification 54 Figure 3.23: Frequency and Phase Response for the Identified Actuator Responses Using Model #2. Blue Line: Average Response; Red Lines: All Bumped Actuator Responses. Chapter 3. Cross Directional Actuator Response Identification Parameters Response Model #1 Response Model #2 (p(x) = \x\) hdB (Spatial 3dB Frequency in mm~ ) M B (3dB Magnitude in dB) M (Magnitude at Nyquist Frequency in dB) f (Crossover Frequency in mm' ) 6 (Phase Margin in Degrees) l M ny 1 c m 55 0.0018 40.34 10.5 0.010 72 {p{x) — X) 2 0.002 43.75 17 0.007 101 Table 3.5: Important Data for Frequency and Phase Plots for Averaged Identified Actuator Responses. (3.22) and (3.23), a reasonable choice for control bandwidth is within a range between 0.0045 and 0.0048 mm -1 or 2.54 and 2.71 actuator spacings since the model uncertainty and the amplification gain required are around 10 to 12dB which is a more appropriate range. Therefore, the choice of the control bandwidth or control cutoff frequency can be a complicated task because all the constraints imposed by the process model have to be considered. The control bandwidth should also be chosen based on the power spectral densities of the CD profile. The computation of the power spectral densities is explained in Appendix (A.3). Figure (3.24) shows the power spectral densities for raw, wavelet and residue profiles in logarithmic scale. The method for computing wavelet and residue profiles will be described in the next section. The advantage for using wavelets for analyzing paper machine data is that white noise can be eliminated by decomposing the raw data into several wavelet levels followed by denoising. With the help of the power spectral densities plot, the chosen range of cutoff frequencies can be verified. It is assumed that white noise has spatial frequencies larger than the Nyquist frequency. From Figure (3.24), raw CD profile is decomposed up to wavelet level 3. Next section has details on how to choose the right level of wavelet decomposition. The analysis of Bode plots above concludes that the cutoff frequency for separating CD profile into controllable and non-controllable components is between 0.0045 and 0.0048 Chapter 3. Cross Directional Actuator Response Identification Power Spectral Density for Raw. Wavelet and Residue Profiles {Basis Weight) 140 I 10 56 1 . . 10" 5 . . , 10~ Spatial Frequency (mm~1) 3 , 10~ 2 . . . . . 10"' Figure 3.24: Power Spectral Densities for Basis Weight. mmr or 2.54 and 2.71 actuator spacings. When referred to the power spectral density 1 plot, this range has a good agreement with the discussions above since around that range the residue has little power densities and the wavelet profile covers over an estimation of 90% of the area under the raw profile curve. This confirms that the range of the cutoff frequencies has been reasonably chosen. 3.5 Summary In this chapter, the CD actuator response identification using continuous wavelet transform was considered. The formulation for the overall identification problem was first developed. The identification procedure was then explained and demonstrated by using industrial data obtained from Canfor pulp and paper mill in Prince George. Two different CD response shapes were used as mother wavelets. The estimation results using each response model were illustrated and compared. The second response model (p{x) = x ) was found to 2 Chapter 3. Cross Directional Actuator Response Identification 57 be a better response model but it lacks physical validation. Therefore, a response model derived from physical principles for paper machine equipped with dilution flow actuators is required. Finally, the choice of response width and the significance of the identified CD response shape were discussed. The control cutoff frequency that separates the CD profile into controllable and non-controllable components falls into the range between 0.0045 and 0.0048 mm -1 or 2.54 and 2.71 actuator spacings. Chapter 4 Real Time Paper Machine Data Wavelet Analysis 4.1 Introduction The development of an online automated diagnostics tool for analyzing paper machine data is useful to mill operators when monitoring the performance of a cross directional (CD) control system. The use of Object linking and embedding for Process Control (OPC) interface as a medium for client-server based data acquisition provides an efficient, robust and accurate data transfer between the CD control system and a monitoring workstation. Furthermore, there is a need to develop a real time diagnostics toolbox to help paper mill companies reduce costs and save time on troubleshooting. This chapter describes a real time wavelet toolbox known as Waveplot which has been under development at UBC for several years [62]. The outline of this chapter is to first give a brief overview on OPC. Then a detailed description on the main features of Waveplot such as color maps display, performance indices as well as controllable and non-controllable standard deviations will be given. Specifically, the performance indices and standard deviations will be re-calculated according to the knowledge of the identified CD actuator response shape discussed in Chapter 3 to reflect recent views on CD controllability. Moreover, special features such as turn-up marks and real-time display on workstations will be described. The displays are demonstrated by collecting paper profiles through the OPC server from a multi-grade, operating paper machine. 58 Chapter 4. Real Time Paper Machine Data Wavelet Analysis 4.2 59 Overview of Object Linking and Embedding for Process Control (OPC) The use of application software and microprocessors in manufacturing plants could increase productivity given that these digital devices can communicate with each other, but it is not often the case. Before OPC was adopted by the automation industry, application software suppliers had to develop their own interfaces for each of their applications which is very time consuming and costly. OPC opens the door for open connectivity with open standards [17, 16]. With a series of standards specifications, a standard set of objects, interfaces and methods used in process control and manufacturing automation applications can be defined so that interoperability between these applications can be guaranteed. The best analogy for needing a standard would be printer drivers in DOS and then in Windows [16]. In the past software application suppliers had to write their own printer driver for every printer. For example, software developers who wrote AutoCAD had to write the printer drivers besides writing the AutoCAD application. Similarly, developers who wrote WordPerfect had to write software to interface the application with the printers such as HP LaserJet or Epson Bubble Jet printers in addition to writing their core application and so on. In industrial automation industry, for instance, GE Fanuc wrote their Human Machine Interface (HMI) and a list of drivers for the industrial devices such as every brand of Programmable Logic Controller (PLC) [16]. The time spent on writing printer drivers or drivers for industrial devices was endlessly long so that application software suppliers added the cost of labour on writing the drivers to end-users. Windows solved the printer driver problem by including the driver support in their operating systems. With the standardization of printer driver, printer manufacturers not the application developers write the printer drivers for their products by following the printer driver standard. Similarly, under the OPC standard, the industrial device manufacturers can write their own OPC Data Access Servers and Clients without worrying about the interoperability with other industrial devices. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 60 The OPC client/server architecture is shown in Figure (4.1) [60]. Different vendors such as Foxboro can provide OPC servers. Vendors supply code that determine which server has access to the devices and data [18]. They also determine the data names and other details about how the server physically accesses the data. The OPC clients can have access to different OPC servers provided by the vendors. The clients could be applications such as Matlab or other data acquisition devices. The application developers write code to instruct the clients to write or read data to or from the servers. 4.3 Waveplot for Real Time Paper Machine Data Analysis Tool 4.3.1 R e a l T i m e D a t a Access Methodology The use of OPC provides efficient data transfer so that it is suitable for real time data access from the paper machine and the Distributed Control Systems (DCS). The simple data access architecture is shown in Figure (4.2). The raw basis weight and moisture data as well as other miscellaneous data are collected using the ABB Smart Platform scanner and are stored in the Foxboro IA DCS. An OPC server is installed to allow data delivery to the clients from the server or writing data from the clients to the DCS. Waveplot is developed and implemented using the Matlab Graphical User Interface Development Environment (GUIDE) as it can provide satisfactory graphics illustrations and straightforward as well as intuitive interfaces for users. The OPC Client for Matlab developed by IPCOS is used as the interface between the OPC server and Matlab. This Matlab OPC client software enables users to transfer data to or from the OPC server through direct transfer with the read and write commands or transfer via the cache of the OPC client. Figure (4.3) illustrates the two data transfer principles [57]. As shown in the figure, users can either use read/write commands to retrieve/write the Matlab variables or use readcache/writecache commands to retrieve/write the Matlab variables indirectly Chapter 4. Real Time Paper Machine Data Wavelet Analysis Figure 4.1: OPC Client/Server Relationship. 61 Chapter 4. Real Time Paper Machine Data Wavelet Analysis 62 Foxboro IA D C S O P C Server O P C Interlace Client #1 Client #2 Figure 4.2: Real Time Data Access Diagram. Client #3 Chapter 4. Real Time Paper Machine Data Wavelet Analysis 63 OPC Server Read Write Matlab ReadCache WriteCache OPC Clidnt Cache Read/Write Matlab Variable Figure 4.3: Diagram for Data Transfer Between OPC Server and Matlab Client. to or from the OPC server. One of the advantages of data transfer via cache compared with the direct transfer is a faster transfer. For instance, if one wants to transfer over 1000 tags or variables to the OPC server, it only takes less than 1 second to transfer the data via cache while it could take more than 30 seconds by direct transfer. To read data through cache from the OPC server, only one readcache command is used to retrieve all the required data from the server to the cache. Then with the normal read commands the data can be retrieved from the cache to the Matlab environment. For writing data to the OPC server, Matlab variables are first loaded into the cache by using regular write commands. All the variables stored in the cache are then written to the OPC server by using only one writecache command. In the case where the data from the paper machine are acquired from the DCS, only readcache and read commands are used. The required data that are needed to obtain from the OPC server for wavelet data analysis in Waveplot are databox values for basis weight and moisture of each scan, setpoint values for basis weight and moisture actuators as well as the MD stockflow,reel turn-upflagand sheet break flag. Table (4.1) lists the required Chapter 4. Real Time Paper Machine Data Wavelet Analysis 64 data and the number of variables which are retrieved from the OPC server. There are a total of 1321 tags that are needed to retrieve from the OPC server and therefore, using the cache for data transfer with one readcache command is chosen. Note that each databox or each actuator setpoint for both paper properties is assigned with a tag name. Number of Paper Machine Data Tags Basis Weight 600 Moisture 600 Basis Weight Actuators 78 Moisture Actuators 40 MD Stock Flow 1 Turn-up Flag 1 Sheet Break Flag 1 Total 1321 Table 4.1: Required Paper Machine Data and the Number of Tags. Figure (4.4) illustrates a flowchart for the data retrieval process. The first three steps consist of first connecting to the OPC server (local or remote). The tags are then created based on the syntax according to the OPC standards specifications. The actions (read or write) of the tags are then defined. The loop after the data definition stage is to read the paper machine data from the OPC server and then perform data analysis using wavelets. Since a new set of data arrives the OPC server every 30 seconds, a command for synchronizing the Matlab client with the server is therefore required. A counter variable is also created to record the number of scans of paper machine data collected. In Waveplot, the graphics including the color maps and performance indices are refreshed every 300 scans. Chapter 4. Real Time Paper Machine Data, Wavelet Analysis 65 Connect to the O P C server Build tag names Define read/write attributes for the tags # of scans: 300? Yes Refresh the Waveplot graphics Reset the scan count No Synchronize the Matlab environment with the O P C server 1 Read data from O P C server to O P C client Perform wavelet data analysis t Increment scan count Figure 4.4: Flowchart for the Data Retrieval Process. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 4.3.2 66 Overview of M a i n Features of Waveplot This section gives a detailed description on Waveplot as real time diagnostics tool. As mentioned above, the application is created using the Matlab Graphical User Interface Development Environment (GUIDE). Figure (4.5) illustrates the GUI window of the toolbox. Waveplot allows users to choose the options within the 7 categories for analyzing paper properties. It gives users flexibility in choosing wavelet parameters to observe different characteristics of the paper profiles. These options are: 1. Display • Basis Weight • Moisture 2. Wavelet Type • Daubechies • Symmlet • Coiflet 3. Thresholding Method • Hard • Soft 4. De-noising Method • Std • MAD • MultiMAD 5. Wavelet Level • From 1 to 4 Chapter 4. Real Time Paper Machine Data Wavelet Analysis 67 6. Wavelet Length • Daubechies - 2,4,6,8,10,12,14,16 • Symmlet - 6,8 • Coiflet - 2,3 7. Plot • Raw CD Profile • Wavelet CD Profile • Residue Profile • Standard Deviations All the features of Waveplot will be described in the next section. 4.3.2.1 Color Maps Display The main goal for CD control is to achieve uniformity of the paper produced. Color maps provide mill operators visual observations of paper properties. In Waveplot, there are three types of color maps namely raw, wavelet and residue profiles. Raw CD profiles pit) are computed from the pre-processed data collected from the scanners which are given by p(t) = y(t)-y(t), ten (4.1) where y(t) is the pre-processed paper properties collected from scanning devices. y{t) is the mean value of y{t) at scan t. The CD wavelet profile is computed according to multiresolution analysis and filter banks theories discussed in Chapter 2. The de-noising and thresholding methods described in Chapter 2 are used for noise reduction. The residue of a scan of data is computed as the difference between the wavelet and the raw profile. Mathematically, Figure 4.5: Waveplot GUI Window. Chapter 4. Real Time Paper Machine Data Wavelet Analysis m Raw Image for Basis Weight on Nov 12 23:38pm 500 600 Figure 4.6: Raw Image for Basis Weight Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) Pr(t)=p(t)-P(t), (4.2) where pit) denotes the wavelet CD profile and p {t) is the residue profile. The residue r profiles represent variations in high frequencies and are used to observe non-controllable patterns. Fi gure (4.6), (4.8) and (4.10) illustrate the raw, wavelet and residue basis weight profile respectively. Figure(4.7), (4.9), (4.11) show the profiles for the moisture property. The images were all collected on November 12, 2003 at 23:38pm on a 74.5 g/m grade of paper 2 using Waveplot. The moisture wavelet profile shows more non-uniformities of the moisture content with apparent streaks at the paper edges. Both residue profiles for basis weight and moisture show a few streak lines around databox location 150. The streaks imply there are larger high frequency variations around that area. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 70 Raw Image for Moisture on Nov 12 23:38pm 100 200 300 400 500 Figure 4.7: Raw Image for Moisture Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) Wavelet Image for Basis Weight on Nov 12 23:38pm 100 200 300 400 500 Figure 4.8: Wavelet Image for Basis Weight Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) Chapter 4. Real Time Paper Machine Data Wavelet Analysis 71 Wavelet Image for Moisture on Nov 12 23:38pm 250 100 200 300 400 500 600 Figure 4.9: Wavelet Image for Moisture Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) Residue Image for Basis Weight on Nov 12 23:38pm 100 200 300 400 500 600 Figure 4.10: Residue Image for Basis Weight Collected on December 2, 2003 at 23:38pm. (Scan Vs Databox) Chapter 4. Real Time Paper Machine Data Wavelet Analysis 72 Residue Image for Moisture on Nov 12 23:38pm 0 100 Figure 4.11: Residue Image for Moisture Collected on November 12, 2003 at 23:38pm. (Scan Vs Databox) 4.3.2.2 M D and CD Performance Indices One of the most important applications of wavelet multiresolution analysis is performance monitoring and control performance assessment. The parameter that indicates the quality of control performance is the performance index. With the help of multiresolution analysis, CD profile can be separated into controllable and non-controllable components. The CD performance index can be calculated and is given by [63] °CD where °contr is the variance of the CD wavelet profile and a 2 <J CD ontr is the variance of the controllable components. In order to determine the wavelet level that separates the controllable and non-controllable components for calculating the CD performance index, the range of wavelengths for approximation and detail of each level in terms of actuator spacing is calculated. The spatial Chapter 4. Real Time Paper Machine Data Wavelet Analysis 73 wavelength A is defined as ? n i+1 < 4 ' 4 ) where i is the wavelet level, n is the number of actuators and n is the number of databoxes. d a For this specific paper machine in study, there are 600 databoxes and 78 actuators for basis weight as well as 600 databoxes and 40 actuators for moisture property. Therefore, the wavelengths (in actuator spacing) for basis weight and moisture properties are BasisWeight : A = 0.13 • 2 Moisture : A = 0.067 • 2 i+1 i+1 (4.5) (4.6) For discussion convenience, only the CD variation separation for basis weight is considered. The same concept applies for analyzing CD moisture content. Figure (4.12) illustrates the wavelet decomposition tree and the range of wavelengths for each node. In Chapter 3, the.control cutoff frequency is calculated and it falls within a range between 2.54 and 2.71 actuator spacings. From Figure (4.12), the level that separates controllable and noncontrollable variations would be 3 since Node(3,0) in the decomposition tree has the cutoff wavelengths. Therefore, components above this level are regarded as controllable. However, the separation using merely wavelet decomposition is still not accurate since some components which are not controllable (with wavelengths between 2.08 and 2.54 actuator spacings) are improperly regarded as controllable. Wavelet packets could be used to further divide the wavelengths within one wavelet level for finer frequency resolution [34]. Despite the inaccuracies, the basic wavelet decomposition used in this case provides a good approximation of the "true" CD performance indexes obtained if wavelet packets were used instead. The algorithm for calculating MD performance index is shown in Appendix A.2. Figure (4.13) shows the Waveplot GUI while it was collecting data for wavelet analysis. The Chapter 4. Real Time Paper Machine Data Wavelet Analysis 74 (0,0) 0.26 - Inf (1,0) 0.52-Inf (2,0) 1.04- Inf (3,0) 2.08 - Inf (4,0) 4.16- Inf (1,1)0.26-0.52 (2,1)0.52-1.04 (3,1) 1.04-2.08 (4,1)2.08-4.16 Figure 4.12: Wavelet Decompostion Tree and the Wavelength for each Node. color map of basis weight is shown in the middle plot. The wavelet image gets incremented every 30 seconds when a new scan is collected from the OPC server. The plot on the right hand side shows the CD performance index. It also gets updated the same time as when the color maps get updated. The MD performance index is shown on the left hand side of the Waveplot window. 4.3.2.3 Controllable and Non-controllable Standard Deviations Besides the CD performance index, the standard deviations of controllable and noncontrollable components can also provide good control performance assessment. These components are determined by the CD profile separation using wavelet decomposition described above. By plotting the two components on a same graph, two sets of standard deviations can be visually compared. Normally, the standard deviation of the controllable components should be much lower than that of the non-controllable components for good control performance. Figure (4.14) illustrates the standard deviations plot in the middle of Chapter 4. Real Time Paper Machine Data Wavelet Analysis Figure 4.13: Wavelet Color Map and MD and CD Performance Indices. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 76 the GUI once users selected the standard deviations option under the "Plot" popup menu located at the top of the window. 4.3.3 Special Features 4.3.3.1 Time Stamps Time stamps are used to give time information to mill operators when certain events occurred. They are beneficial to help mill engineers trace back any control problems that occurred in the past. The events for which time stamps are recorded include: • current scan received (located at the second edit box below the push buttons) • every 50 scans collected (located between the MD performance index axes and the r color map axes) • last 300 scans collected (located at the first edit box below the push buttons) • last scan number that a turn-up event happened (located at the last edit box) Figure (4.15) illustrates the Waveplot GUI with all the time stamps shown in all the edit boxes. Special attention is paid to the 3 edit boxes below the push buttons. In this example, the last 300 scans were collected and saved in the historical database around 10:18am on November 12, 2003. The current scan number was at 295 which was received at 13:09pm. Furthermore, the time stamps recorded for every 50 scans were at 10:47, 11:09, 11:44, 12:13 and 12:43pm. 4.3.3.2 Menu Bar In addition to the color maps, plots of performance indices and standard deviations, Waveplot also includes other miscellaneous yet important plots to support the main features. After pressing Start to run the application, there will be a menu bar displayed at the top of the GUI. The menu bar contains several options for mill operators to choose. The available plots for basis weight and moisture properties in the menu bar are: Chapter 4. Real Time Paper Machine Data Wavelet Analysis Figure 4.14: Standard Deviations. 77 wavejilntjnain Cross-Direction Machine-Direction Display || Basis Weight Wavelet Type ^\ Wavelet Length s | Syrnrniet . . ' Actuator Profiles " • .. Plot ^V/avete* CD P r o f i l e " Misc. Plots About Method MI. MAD j . IMuMMAD Thresholding Method H 1 Wavelet Level Denoisirg Quit Ii Crosshair 1 Axes Setup] _Jb!*gBii2ir^^ Cuient Scan Numbei:295 TimeFleceived:1309:56 Datall/KWCKra Last Turnup At Scan ft 295 Paper M achine S tatus: RUNNING Wavelet B a s i s Weight Image Chapter .4. Real Time Paper Machine Data Wavelet Analysis 79 • Cross Direction - CD profiles (raw, wavelet and residue) • Machine Direction - MD profiles (scan average) • Actuator Profiles - actuator setpoints in bar form • Misc. Plots - power spectral densities and 2-sigma plots The power spectral densities of the CD profiles are computed using Welch's averaged periodogram method [43]. The mathematics for this method is explained in Appendix A.3. The 2-sigma is twice the standard deviation of the CD profiles (raw, wavelet and residue). Figure (4.16) shows a plot of basis weight actuator setpoints after a user clicked on the options in the menu bar located at the top. 4.3.3.3 Crosshair . The crosshair function allows mill operators to obtain values at specific coordinates on the color maps. This gives the operators quick reference quantitatively besides color visualizations. After clicking the crosshair push button, a cross will appear on the GUI to help users locate the desired coordinate. A message box will pop up to inform users the corresponding value of the location chosen. Figure (4.17) illustrates the message box. 4.3.3.4 Turn-up Marks and Machine Status Reel turn-up refers to the event when a new reel is started at the end of the machine after collecting 10 tons of paper (approximately 120 scans of data). The turn-up event is detected by a turn-up flag. The turn-up is started by first toggling the flag from 0 to 1. By the time when the turn-up finishes, the flag then switches from 1 back to 0. This event is recorded as tick marks located at the right hand side of the color map axes in the GUI. As shown in Figure (4.15), the turn-up events happened at scan number 25,137 and 205 which are displayed on the axis of the color map. Moreover, when the paper machine stopped running, a message will be displayed to notify mill operators and all the computations will Figure 4.16: Basis Weight Actuator Setpoints. Figure 4.17: Crosshair. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 82 Figure 4.18: Mill Workstation. be halted. Plots such as color maps, performance indices and 2-sigma will not be updated until the machine starts running again. 4.3.3.5 Real Time Display on Workstations Mill operators may need the most updated information on the analyzed machine data for online control performance assessment. Therefore, it is necessary to display the updated color maps and plots on the mill workstations like the one shown in Figure (4.18). Important information such as wavelet color maps and performance indices are displayed on a webpage for the operators. The webpage is refreshed whenever the data are updated. 4.3.3.6 Historical Storage In order for the mill engineers to trace back any control problems that occurred in the past, a database has been built for offline diagnostics. For every 300 scans of analyzed data received, the graphics and data with time stamps will be saved in the database. This database is beneficial to customer service since when customers filed complaints to the Chapter 4. Real Time Paper Machine Data Wavelet Analysis 83 pulp and paper mill about the poor paper quality produced by the paper machine, mill engineers need to know if the reason for producing such a poor quality of paper is actually caused by poor control performance or other post paper machine processes. By reviewing the related color maps and performance indices on specific time frame, mill engineers can determine if the complaints from the customers are justifiable. Furthermore, the database can also help engineers improve control performance by comparing the data from the past with the recent ones. 4.3.4 A R e a l T i m e Diagnostics Example A real life diagnostics example from the Canfor pulp and paper mill in Prince George is presented to show the effectiveness of Waveplot as online diagnostics tool [32]. When the paper machine started running again after the maintenance shutdown on December 1, 2003, the machine had very poor CD performance. A series of bump tests showed that there were weak or no responses in the middle of the headbox. Moreover, color maps collected using Waveplot showed that the control performance at the edges of the CD basis weight was very poor. Figure (4.19) and (4.20) illustrate the wavelet color maps of basis weight collected on November 10 and December 2 ,2003 respectively and Figure (4.21) and (4.22) show the wavelet color maps for moisture on the same dates. The wavelet color maps collected on December 2, 2003 obviously show that the control performance was poorer especially at the edges indicated by the red streaks shown in the basis weight color map. Mill engineers got notified immediately after the results displayed on Waveplot so that investigation into the problem can be launched. However, they did not know if the cause of the problem was due to CD or MD control, or both. Therefore, mill engineers suggested comparing the CD and MD performance indices collected recently with the ones collected a month ago which have been stored in the database. Figure (4.23) and (4.24) show the MD performance indices while Figure (4.25) and (4.26) illustrate the CD performance indices collected on the same dates. The MD performance indices shown in the figures are about Chapter 4. Real Time Paper Machine Data Wavelet Analysis 81 Wavelet Image for Basis Weight on Nov 10 0 8 : 1 0 a m Figure 4.19: Wavelet Color Map for Basis Weight Collected on November 10, 2003 at 08:10am. (Scan Vs Databox) Wavelet Image for Basis Weight o n Dec 2 17:23pm T W so H i I I Figure 4.20: Wavelet Color Map for Basis Weight Collected on December 2, 2003 at 17:23pm. (Scan Vs Databox) Chapter 4. Real Time Paper Machine Data Wavelet Analysis 85 Wavelet Image for Moisture on N o v 10 0 S : 1 0 a m 400 500 600 Figure 4.21: Wavelet Color Map for Moisture Collected on November 10, 2003 at 08:10am. (Scan Vs Databox) Wavelet Image for Moisture on D e c 2 17:23pm 13 1 *> ( * | ( '• 0 1. •- 100 • 1 200 j 1 -J 300 400 :I : 500 600 Figure 4.22: Wavelet Color Map for Moisture Collected on December 2, 2003 at 17:23pm. (Scan Vs Databox) Chapter 4. Real Time Paper Machine Data Wavelet Analysis 86 Figure 4.23: MD Performance Indices for Basis Weight Collected on November 10 and December 2, 2003. the same but the CD performance indices collected on December 2, 2003 is at least three times higher than that of collected the month before. From the performance indices, it is deduced that the poor performance occurred was caused conclusively by a CD problem. From the mill engineers' experience, they suspected that the problem was originated from the hardware of the headbox. After a thorough investigation, they found that the flow sheets were much too short with the confirmation from Mitsubishi Heavy Industries that a wrong type of flow sheets had been shipped. As mentioned in the introduction (refer to Figure (1.2)), the stock is discharged from the tube bank and is injected into the thin converging channels in the MH Headbox-C IV Chamber. The installation of flow sheets is used to prevent excessive stock mixing in the cross direction coming from the upstream so that the stock can be evenly distributed across the machine. The length of theflowsheets affects the degree of turbulence of the stockflow.With theflowsheets being shorter than designed, the degree of turbulence was not in a right level which caused excessive stock mixing at the edges resulting in the saturation of CD actuators in that area. Consequently, Chapter 4. Real Time Paper Machine Data Wavelet Analysis 87 M D Pen". Index f o r M o i s t u r e o n N o v 1 0 0 6 : 1 0 p m M D Pen*. I n d e x f o r M o i s t u r e o n D e c 2 1 7 : 2 3 p m Figure 4.24: MD Performance Indices for Moisture Collected on November 10 and December 2, 2003. Figure 4.25: CD Performance Indices for Basis Weight Collected on November 10 and December 2, 2003. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 88 C D Perf. Index f o r M o i s t u r e o n N o v 1 0 0 S : 1 0 a m Figure 4.26: CD Performance Indices for Moisture Collected on November 10 and December 2, 2003. the CD edge control was nullified. This in turn affected the middle actuators and that was why weak or no actuator responses in the middle area of the machine were observed when bump tests were conducted. From this diagnostics example, Waveplot helped Canfor's mill engineers in a number of ways. Firstly, the real time displays of color maps, performance indices and other plots helped them discover the problem right away so that investigation could be launched in an efficient fashion. The historical database also helped the engineers retrieve information from the past as benchmarks so that comparisons with the recent color maps and performance indices could be possible. Last but not the least, Waveplot helped the engineers narrow down the cause of the problem into a CD control problem. This considerably saved lots of time in troubleshooting. Chapter 4. Real Time Paper Machine Data Wavelet Analysis 4.4 89 Summary In this chapter, a detailed description of a real time diagnostics tool known as Waveplot using OPC was given. A brief introduction to OPC and the basic data communication methodology were also given. The main features such as color map displays and performance indices as well as other supporting features were then described and shown by screen captures of the GUI window. A diagnostics example was demonstrated to show the effectiveness of implementing Waveplot in the Canfor pulp and paper mill. Chapter 5 Conclusions and Future Work 5.1 Conclusions Results of the CD actuator response identification and detailed descriptions on Waveplot as online automated diagnostics tool were presented. For CD actuator response identification, a parametric model for CD actuator response developed in [23] was used. The approach for the identification problem using continuous wavelet transform was explained. The procedure to identify the CD response shape by using two different CD response models was then described. The results of using both models were illustrated and compared. For real time diagnostics, Waveplot has been developed to help troubleshoot paper machine MD and CD control systems and to implement at the computer workstations from the Canfor pulp and paper mill in Prince George. Data acquisition using OPC was first explained. Main and other supporting features in the online automated diagnostics tool were then described. A real time diagnostics example using Waveplot was also presented to demonstrate the effectiveness of the application. The highlights of the major results of this thesis are in the following: • From the alignment plot illustrated in Figure (3.10), although the edge effect exists, it was found that the edge effect for this particular paper machine is not strong. The deviations of the plotted mapping line at the edges from the nominal alignment line are relatively low. • CD actuator response identification using continuous wavelet transform with two different response models (Response Model #1: r(x) = e~ ^cos(-Kx) and Response a 90 Chapter 5. Conclusions and Future Work Model #2: r(x) = e~ ax COS(TTX)) 91 has been demonstrated. It was found that Response Model #2 gives a better performance by using visual judgement and comparing estimation errors. However, the response model lacks physical validation for this particular paper machine equipped with dilution flow CD actuators. • The general assumption that the controllable components with frequencies lower than the Nyquist frequency or with wavelengths longer than 2 actuator spacings does not consider the model uncertainty, robustness as well as the actuator limit. To determine of the control bandwidth, Bode plots of the bumped actuator responses and the power spectral densities for raw, wavelet and residue profiles have to be considered. The reasonable control bandwidth for this particular paper machine in study was found to fall in a range between 0.0045 and 0.0048 mm or 2.54 and 2.71 1 actuator spacings. This revised control cutoff frequency has a considerable impact on recent views for defining controllable and non-controllable components separated from CD profiles and how the performance indices are calculated. • The use of Waveplot as an online automated diagnostics tool provides color maps for raw, wavelet and residue profiles for basis weight and moisture paper properties. The OPC interface as a medium for client-server based data acquisition provides efficient data transfer to make real time wavelet data analysis possible. Other important features in the application such as performance indices, power spectral densities and standard deviations for controllable and non-controllable components give mill operators a complete assessment of control performance. • The real time diagnostics example provided in this thesis demonstrated that Waveplot can help paper mill companies save time and costs on troubleshooting. Chapter 5. Conclusions and Future Work 5.2 92 Future Work There is some future work that can be carried on after this research. These tasks are in the following: • Dilution Headbox Modelling: Although identification results suggested that the Response Model #2 gives a better fit for the actuator responses, the model needs to be validated physically. With the paper machine equipped with consistency profiling headboxes, physically based modelling of the actuator response could provide a more accurate actuator response identification. • Verification of the Superposition Assumption: In Chapter 3, it is assumed that the ' process dynamics is linear so that the the system response can be viewed as the superposition of all the actuator responses. However, this assumption has yet to be verified by using industrial data. • Real Time Implementation of CD Variation Separation Using Wavelet Packets: The use of wavelet packets for separating CD profiles into controllable and non-controllable components is beneficial since it givesfinerfrequency resolution. Although Jiao [34] already suggested using wavelet packet analysis for CD variation separation, the algorithm has yet to be implemented in real time with additional computation costs. Appendix A A.l Proof of Admissibility of Mother Wavelets The admissibility condition for a mother wavelet function tp(x) was defined in Chapter 2 as: f \lb(LU)\ + c o 2 ^f^(iw<oo (A.l) M V-oo In this thesis, two CD response models have been used in Chapter 3 as mother wavelets. Admissibility has to be proved for each of the response model. A.1.1 ip(x) = e~ ^cos(-Kx) a The Fourier transform for ip(x) is given by: +oo / e- ^cos{nx)e- dx a (A.2) juJX • CO />0 /•+oo = - / e- cos(irx)e- dx ax a + , jux u (a + joj) + 7T 2 2 + ( . + / _, (a — jto) + 7T 2 Next step is to substitute Equation (A.4) into (A.l): 93 e cos(irx)e- dx ax juJx (A.3) (a.4) Appendix A. 94 +oo / ,.,2 8 a ( a +, w + 7r f+OO / 2 (A.5) -du M •oo 2 2 2\2 V w ( ( a + w + vr ) + 27r w ) 2 2 2 2 2 ' 8a (a +w lim w((a + w + 7 r ft—» + 00 J Q rh 2 2 2 + 2 2 2 ) 2 7r ) 2 (A.6) 2 2 + 2^ a; ) 2 2 2 (A. 7) 2 From Equation (A.7), since the integrand of the integral is of the form ^, L'Hospital's rule can be applied. Equation (A.7) becomes: C * A . 1.2 = lim f[ £ [ 0 ijj{x) = e a x 8 a I£[w((a h-*+ooJ V + ^ + ^) l du < + w + 7r ) + 27r w ) 2 2 2 2 2 2 2 oo (A.8) cos(irx) From the Fourier Transform table, 2 0 x 2 2 u (A.9) where T(.) is the Fourier operator. By duality, _ Tie By setting a = x2 i 2^COS(TTX)) <r (uj-7r) 2 = av27re 2 2 , (A. 10) the Fourier transform of ip(x) becomes: , T f , s /7T V a _(u,-„) 2 (A.ll) For a — 7 which is used as a fixed value in the parametric model established in Chapter 3, Equation (A.ll) can be written as: Appendix A. 95 7T (»-*) 2 (A. 12) 28 By using Equation (A.l), the admissibility function is: ,2 |2 (A. 13) Equation (A. 13) shows that the mother wavelet function is not strictly admissible. In fact, this function is similar to the Gaussian wavelet with center frequency of w designed by 0 Morlet which is given by: [3] h(x) — 7r e -JUQX ~\ (A.14) 4 The Morlet wavelet is not strictly admissible since the DC value is not zero. However, the function can be regarded as practically admissible if the DC value is small enough or close l to zero [45]. The condition can be applied to Equation (A. 12) where the DC value is: (A.15) Therefore, the mother wavelet used in Response Model #2 is said to be practically admissible. A.2 Computation of MD Performance Index Suppose that a scan of raw measurements is denoted as y {t). raw The MD component ymd{t) for that scan of measurements is defined as the mean value of the measurements denoted as y {t). raw For computing the MD performance indices, a sliding window of 80 elements is used with 50 minimum elements. The performance index is arbitrarily set to 1 (i.e. M — 1) Appendix A. 96 until 50 scan averages are accumulated. The variance of elements contained in the sliding window is denoted as a^. The discrete linear model that is used consists of the plant dynamics term and the disturbance term which is given by: » -|P«^ " ^ 5 5 " W (,) ( + ( A 1 6 ) where A (q~ ) and B(q^) are both polynomials of degree n and m respectively with m <n. l x q~ is the delay operator defined as q~ u{t) = u(t — 1). The input e{t) is a sequence of l l uncorrected random variables with a zero mean and constant variance, N(0,a ). The 2 operator A is defined as A = 1 — q~ and the term A allows the mean of the disturbance l d output exhibiting nonstationary behaviour [5]. The control signal which results in the minimum achievable variance in the output is [36]: The minimum variance control law and polynomials F(q~ ) and G(<7 ) can be found by x _1 solving the Diophantine equation: Ciq- ) = F( - )A (, - )A + q^Giq- ). 1 1 (? 1 2 (A.18) 1 d 3 As a result, the process output minimum achievable variance is calculated as the following: a 2 = E[y (t + k)} = E[F(q- )e(t + k)} (A.19) = ^(l +/i+/ (A.20) 2 mv l 2 2 2 + '-- + A -i) 2 The algorithm for computing the MD performance index for a scan of data is based on the use of the Laguerre network and Recursive Extended Least Square (RELS). The goal is Appendix A. 97 to compute the minimum estimated achievable output variance suggested by [5]. In this thesis, only the necessary equations describing the Laguerre network and the RELS are given. The properties of the Laguerre network can be found in [58, 6, 52] while the details of the RELS can be found in [38]. Any stable sampled linear process h(t) can be represented by an infinite expansion of Laguerre functions li(t) [5]: oo Ht) = J29Mt), i=i (A.21) where & is the ith Laguerre gain and kit) is the output of the ith Laguerre filter which is defined as [58]: where a is the time scale of the filter and q is the discrete time operator. The calculation of a can be found in [58]. In practice, the infinite expansion stated in Equation (A.21) can be truncated to N summing elements as: N • = i=l (A.23) In this case, the expansion is truncated to having 10 Laguerre filters (N = 10). The discrete Laguerre network in Equation (A.21) can be represented in state-space form: l{t + l) = Al{t) + be{t) (A.24) y(t) = c l(t) + e(t) (A.25) T where elements of the A matrix are given by: Appendix A. 98 Clij =a i=3 Clij =0 i <j Clij = (- -ay-'-^l-a ) (A.26) i>j 2 The b matrix elements are defined as: bi = ( - a ^ v T ^ o * i =1 (A.27) N The c and I vectors are respectively defined as: 9i l(t) = [l (t) l 92 (A.28) 9N (A.29) l (t) ••• l (t) P 2 N As a result, the minimum achievable output variance can be found by using Equation (A.20) and is given by: (A.30) crL = <\\ + (c b) + (c Abf + --- + {c A ~ b)% 1 2 1 1 k 2 where the terms c b, c Ab, c A b,... are the Markov parameters of the process. It can be T T T 2 seen from Equation (A.25) that the white noise sequence e(t) is unknown. The recursive extended least squares (RELS) [38] can then be applied to estimate c {t) and a . T 2 e l{t) = Al(t-l) + bn(t-l) (A.31) c{t) = c{t-l) + P{t)l{t)[y(t)-c {t-l)l{t)} (A.33) n(t) = y(t)-c (t)l(t) (A.34) T T The white noise e(t) can be estimated using the residual rj(t). Therefore, the estimated Appendix A. 99 minimum achievable output variance is: at, = $[1 + (c bf + {<?Ab? + ••• + (c A -*bn T T (A.35) k Finally, by using the previous equation and the MD variance of the sliding window, the MD performance index is calculated as: Af=^. A.3 (A.36) mv ^ Computation of Power Spectral Densities of CD Profiles The objective for estimating power spectrum of a finite signal is to observe the distribution over frequency of the power contained in a signal. The following is the summary of the material found in the documentation of the Signal Processing Toolbox in Matlab [39]. The power spectrum S (u) of a sequence which is, in this case, CD profiles y(x) is the xx discrete-time Fourier transform of the autocorrelation of the sequence known as Rx {m): X + OC S (u) = Y, Urn)e- " R xx j m (A.37) 771—— OC The power spectral density (PSD) of the CD profile sequences can be expressed as: PM = (A.38) The unit of the PSD is watts/rad/sample or simply watts/rad. In this thesis, the estimation of the power spectral density is achieved by using Welch's method. For details about this method, refer to [43, 26, 47]. The following is the procedure for computing the PSD using Welch's method for an input vector y [39]: Appendix A. 100 1. The vector y is first segmented into eight sections with equal length, each with 50% overlap 2. Hanning window with the same length as the each dividing section is applied to each segment of y 3. The sequence is then zero padded to the nearest power of 2. Fast Fourier Transform (FFT) is then applied to the windowed, zero padded data 4. The periodogram (power spectral densities) of each windowed segment is computed according to the mathematics described above 5. A set of periodograms is averaged to form the final spectrum estimate S (to) xx Bibliography [1] A.Abbate, C.M.DeCusatis, and P.K.Das. "Wavelets and Subbands". Birkhauser, Boston, 2002. [2] A.Cohen, I.Daubechies, and J.C.Feauveau. "Biorthogonal Bases of Compactly Supported Wavelets". Comm. Pure Appl. Math., v.45:pp.485 - 560, 1992. [3] A.Grossman and J.Morlet. "Decomposition of Hardy Functions into Square Integrable Wavelets of Constant Shape". SIAM Journal on Mathematical Analysis, v. 15 (4): pp. 723 - 736, 1984. [4] B.Jawerth and W.Sweldens. "An Overview of Wavelet Based Multiresolution Analyses". Technical report, Technical Report 1993:1, Industrial Mathmematics Initiative, Department of Mathematics, University of South Carolina. [5] C.B.Lynch and G.A.Dumont. "Control Loop Performance Monitoring". IEEE Transactions on Control Systems Technology, 4(2):pp.l85 - 192, March 1996. [6] C.C.Zervos and G.A.Dumont. "Deterministic Adaptive Control Based on Laguerre Series Representation". Int. J. Contr., 48(6):pp.2333 - 2359, 1988. [7] C.S.Burrus, R.A.Gopinath, and H.Guo. Introduction to Wavelets and Wavelet Trans- forms". Prentice Hall, Upper Saddle River, New Jersey 07458, 1998. [8] D.Gabor. "Theory of Communication". Inst. Elec. Eng., v.93:pp.429 - 457, 1946. [9] D.Gorinevsky, M.Heaven, C.Hagart-Alexander, M.Kean, and S.Morgan. "New Algorithm for Intelligent Identification of Paper Alignment and Non-linear Shrinkage". Pulp and Paper Canada, pages 76-81, June 1998. [10] D.L.Donoho. "Nonlinear Wavelet Methods for Recovery of Signals, Densities and Spectra from Indirect and Noisy Data". In Proceedings of Symposia in Applied Mathematics, volume 00, pages 173-205. American Mathematical Society, 1993. 101 Bibliography 102 [11] D.L.Donoho. "De-noising by Soft-thresholding". IEEE transactions on Information Theory, 41(3):pp.613-627, May 1995. [12] D.L.Donoho and I.M.Johnstone. "Miniraax Estimation via Wavelet Shrinkage". Technical report, Technical Report 402, Stanford University, Department of Statistics, July 1992. [13] D.L.Donoho and I.M.Johnstone. "Adapting to Unknown Smoothness via Wavelet Shrinkage". Technical report, Technical Report 405, Stanford University, Department of Statistics, June 1993. [14] D.L.Donoho and I.M.Johnstone. "Ideal Spatial Adaptation via Wavelet Shrinkage". Biometnka, 81:pp.425-455, July 1994. [15] E.M.Heaven. "Application of System Identification to Paper Machine Model Development and Simulation". Pulp and Paper Canada, pages pp.49 - 54, April 1996. [16] OPC Foundation. "About OPC - What Is OPC?". [17] OPC Foundation. "OPC Technical Overview". [18] OPC Foundation. "OPC Data Access Custom Interface Specification 3.0", March 4, 2003. [19] G.A.Dumont, J.Ball, G.Emmond, and M.Guillemette. "Paper Machine CD Control using Wavelet Based MD/CD Separation". In Advanced Process Control Application for Industry Workshop. IEEE Industry Applications Society, 2003. [20] G.A.Smook. "Handbook for Pulp and Paper Technologists". Angus Wilde Publications Inc, Vancouver, 2nd edition, 1992. [21] G.E.Stewart. "Two-Dimensional Loop Shaping: Controller Design for Paper Machine Cross-Directional Processes". PhD thesis, The University of British Columbia, July 1997. Bibliography 103 [22] G.E.Stewart, J.U.Backstrom, P.Baker, C.Gheorghe, and RN.Vyse. "Controllability in Cross Directional Process". Pulp and Paper Canada, v.l03(8):pp.32 - 38, August 2002. [23] Dimitry Gorinevsky and Cristian Gheorghe. " Identification Tool for Cross-Directional Processes". IEEE Transactions on Control Systems Technology, v. 11 (5), September 2003. [24] G.Strang. "Wavelets". American Scientist, v.82:pp.250 - 255, 1994. [25] G.Strang and T.Nguyen. "Wavelets and Filter Banks". Wellesley Cambridge Press, Cambridge, 1996. [26] M. Hayes. "Statistical Digital Signal Processing and Modeling". John Wiley and Sons, 1996. [27] I.Daubechies. "Orthonormal Bases of Wavelets with Finite Support - Connection with Discrete Filters". In J.M.Combes, A.Grossman, and P.Tchamitchian, editors, Wavelets, Time-Frequency Methods and Phase Space, Marseille, France, December 1987. Proceedings of International Colloquium on Wavelets and Applications. [28] I.Daubechies. "Orthonormal Bases of Compactly Supported Wavelets". Comm. Pure Applied Math., XLI(41):pp.909 - 996, November 1988. [29] I.Daubechies. "The Wavelet Transform, Time-Frequency Localization and Signal Analysis". IEEE Trans. Information Theory, v.36(5):pp.961 - 1004, 1990. [30] I.Daubechies. "Ten Lectures on Wavelets". Lowell, MA, 1992. CBMS-NSF Conference on Wavelets and Applications, SIAM, Philadelphia, PA. [31] Mitsubishi Heavy Industries. "ME Headbox-C IV Operation and Maintenance Man- ual", March 2001. [32] J.Ball, January 14, 2004. Email Conversation. Bibliography 104 [33] J.Ghofraniha. "Cross Directional Response Modeling, Identification and Control of Dry Weight Profile on Paper Machine". PhD thesis, The University of British Columbia, July 1997. [34] X. Jiao. "Paper Machine Data Analysis and Optimization using Wavelets". Master's thesis, The University of British Columbia, January 1999. [35] K.Cutshall. "Pulp and Paper Manufacture", volume 7, chapter Paper Machine Operations - Cross-direction Control, pages 472 - 506. 3rd edition, 1991. [36] K.J.Astrom. Introduction to Stochastic Control Theory". Academic, New York. [37] K.Kristinsson. "Cross Directional Control of Basis Weight on Paper Machine using Gram Polynomials". PhD thesis, The University of British Columbia, January 1994. [38] L.Ljung and T.Soderstrom. Theory and Practice of Recursive Identification". MIT '. Press, Cambridge, MA. [39] The Mathworks. "Signal Processing Toolbox User's Guide". [40] M.V.Wickerhauser. "Smooth Localized Orthonormal Bases". Comptes Rendus de I'Academie des Sciences de Paris, 316:pp.423-427, 1993. [41] O.Rioul and M.Vetterli. "Wavelet and Signal Processing". IEEE Signal Processing Magazine, 8(4):pp.l4 - 38, October 1991. [42] P.C.Lemarie and Y.Meyer. "Ondelettes et bases Hilbertiennes". Rev. Mat. IberoAmer, pages pp.1 - 18, 1986. [43] P.D.Welch. "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms". IEEE Trans. Audio Electroacoustics, v.AU15:pp.70 - 73, June 1967. [44] P.E.Wellstead, W.P.Heath, and A.P.Kjaer. "Identification and Control for Web Forming Processes". In 13th World IFAC Congress, volume 7, San Francisco, CA. Bibliography 105 [45] P.Flandrin. "Time-Frequency/Time-Scale Analysis". Academic Press. USA, 1999. [46] P.Heller and R.Wells. "Sobolev Regularity for Rank M Wavelets". Technical report, CML Technical Report TR 96-08, Rice Univerisity, Computational Math. Laboratory, 1996. [47] PStoica and R.L.Moses. "Introduction to Spectral Analysis". Prentice-Hall, Englewood Cliffs, NJ, 1997. [48] R.D.Braatz, A.P.Featherstone, and B.A.Ogunnaike. "Identification, Estimation and Control for Sheet and Film Processes". In 13th World IFAC Congress, volume 7, San Francisco, CA. [49] S.C.Chen. Actuation Cell Response and Mapping Determinations for Web Forming Machines, 1992. US Patent 5,122,963. [50] S.C.Chen. "Full-width Sheet Property Estimation from Scanning Measurements". In Control Systems '92, pages pp.123-130, Whistler, B.C., Canada, September 1992. [51] S.C.Chen and R.Subbarayan. "Identifying Temporal and Spatial Responses of Cross Machine Actuators for Sheet-forming Processes". IEE Proc. - Control Theory Appi, 149(5):pp.441 - 447, September 2002. [52] S.Gunnarsson and B.Wahlberg. "Some Asymptotic Results in Recursive Identification using Laguerre Models". Int. J. Adaptive Contr. and Signal Processing, 5:pp.313 333, 1991. [53] S.Mallat. "A Theory of Multiresolution Signal Decomposition: The Wavelet Representation". IEEE transactions on pattern analysis and machine intelligence, ll:pp.674- 693, 1989. [54] S.Mallat. "A Wavelet Tour of Signal Processing". Academic Press, San Diego, 2nd edition, 2001. Bibliography 106 [55] S.Palavajjhala, R.L.Motard, and B.Joseph. "Wavelet Applications in Chemical Engi- neering", chapter Computational Aspects of Wavelets aned Wavelet Transform. Ed. Motard and Joseph, Kluwer Academic Publishers. [56] S.R.Duncan. "Estimating the Response of Actuators in a Cross-Directional Control System". In Control Systems '96, pages pp.19-22, Halifax, N.S., Canada, May 1996. [57] IPCOS Technology. "OPC For Matlab User Manual 3.0.0.3". [58] Y.Fu and G.A.Dumont. "Optimum Laguerre Time Scale and its Online Estimation". IEEE Transactions on Automatic Control, 38(6):pp.934 - 938, June 1993. [59] Y.Meyer. "Ondelettes. Filtres Miroirs En Quadrature Et Traitement Numerique De L'image". Technical report, Technical Report, Cahiers de Mathematiques de la Decision, Paris, France, 1989. [60] Y.Shimanuki. "OLE for Process Control (OPC) for New Industrial Automation Systems", volume 6, pages pp.1048 - 1050. IEEE SMC '99 Conference Proceedings, IEEE International Conference on Systems, Man, and Cybernetics, 1999. [61] Z.Nesic. "Paper Machine Data Analysis using Wavelets". Master's thesis, The University of British Columbia, February 1996. [62] Z.Nesic, M.S.Davies, and G.A. Dumont. "Paper Machine Data Analysis and Compression Using Wavelets". TAPPI Journal, v.80(10):pp.l91 - 204, 1997. [63] Z.Nesic, M.S.Davies, G.A. Dumont, and D.Brewster. "CD Control Diagnostics Using a Wavelet Toolbox". International CD Symposium'97, Finland, June 1997.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Real time paper machine data wavelet analysis and cross...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Real time paper machine data wavelet analysis and cross directional actuator response identification Kan, Ka Keung 2004
pdf
Page Metadata
Item Metadata
Title | Real time paper machine data wavelet analysis and cross directional actuator response identification |
Creator |
Kan, Ka Keung |
Date Issued | 2004 |
Description | The development of an online automated diagnostics tool for analyzing paper machine data is useful to mill operators when monitoring the performance of a cross directional (CD) control system. The use of Object linking and embedding for Process Control (OPC) interface as a medium for client-server based data acquisition provides an efficient, robust and accurate data transfer between a CD control system and a monitoring workstation. Therefore, a real time diagnostics toolbox has been developed to help paper mill companies reduce costs and save time on troubleshooting. The Real Time Wavelet Toolbox (known as Waveplot and developed over several years at UBC, [62]) provides color maps of raw, wavelet and residue profiles for both basis weight and moisture. Moreover, machine directional (MD) and cross directional (CD) control performance indices as well as the standard deviations for controllable and non-controllable CD components are calculated and displayed. Other important plots such as power spectral density and two-sigma plots are also included in the toolbox for further analysis. Waveplot also helps mill operators identify the high-resolution CD actuator response shape and CD mapping through bump tests. The identification is achieved by using a continuous wavelet transform with the chosen CD actuator response shape model as the mother wavelet. The control performance index reflects recent views on CD performance index in terms of CD controllability given the knowledge of the identified CD actuator response shape. Furthermore, the identified actuator response centers could help mill operators investigate any problems due to CD mapping misalignment and shrinkage, particularly at the paper edges. The displays and identification results are validated and demonstrated by collecting paper profiles through the OPC server from a multi-grade, operating paper machine. |
Extent | 17915701 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065597 |
URI | http://hdl.handle.net/2429/15269 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2004-0242.pdf [ 17.09MB ]
- Metadata
- JSON: 831-1.0065597.json
- JSON-LD: 831-1.0065597-ld.json
- RDF/XML (Pretty): 831-1.0065597-rdf.xml
- RDF/JSON: 831-1.0065597-rdf.json
- Turtle: 831-1.0065597-turtle.txt
- N-Triples: 831-1.0065597-rdf-ntriples.txt
- Original Record: 831-1.0065597-source.json
- Full Text
- 831-1.0065597-fulltext.txt
- Citation
- 831-1.0065597.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065597/manifest