Online Detection of Picketing and Estimation of Cross Direction and Machine Direction Variations using the Discrete Cosine Transform by Soroush Karimzadeh B.Sc., Sharif University of Technology, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Electrical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December, 2008 © Soroush Karimzadeh 2008 Abstract A scanning sensor at the dry end of a paper machine samples a noisy mix ture of both cross direction and machine direction variations. Although many techniques have been developed to separate these two variations and estimate a filtered profile, industrial practice has changed little over the years. In this work, a novel online Cross Direction (CD) /Machine Direction (MD) separation approach is developed. The method uses a fundamental prop erty of the CD variations to separate them from the MD variations. It is shown that the scanned CD variations build an even periodic function in the time domain and appear only as integer multiples of twice the scanning frequency in the frequency domain. Based on this property, the Discrete Cosine Transform (DCT) is utilized to separate the CD profile from noise and the MD variations. The performance of this method is verified through comprehensive simulation tests and is compared with existing wavelet and exponential filters. It is shown that the suggested method has better per formance in simulation scenarios. This is further validated by applying the method to industrial data. In the second part of this thesis, a novel picketing detection method is devel oped. The picketing pattern is associated with growing components in the spectrum of actuator profiles. A conventional change detection algorithm, CUSUM, is applied to detect these growing components. The method is verified in two simulation scenarios. It is shown that this method can detect and predict a picketing pattern before any picketing pattern is visible on the raw profile. 11 Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures . . . . . . vi Acknowledgements viii Dedication ix 1 Introduction and Background 1 1.1 Review of Existing MD/CD Separation Methods 3 1.2 Motivation 4 1.3 Outline of Thesis 4 2 Fundamental Properties of Paper Machine Variations . . 6 2.1 General Assumptions about CD and MD Variations 6 2.2 The Fundamental Property of CD Variations 7 3 The Discrete Cosine Transform 11 3.1 Definition of DCT 11 3.2 Relationship between DFT and DCT 13 3.3 CD/MD Separation Using DCT 13 4 Applying the DCT Method to Simulated and Industrial Data 16 4.1 Wavelet and Exponential Filtering 16 4.2 Edge Effects and Step Changes 18 4.3 Smooth Surfaces 23 4.4 Grade Profile Change 28 4.5 Data Created by Moisture Variation Model 33 111 Table of Contents 4.6 Bumps and Streaks 39 4.7 Spanning the MD Frequencies 44 4.8 Comparison of the Estimated Profiles Using Industrial Data 46 4.9 Conclusion 48 5 A Picketing Detection Algorithm 51 5.1 The Idea 51 5.2 CUSUM Algorithm 53 5.3 Simulation Results 54 5.3.1 Misalignment Scenario 54 5.3.2 No Spatial Filtering Scenario 54 5.4 Conclusion 58 6 Conclusions and Future Work 62 6.1 Conclusions 62 6.2 Future Work 63 Bibliography 64 7 Appendix A - Matlab Programs for Simulations in Chapter 4 67 7.1 The Main Program 67 7.2 The DCT Program 77 8 Appendix B - Matlab Programs for Simulations in Chapter 5 81 iv List of Tables 4.1 The J error for different methods in the edge effects scenario 20 4.2 The J error for different methods in the smooth surfaces sce nario 24 4.3 The J error for different methods in the grade change scenario 30 4.4 The J error for different methods in the moisture variation model scenario 35 4.5 The J error for different methods in the bumps and streaks scenario 41 V List of Figures 1.1 Simplified paper machine (J. Ghofraniha) [1] 1 1.2 Scanner path on the paper 2 2.1 A sample time invariant CD profile 7 2.2 A sample profile with no MD variation showing the scanner path on the profile 8 2.3 The scanned data versus time 8 2.4 The Fourier series coefficients a and b for the signal in Fig ure2.3 9 3.1 The four different forms of DCT 12 4.1 Edge Effects: raw scan without noise 18 4.2 Edge Effects: raw scan with noise 19 4.3 Edge Effects: estimated profiles 21 4.4 Edge Effects: estimated profiles 22 4.5 Edge Effects: estimated error versus time (scan number) . 23 4.6 Edge Effects: estimated MD profiles 24 4.7 Smooth Surfaces: raw scan without noise 25 4.8 Smooth Surfaces: raw scan with noise 25 4.9 Smooth Surfaces: estimated profiles 26 4.10 Smooth Surfaces: estimated profiles 27 4.11 Smooth Surfaces: estimated error versus time (scan number) 28 4.12 Smooth Surfaces: estimated MD profiles 29 4.13 Grade Profile Change: raw scan without noise 29 4.14 Grade Profile Change: raw scan with noise 30 4.15 Grade Profile Change: estimated profiles 31 4.16 Grade Profile Change: estimated profiles 32 4.17 Grade Profile Change: estimated error versus time (scan num ber) 33 4.18 Grade Profile Change: estimated MD profiles 34 4.19 Moisture Variation model: raw scan without noise 35 vi List of Figures 4.20 Moisture Variation model: raw scan with noise 36 4.21 Moisture Variation model: estimated profiles 37 4.22 Moisture Variation model: estimated profiles 38 4.23 Moisture Variation model: estimated error versus time (scan number) 39 4.24 Bumps and Streaks: raw scan without noise 40 4.25 Bumps and Streaks: raw scan with noise 41 4.26 Bumps and Streaks: estimated profiles 42 4.27 Bumps and Streaks: estimated profiles 43 4.28 Bumps and Streaks: estimated error versus time (scan number) 44 4.29 Bumps and Streaks: estimated MD profiles 45 4.30 The target CD profile 45 4.31 The raw scan without noise 46 4.32 The raw scan with noise 47 4.33 The CD estimation Error for different MD frequencies . . 47 4.34 Raw Scans 48 4.35 Estimated profile using exponential filtering 49 4.36 Estimated profile using wavelet 49 4.37 Estimated profile using DCT 50 5.1 A raw profile which has picketing 525.2 The actuator profile corresponding to the raw scan 52 5.3 The raw profile 55 5.4 The actuator profile 55 5.5 The power of highest frequency component versus scan num ber 56 5.6 The g values for P.j versus scan number 562Xa5.7 P. versus scan number 57 5.8 The g values for P. versus scan number 57 5.9 The raw profile 58 5.10 The actuator profile 59 5.11 The power of highest frequency component versus scan num ber 59 5.12 The g values for Pj versus scan number 602Xa5.13 P. versus scan number 60 5.14 The g values for Pr versus scan number 61 vii Acknowledgements My deepest gratitude to my thesis supervisor Dr.Guy Dumont for giving me the opportunity to make this accomplishment and guiding me through the path of our research. I would also like to thank the control group at Pulp and Paper cen ter, Drs Michael Davies, Philip Loewen, Bhushan Gopaluni and Jun Yan for their insightful discussions and advice. My thanks and best wishes to my colleagues and friends Mohammad Ammar, Fazel Farahmand and Au Soltanzadeh who accompanied me during this journey. My sincere thanks to Hassan Masoom for reviewing this thesis. The industrial results presented in this thesis were produced in Canfor Paper mill in Prince George. Special Thanks to Grant Emmmond, Michel Guillemette and Alexandre Rousseau. My gratitude to Lisa Hudson, Ken Wong and Doris Metcalf for their assistance whenever I encountered a problem. This work was performed as a collaboration between Canfor Pulp and Paper, Honeywell Process Solutions and the University of British Columbia. viii Dedication This thesis is dedicated to my family, Sanaz Karimzadeh, Majid Karimzadeh and Farimah Arzpayma, for their continuous support and their everlastinglove. ix Chapter 1 Introduction and Background The final stage in the production of paper is the paper machine. Here the stock, a suspension of fibre in water, is deposited on a web that carries the stock over a number of drainage elements such as foils and blades and through a number of pressing rolls and driers to dry the paper ([2]). A simple schematic for a paper machine is shown in Figure 1.1. Figure 1.1: Simplified paper machine (J. Ghofraniha) [1] S 1t zcw D iy )±e secLb n C aZnders tack 0 —fr3a Scanner R eel 1 Chapter 1. Introduction and Background The closed loop control of important paper sheet properties is carried out by two controllers which are known as the cross direction (CD) controller and the machine direction (MD) controller. The cross direction controller is a regulator which minimizes the variations across the paper, while the machine direction controller tracks a setpoint. The CD basis weight profile is controlled by a number of slice lip or dilu tion actuators in the headbox. The CD moisture content profile is controlled by rewet showers, steam boxes, or infrared heating boxes. The feedback is completed by a scanner which scans across the paper to measure its basis weight and moisture. The measured profile is then sent back to the corre sponding controller to generate the new instructions for the actuators. The controllers need to be fed by accurate measurements to be able to control the output profile correctly. This, however, has posed a challenge ever since the 0-frame scanners were built [3]. The challenge is that since it takes between 20 and 30 seconds for the scanner to traverse the sheet and the paper is produced at a speed of upto 120 km per hour, the scanner only scans the web along a zigzag path on the paper (Figure 1.2). This implies neither the true CD variation nor the true MD variation is available for control. This problem has been referred to as MD/CD separation in the paper production literature. A similar issue is present in other sheet and film processes as well [4]. Figure 1.2: Scanner path on the paper Different solutions have been developed over the years to address this problem. Most of the mechanical solutions are considered to be costly, com plex and impractical [3]. Thus, software solutions have been developed to Machine Direction 2 1.1. Review of Existing MD/CD Separation Methods avoid expensive mechanical redesign. This means a data processing block isneeded to estimate the CD and MD variations before they are sent to thecorresponding controllers. 1.1 Review of Existing MD/CD SeparationMethods The conventional method used in industry to separate MD and CD variations is an exponential filter which simply uses a first order filter in themachine direction to estimate the CD profile. The mean of each scan isthen extracted as the MD estimate [5]. This method, however, does notconsider the scanning geometry and is not based on the properties of thevariations. Chen suggested the use of a dual Kalman filter to estimate both the CDand MD variations using a linear model([6]). Using Lindeborg’s nonlinearmodel (see [7]), the EIMC algorithm discussed in [8] combines a Kalman filter to estimate the CD moisture profile and a modified least squares parameter identifier to estimate the MD variation. Several machine and operatingpoint dependent parameters must be tuned for this method which makesits online implementation challenging. In [9] an ARMA MD model is usedto estimate the MD basis weight variations with an identification methodsimilar to the one used in [8]. In [10] the discrete wavelet transform was applied for estimation andbetter visualization of MD and CD variations. Furthermore, significantcompression of the process data was achieved. The wavelet properties wereused in [11] and [12] to compress and store the data. In [13] it was shownthat the wavelet transform can compress the process data efficiently. Inaddition wavelet and wavelet packets were used to separate the controllableand uncontrollable variations. The main application of wavelets to separate MD and CD variations waspublished in [14] and [15]. There, it was shown that wavelets can achievebetter estimates for CD variations compared to conventional methods. Also,in this work, the focus has been towards filtering the profile such that onlythe controllable frequencies are passed to the controller. The method owesits good performance to the strong denoising ability of wavelet filters in twodimensions, but still did not utilize the scanning geometry.In [16] and [17] a method is suggested that uses wavelet filtering andexponential filtering in a bootstrap recursive scheme to separate the MDand CD variations. By using a variable scan rate in [18], a method to 3 1.2. Motivation identify the aliased components of MD variations from the estimated CDprofile was suggested. But no systematic way to exploit this informationhas been developed so far. The more recent work has been increasingly focused on the scanninggeometry. Specifically, Duncan and Welistead in [19] and [20] used thegeneralized sampling theorem to estimate the MD variations and separatethem from CD variations. This work is based on the fact that at differentCD locations, the MD signal is being sampled at different rates. Thus,by using reconstruction filters with weighting coefficients dependent on theCD position, they were able to reconstruct the true MD variation. This,however, was only achieved for offline data and can not be implementedonline since the filters are noncausal. Another recent work by Duncan and et al. [21] utilizes one scanninggauge and one fixed gauge. This solution is not desirable in many millsas the method operates only at the cost of a new scanner, which is a veryexpensive proposition. 1.2 Motivation In this work, a novel method for CD! MD separation is proposed and evaluated based on the scanning geometry. The scanning geometry creates atypically overlooked property in the CD profile, which can be used to separate the CD profile from MD variations. The simplicity of the approach inthis work that exploits this property allows an online implementation withonly a software upgrade. 1.3 Outline of Thesis The fundamental properties and assumptions about CD and MD variationsare presented in chapter 2. In addition, a new perspective into propertiesof CD profiles is presented as the main building block for a new CD/MDseparation technique. In chapter 3, the Discrete Cosine Transform is introduced. Some of theproperties and similarities to the discrete fourier transform are investigated.Based on the properties of the DCT, a novel online CD/MD separationtechnique is presented in section 3.3. Through various simulation scenarios in chapter 4, the performance ofthis technique is compared with the existing CD/MD separation methods. 4 1.3. Outline of Thesis The method is also verified by applying it to real measurements from a papermachine. In the second part of this thesis, chapter 5, a novel online picketingdetection method is introduced. The performance of this method is alsoevaluated though different simulation scenarios. 5 Chapter 2 Fundamental Properties of Paper Machine Variations In this chapter, a fundamental property of the scanned data is explored.First, the main assumptions about CD and MD variations are discussed.Then, a lesser known property of CD variations, which form the basis of theapproach taken in this work, is explained. 2.1 General Assumptions about CD and MD Variations In this work it is assumed that the total scanned profile is p(x,t) Pcd(X,t) +Pmd(t) (2.1) r(x, t) = p(x, t) + w(x, t) (2.2) Here Pmd(t) denotes the true MD variation, pcd(x, t) denotes the crossdirection variation, p(x, t) the exact sum of the two variations with no noise and w(x, t) represents the measurement noise present in the measured 2D sigiial r(x, t). Since the MD variations are created by disturbances prior tothe headbox and before the distribution of stock on the sheet, MD variations are assumed to be independent of the CD position. This implies that Pmddoes not depend on position x. It is also reasonable to assume that the crossdirectional variation, which is created by the flow patterns inside the head-box, changes much slower with time in comparison to the MD variations.Thus, for piece-wise analysis, we can safely assume pcd(X, t) Pcd(X) whichimplies that the CD profile does not depend on time. Another important observation from the above model is that the esti mated profiles should be filtered to eliminate the noise n(x, t). In case ofCD profile, the estimated profile should also be filtered so that only the controllable variations are passed on to the CD controller. 6 2.2. The Fundamental Property of CD Variations 2.2 The Fundamental Property of CD Variations In this section the unique effect of CD variations on the power spectrum of data is discussed. Assume a sheet of width 50, in which there are noMD variations. A sample CD profile is shown in Figure 2.1. The scanner measures the profile along the path shown in Figure 2.2. Figure 2.1: A sample time invariant CD profile If the data is graphed, as it was measured by the scanner, Figure 2.3is the result. Here, the sampling rate is 1 sample per second. Thereforeit takes T8 = 50 seconds to scan the width of the sheet (T5 denotes the scanning period). Figure 2.3 shows that the scanned data is an even periodic signal. Scans one and three are equal to each other. Scans two and four are also equal.In addition, scans two and four are equal to the reflections of scans one and three. Alternatively, this can be seen as scans one and two being repeatedperiodically with period 2T3. A possible problem is the zeros padded on both sides of the measured scan in the signal coming from the scanner. This is due to the fact that the scanner goes off the sheet after it reaches the edge of paper before comingback. As a result, the measurements made in this interval are recorded asNAN ( not a number ) and later replaced by zeros. Also, the number of zeros can change with time but this happens very slowly. Interestingly, this fact does not affect the periodic and even property of the CD profile because equal number of zeros are padded on the side of adjacent scans. CD position 7 2.2. The Fundamental Property of CD Variations Figure 2.2: A sample profile with no MD variation showing the scanner path on the profile Figure 2.3: The scanned data versus time 1 scan 2 scan 3 scan 4 Time 8 2.2. The Fundamental Property of CD Variations Using the Fourier series ([22]) any periodic signal f(x) can be writtenas: f(x) = ao + [a cos(nwox) + b sin(nwox)J Also, it is known that for an even periodic function, the Fourier series hasonly cosine terms and the coefficients for all sine functions are zero. Thismeans, in absence of any MD variation, the scanned signal can be expressedas a sum of cosine functions with frequencies equal to integer multiplesof (2T8)’. The Fourier series coefficients for the signal in Figure 2.3 areshown in Figure 2.4. The Fourier sine coefficients are all zero. The Fouriercosine coefficients are also zero for every second coefficient. This is due tothe internal even symmetry of one period of the signal in Figure 2.3. Thisproperty can be exploited to extract the CD variations spectrum from thescanned signal. 200Cr’ ISOH 500 IS 30 40 & 10 10 2030 40 51 0.012000 1’ I 21 1500 ooos .21 1000 500 a I 00 5 10 15 5 10 15n Figure 2.4: The Fourier series coefficients a and b for the signal in Figure2.3 Relaxing the assumption about the MD variations, it is noticed that theMD variations can only interfere with the CD profiles spectrum if they havecomponents close to multiple integers of - . This effect will be simulatedin section 4.7. The fundamental property discussed in this section will be exploited to sep 9 2.2. The Fundamental Property of CD Variations arate the CD profile from the MD variation in the next chapters. For more information about this property interested reader is referred to [23]. 10 Chapter 3 The Discrete Cosine Transform In this chapter, the Discrete Cosine Transform (DCT) is first defined. In section 3.2, the relationship between the Discrete Fourier Transform and DCT is investigated. The main contribution of this work, a new algorithm for CD/MD separation, is presented in 3.3. 3.1 Definition of DCT Assuming an arbitrary N point signal x [n], there are eight different ways to create an even periodic sequence from it. The four more common ones are illustrated in Figure 3.1. In this Figure, the x[n] signal has a length N = 4 corresponding to the first four points in the graph. ij [nJ is periodic with period 2N — 2 = 6 about n = 0 while 2[n] is periodic with 2N = 8 about the half sample point n —. 3[n] has period 4N = 16 and is symmetric about n = 0 and x4[n] has the same period but its symmetric around the half sample point ii = — These four different cases show the implicit periodicity in the 4 well known forms of DCT known as DCT-1, DCT-2, DCT-3 and DCT-4. Among these four forms, DCT-1 and DCT-2 are most used in practice. It is also worth mentioning that there are 8 ways to build odd peri odic functions from x[n] which make the building blocks for Discrete Sine Transform (DST). The Type-i DCT is defined by the forward/inverse transform pair ([24] and [25]): xcl[k] = 2a[nx[n]cos(1), k 0 k N —1 (3.1) x[n] = N — 1 [k]Xd1[k]cos(T), 0 n N —1 (3.2) ii 3.1. Definition of DCT ( 5 10 15 n (a) w Icn (1) X4[F1] :1Ll -4iflI ii []{2, n=OandN—1 1, OnN—2 N-i rk(2n-j- 1)Xc2[kj = 2 x[n]cos( 2N n=O N-i1 irk(2n+1)x[n]=/3[k]Xc[k]cos( 2N k=O 4 I ci I I li T 11111 IT I fiJi 2I1111 (c) I it (d) Figure 3.1: The four different forms of DCT where: whereas the Type-2 DCT is: OkN—1 (3.3) OnN—1 (3.4) where: k = 11, 1kN—1 The cos(-) for DCT-1 arid cos(7) for DCT-2 form a collection of orthogonal basis functions. This can be proven by showing that these functions are the eigenvectors for a real symmetric matrix (see [26]). 12 3.2. Relationship between DFT and DOT 3.2 Relationship between DFT and DCT There is a very close relationship between various forms of the DCT and the Discrete Fourier Transform (DFT). In this section, the relationship between DCT-1 and the DFT is shown. The x [n] is one period of the periodic signal 5 [nJ in the previous section and is related to xn] by xi[n] = Xcx[((n))2N_2] + xa[((n))2N_2] = i[nj, ii = 0,1, ...,2N —3 where xe[n] = c[n]x[n]. Here the notation ((n))N is used for (n modulo N) [22]. Now if X1 [k] denotes the (2N — 2) point DFT of the (2N — 2)- point sequence xi [n] then: Xi[k]=X[kj+X[k]=2Re{Xa[k]}, k=0,1,...,2N3 where X[k] is the (2N — 2) -point DFT of the N-point sequence c[n]x[n] computed by padding N — 2 zero samples. Using the definition of DFT: X1 [k] = 2Re{X [k] } = 2 [n] x [n] cos( 2N—2 = Xc1 [k] This means that the Discrete Cosine Transform for an N-point signal x[n] is equal to the (2N — 2)-point Discrete Fourier Transform of one period of the signal, created by symmetrically extending x[nj according to Figure 3.1(a). 3.3 CD/MD Separation Using DCT In section 2.2, it was shown that, regardless of their position, if the CD variations are recorded in the time domain then the resulting signal is even and periodic. It is known that for a periodic signal with period T and wo = the Fourier series is: x(t) = ao + [Bk cos(kwot) + Ck sin(kot)J It can easily be shown that for an even periodic signal, the 0k coefficients are zero and the signal can be expressed as the sum of the DC component a0 and 13 3.3. CD/MD Separation Using DOT the cosine coefficients Bk. Since the CD variations make an even periodicsignal with period 2T3, their Fourier series expansion consists entirely ofcosine terms with frequencies equal to integer multiples of 4-. Moreover,the CD variations are considered zero mean. As a result, the DC componentwould be zero. So Pcd(t) = Bkcos() where pcd(t) denotes the measured CD profile versus time. In other words,the CD profile can be entirely represented by cosine functions at frequenciesk 2T Thus, if the raw data is expressed in terms of cosines, by keeping theterms corresponding to these frequencies and setting the rest to zero, theCD profile can be estimated. The inevitable error in this estimation will bedue to the presence of MD variations and noise at these frequencies.After estimating the CD profile, the MD variation can be estimated bysubtracting the CD profile from the raw scan. Thus, an MD estimate can befound at each sampling point. To eliminate the noise, this MD estimate hasto be filtered. In this work, the final MD estimate is set equal to the meanof MD estimates at each sampling point over each scan. However, betterestimates can be achieved if a different filter is used. Furthermore, this canincrease the MD controller’s bandwidth. The proposed separation method is as follows: Step 1 - A window consisting of the last 32 scans is updated by addingthe new scan as the last scan to it and shifting the previous scans. Step 2 - The DCT coefficients for this updated window is found. Step 3 - The coefficients with frequencies corresponding to those of the CDvariations are extracted. Step 4 - The inverse DCT of these extracted coefficients is used to reconstruct the CD profile. Step 5 - The reconstructed profile is filtered by a wavelet filter to separate the controllable and uncontrollable variations. Step 6 - The estimated CD profile is subtracted from the raw profile tofind the MD variations. 14 3.3. CD/MD Separation Using DCT Step 7 - The MD variations are filtered to remove the noise. The length of the window has been chosen as 32 because the DCT method was found to achieve better results with this window length in simulations. The performance of this method in simulation and in practice is shown in Chapter 4. 15 Chapter 4 Applying the DCT Method to Simulated and Industrial Data In this chapter, the DCT method is applied to both simulated and industrialdata. The wavelet based CD/MD separation method and the exponential filtering are explained in section 4.1. In sections 4.2 to 4.7, different simulationscenarios are utilized to evaluate the performance of the DCT algorithm. In4.8, the results from applying this method to industrial data are presented. 4.1 Wavelet and Exponential Filtering The model described by Equations 2.1 and 2.2 is continuous. In practice,however, only the sampled form of the scanned signal r(x, t) is available.Thus, the variables r(x, t), w(x, t), p(x, t) and Pmd(t) have to be replacedwith their discretized form r(n, m), w(n, m), pj(n, m) and Pmd(m), wheren denotes the CD bin and m denotes the scan number. The widely used industrial algorithm to separate MD and CD variationsis exponential filtering. In this method, a weighted sum of the previous CDestimate and the new CD scan at each CD position is considered as the newCD estimate ([5]). The average of each raw scan is used as an MD estimatefor the whole scan. The equations describing the algorithm are: pd(m) r(n,m) (4.1) pcd(n,m) = p (r(n,m) Pmd(m)) + (1— p) pcd(n,m 1) (4.2) where 16 4.1. Wavelet and Exponential Filtering r(n, m) is the measured paper property at CD bin n in the mth scan, Pmd is the MD estimate which is the same for all the CD bins in one scan, Pcd is the CD estimate, p is the weighting factor which is usually set as p 0.3. This will filter most of MD variations but will allow the CD estimate to change very slowly. In [15] a 2D wavelet denoising filter was devised to estimate CD and MD variations. To do this, a window of 32 scans was fed to a wavelet denoising filter. Furthermore, to emphasize the most recent scan, the last scan was repeated 3 times at the end of the window. The filter used coiflet wavelets with length 2 and hard thresholding with MultiMad noise estimation. In ad dition , for an array of 73 actuators and 685 sensors, the denoised third level approximation were used as basis weight CD estimates of the controllable profile. The MD profile is then computed as the mean value of each scan in the filtered data and the CD profile is found by subtracting the means from each scan. In this work, a similar filter is used for wavelet filtering but without stacking the last scan since doing this would produce artifacts in the estimated CD profile. An instance of this will be shown in the industrial data in section 4.8. In the following sections, the test scenarios developed in [10] are used as benchmarks to compare the performance of different algorithms under different scenarios. To measure the performance of each algorithm, the fol lowing error formulae are utilized. The error at CD position n and scan number m is defined as: Errd(n,m)=Ip(n,m)—.d(n,m)l, n=1,...,N; mr=1,...,M Thus, the error for each scan is defined by: Errd(m) Err,j(n, m) and the total error for m scans is J = Err(n, m) m n In all of the simulation sections 58 scans and 120 CD bins are assumed. 17 4.2. Edge Effects and Step Changes 4.2 Edge Effects and Step Changes In this scenario the effects of step changes in the CD profile and the effects ofprofile changes on the paper edges is investigated. The test profile is defined as: CD profile: pcd(n,m) = MD variations: —1 for1m29, 1n60 1 forl<m<29, 61<n<120 — for3om58, 1n60 for3om58,61<n<120 pmd(m) = sin(O(m)), 8(m) m(2) m l,2,...,60. The noise is white gaussian noise with zero mean and (cr)2 = (0.2)2. Figures 4.1 and 4.2 show the raw scans without and with noise. Ikrget Profile Figure 4.1: Edge Effects: raw scan without noise Figures 4.3 and 4.4 illustrate the estimated profiles produced by differentfilters. It can be seen that the estimated profile by DCT and Wavelets are closer to the target profile in comparison with the exponential filter. In these 1.5 1 2,, -0.5 18 4.2. Edge Effects and Step Changes Raw Piofile I2.5 2 L O5 I-1 I—2 —2.5 Figure 4.2: Edge Effects: raw scan with noise Figures the window length for DCT is 4. It takes longer for the estimated profile to reach a steady state when a DCT algorithm with a window of 32 scans is used. This effect is shown in Figure 4.5. Since it has been assumed that CD profiles are time invariant or slowly time varying, relaxing this assumption in this scenario examines an extreme case. Also note that all the methods suffer from the abrupt change in the CD profile. The errors obtained are shown in Table 4.1. Since the Haar wavelet has the least error in this case [10], the error for it is also included in this table. This low error for Haar wavelets is because of similarity between the square shape of the Haar mother wavelets and the step change. The Haar wavelet used here uses MultiMad, soft thresholding of length 2 and resolution 5 for denoising. 0 19 4.2. Edge Effects and Step Changes Method Jerror DCT4 6.85 DCT32 17.56 Exp 9.60 Wav0flet 11.12 WaVHaar 0.85 Table 4.1: The J error for different methods in the edge effects scenario 20 Ili rg et CD Pr of ile CD P ro fi j 0 E s tj it (W ay ) 0 — 1 60 Sc aji 5 CD Pr of ile E st na te 0 20 40 60 80 io o 12 0 Se ns or s — 1 60 S ca j 0 20 40 60 80 10 0 12 0 S en so rs CD Pr of ile E st jn t 0 (E xp ) t\Z Se rs or s I 0J F ig u 4. 3: Ed ge Ef fe ct s. e st im at ed pr of ile s Se ns or s Fi gu re 4. 4: Ed ge Ef fe ct s: e st im at ed pr of ile s — Ta rg et GD Pr of ile CD Pr of ile Es tim at e (W ay ) 4.3. Smooth Surfaces Figure 4.5: Edge Effects: estimated error versus time (scan number) The estimated MD variations are showii in 4.6. All the three MD esti mates are almost equal since they all use the scan average or average of the difference between the raw scan and the CD estimate as the MD estimate. 4.3 Smooth Surfaces In this scenario, sinusoidal variations are assumed as CD and MD variations to simulate the effect of smooth surfaces. The CD and MD are Pcd(fl) = sin(8(n)), 8(n) = n(6) n = 1,2, ..., 120. Pmd(m) = sin(O(m)), 8(m) m(2) rn = 1,2, ..., 60. with the noise being white gaussian noise with zero mean and (u)2 = (0.5)2. Figure 4.7 shows the CD plus the MD variations without any noise. Figure 4.8 shows the noisy raw scan used for estimation. Figures 4.9 and 4.10 show the target CD profile and the estimated profiles using different methods. In both figures it can be seen that DCT can estimate the target profile more accurately. This is also reflected in Table 4.2. This table also shows the error achieved when the Symmiet wavelets, the best wavelet filter for this Compaeison of Mean Scan Ereors Sean 23 4.3. Smooth Surfaces Figure 4.6: Edge Effects: estimated MD profiles Method Jerror DOT32 5.07 Exp 16.03 Wavcojflet 6.92 Wavsymmjet 6.68 Table 4.2: The J error for different methods in the smooth surfaces scenario case as shown in [10], is used with soft, MultiMad thresholding with length 6 to the third level for denoising. Figure 4.11 shows the estimation error for each scan. It can be seen that the DOT maintains a better performance compared to wavelet and exponential filtering. MD Profile Estimates 5 10 15 20 25 30 35 40 45 50 55 Scan 24 4.3. Smooth Surfaces Lhrget Profile Figure 4.7: Smooth Surfaces: raw scan without noise Figure 4.8: Smooth Surfaces: raw scan with noise Raw Profile 25 ‘ ili rg et CD Pr of ile CD Pr of ile Es tim at e (D CT ) I 1- 0 — 1 60 Sc an s 0 10 0 12 0 20 40 60 80 Se ns or s CD Pr of ile Es tim at e (W ay ) Sc an s 40 60 80 0 20 Se ns or s 12 0 CD Pr of ile Es tim at e (E xp ) 60 12 0 Se ns or s Sc an s 0 60 80 10 0 20 40 Se ns or s Fi gu re 4. 9: Sm oo th Su rfa ce s: e st im at ed pr of ile s Ta rg et CD Pr of ile CD Pr of ile Es tim at e (D CT ) 0. 5 — 0.5 — 1 80 10 0 12 0 I I5 0 40 30 (iD 20 10 I CD Pr of ile Es tim at e (W ay ) CD Pr of ile Es tim at e (E xp ) 20 40 60 Se ns or s Fi gu re 4. 10 : Sm oo th Su rfa ce s: e st im at ed pr of ile s 04.4. Grade Profile Change Comparison of Mean Scan Errors Figure 4.11: Smooth Surfaces: estimated error versus time (scan number) As shown in 4.12, estimated MD profiles are very similar. 4.4 Grade Profile Change In the paper production process, different grades of paper are producedbased on the customer demand. During a change in the grade being produced, the mean value of the raw profile changes rapidly. This scenario is simulated in this section. The CD and MD variations are assumed to be Pcd(Th) = 0.3cos(8(m)), n(2ir)8(n) = 120 n=1,2,...,120. m(2ir) 60 After scan 29, the MD amplitude is increased by 2. Figure 4.13 illustratesthis scenario. The noise added in Figure 4.14 is a white gaussian noise with zero mean and a standard deviation of 0.2. Figures 4.15 and 4.16 show the target and estimated profiles. It is clearthat the DCT can estimate the CD profile with less error. The errors for and Pmd(m) = 0.3sin(8(m)), m=1,2,...,60. 28 4.4. Grade Profile Change MD Profile Estimates 0.5 — .1ID Targct I — tVav —DCT E.sp 5 10 15 20 25 30 35 40 45 50 55 Scan llsrget Profile Figure 4.12: Smooth Surfaces: estimated MD profiles Figure 4.13: Grade Profile Change: raw scan without noise 29 4.4. Grade Profile Change Figure 4.14: Grade Profile Change: raw scan with noise each can be compared in Table 4.3. The symmiet here uses the same methodfor denoising as in the previous section but the length of wavelets is 8 andapproximation level 5 is used. This filter has achieved the best performancein this scenario according to [10]. IViethod Jerror DCT32 1.06 Exp 6.41 Wavc-,jflet 2.57 Wavsymmlet 1.02 Table 4.3: The J error for different methods in the grade change scenario Rw Profile 30 là rg et CD Pr of ile CD Pr of ile Es tim at e (D CT ) L2o • . : 4 0 6 : o i o o i 2o Sc an s 0 20 Se ns or s CD Pr of ile E st im at e (W ay ) Sc an s 0 1 0 — 1 60 10 0 12 0 60 80 20 40 Se ns or s CD Pr of ile Es tim at e (E xp ) c L Ct Se ns or s 10 0 12 0 20 20 40 60 80 Sc an s Se ns or s c A z I— s Fi gu re 4. 15 : G ra de Pr of ile Ch an ge : e st im at ed pr of ile s Ta rg et CD Pr of ile CD Pr of ile Es tim at e (D G) Fi gu re 4. 16 : G ra de Pr of ile C ha ng e: e st im at ed pr of ile s I I50 40 30 C!) 20 10 0 0 0. 5 0 - . 0. 5 — 1 50 40 30 C! ) 20 10 50 10 0 Se ns or s CD Pr of ile Es tim at e (W ay ) 50 40 10 50 40 30 C! ) 20 10 . ( 20 40 60 80 10 0 12 St- H SO LS CD PL ol ile Es tim at e (E xp ) . . 0 20 I 10 0 12 0 20 40 60 Se ns or s 80 10 0 12 0 L” z 4.5. Data Created by Moisture Variation Model Figure 4.17 shows the error for each scan under different methods. The wavelet error corresponds to the coifiet filter. - - Comparison of Mean Scan Errors 0.14 0.12 0.1 0 0 0.08 Q 0.06 ‘0 10 20 30 Scan 40 50 60 Figure 4.17: ber) Grade Profile Change: estimated error versus time (scan num The estimated MD variations are shown in 4.18. The estimated MD variations are close to the target MD variation. 4.5 Data Created by Moisture Variation Model A nonlinear moisture model was developed in [7] to capture the nonlinearity of the drying process. It was shown that there is a coupling between CD and MD variations. In addition, the MD variation is assumed to be Prnd(m) pmd(m) + m) where Pmd(m) is the mean moisture content and (m) is a zero mean stochastic process modeled with a first order dynamics: (m+ 1) = a(m) + w(m) where a is a known constant and w is a zero mean Gaussian white noise. Based on this model, the following CD and MD variations are used to establish a scenario: — Exp IVar - DCT 0.04 . 0.02 . . . . . . .. ... .... —. .-.—. —.— . —.—.—.— —. — 33 4.5. Data Created by Moisture Variation Model MD Profile Estimates iDTart 0.8 0.6 :; 35 40 45 50 55 Scan Figure 4.18: Grade Profile Change: estimated MD profiles Pcd(n, m) = sin(O(n)) * exp(—0.03 * m), (n) = n(6) n = 1,2, ..., 120. Pmd 0.5 + 1 — 0.8z’ w, w N(0, (0.1)2) In Figures 4.19 and 4.20 the raw profile with and without noise is shown.The performance of different filters is shown in Figures 4.21 and 4.22. Sincethe CD profile fades quickly with time, a shorter window length is picked forthe DCT method. It can be seen that both wavelet and DCT, with windowlength 5, can estimate the profile better than the exponential filtering. Table4.4 shows the errors in each case. The error achieved by using a DCT withwindow length 32 is included in this table. The coiflet filter has the leasterror among the wavelet filters [10], the DCT filters and the exponentialfilter. 34 4.5. Data Created by Moisture Variation Model Method Jerror DCT5 9.25 DCT32 14.01 Exp 11.43 Wavcojflet 8.34 Table 4.4: The J error for different methods in the moisture variation model scenario Target Poffle Figure 4.19: Moisture Variation model: raw scan without noise 1.5 40 Scans 20 0 40 Sensors 35 4.5. Data Created by Moisture Variation Model Raw Profile 2 1.5 0.5 —0.5 Figure 4.20: Moisture Variation model: raw scan with noise 0 36 c t I Ta rg et CD Pr of ile 60 CD Pr of ile Es tim at e (D CT ) Sc an s 1 0 — 1 60 Se ns or s CD Pr of ile Es tim at e (W ay ) Sc an s 40 60 80 10 0 0 20 Se ns or s CD Pr of ile Es tim at e (E xp ) 60 Se ns or s 20 Sc an s 0 60 80 20 40 Se ns or s 12 0 Fi gu re 4. 21 : M oi st ur e V ar ia tio n m o de l: e st im at ed pr of ile s Ta rg et CD Pr of ile CD Pr of ile Es tim at e (D CT ) . I: ij 3 0 C, ) 20 10 50 40 30 C ,) 20 10 I 50 S en so rs CD Pr of ile Es tim at e (W ay ) 50 40 30 so 20 10 I. I 60 Se ns on 80 10 0 12 0 Fi gu re 4. 22 : M oi st ur e V ar ia tio n m o de l: e st im at ed pr of ile s 4.6. Bumps and Streaks The error for each scan is shown in Figure 4.23. This Figure confirms thatthe coifiet filter has the best performance in this scenario. In addition, theDCT method with window length 5 can follow the fast changing CD profilemuch better in comparison with the DCT method with window length 32. Comparison of Mean Scan Errors -S , —Exp Was, 015 005 S ,‘. — #1 •. -.5 0 0 10 20 30 40 50 60 Scan Figure 4.23: Moisture Variation model: estimated error versus time (scannumber) 4.6 Bumps and Streaks This scenario simulates very quickly changing local variations in the CDprofile. The CD and MD are defined as: Pcd(fl, m) = A(n, m) cos(O(n)), 8(n) = n(2) n = 1,2,..., 120. where: —1 for30m39, 30n34 1 for2o<m<23, 80<n<85A(n,m)= 0.6 for 1<m<58, 60<n<61 0.2 elsewhere 39 4.6. Bumps and Streaks and m(2ir)Pmd(m) = 0.3sin(&(m)), 0(m) = 60 m= 1, 2, ..., 60. and (u)2 = (0.2)2 for the added white noise. Figures 4.24 and 4.25 show theeffect of noise on the target profile. Target Profile Figure 4.24: Bumps and Streaks: raw scan without noise The estimated profiles are illustrated in Figures 4.26 and 4.27. Theexponential filter has detected both the bumps and the streak but noise isalso present in the estimated profile. The wavelet filter has reduced theeffect of the noise but can not detect the streak . As shown in Figure 4.27,the profile estimated by DCT with window length 32 scans has achieved analmost noiseless CD estimate and has detected the streak. However, only avague image of the big bump is present in the profile and the smaller bumpis completely vanished. This can be explained by the fact that the DCTmethod is vulnerable to fast varying CD profiles like bumps.As shown in Figure 4.28 the error achieved by the DCT method which useswindow length three is the lowest. The optimum window length which givesthe lowest error is achieved by finding the total error for different windowlengths from 2 to 32 and picking the window length which has the minimumerror. The total error is lowest for the DCT method with window length of3 scans as shown in Table 4.5. 0 40 4.6. Bumps and Streaks Raw Profile I I 0.5 0 Figure 4.25: Bumps and Streaks: raw scan with noise Method Jerror DCT3 7.02 DCT32 8.59 Exp 8.45 Wavcojfjet 8.94 Table 4.5: The J error for different methods in the bumps and streaks scenario 1: —0.5 —1 41 CD Pr of ile E st im a te (D CT ) 12 0 CD Pr of ile E st im at e (E xp ) CJD I 12 0 T tr gc t CD Pr of ile 60 80 10 0 12 0 Sc an s 0 20 40 S en so rs CD Pr of ile E st im at e (W ay ) - • . • — V Sc an s 0 20 40 60 80 Se ns or s Se ns or s Sc an s 0 20 40 60 80 Se ns or s Fi gu re 4. 26 : B um ps an d St re ak s: e st im at ed pr of ile s CD Pr of ile Es tim at e (D cr ) I C,) CD Pr of ile Es tim at e (W ay ) a Cl ) a 0 U ) C, ) Se os or s CD Pr of ile Es tim at e (E xp ) 60 Se os or s 10 0 12 0 Fi gu re 4. 27 : B um ps ar id St re ak s: e st im at ed pr of ile s 4.7. Spanning the MD Frequencies 01 C 8 Scan 60 Figure 4.28: Bumps and Streaks: estimated error versus time (scan number) The estimated MD variations are shown in 4.29. Once again, the MD estimates are almost the same for the three methods. 4.7 Spanning the MD Frequencies This scenario is similar to the smooth surfaces scenario but the objective here is to investigate the effect of different MD frequencies on the CD/MD separation methods. The CD profile is shown in Figure 4.30. Figures 4.31 and 4.32 depict the raw profile with and without noise respectively. The MD variation is a cosine function. The frequency of the MD variation is swept from 0.005 Hz to 0.5 Hz and the estimation error at each frequency is shown in Figure 4.33. It can be seen that the DCT method achieves lower errors over the whole range compared to the wavelet and exponential methods. The peak errors in this graph occur at multiples of - where T5 20 seconds. This phenomenon happens since at these frequencies the MD and CD variations are inseparable [16]. — 0.04 0.0k 0 10 20 30 40 50 44 4.7. Spanning the MD Frequencies MD Proflic Estimates Figure 4.29: Bumps and Streaks: estimated MD profiles Figure 4.30: The target CD profile 2 1.5 0.5 0 —0.5 —1 5 10 15 20 25 30 35 40 45 50 55 Scan CD Target Profile —2. 100 0 45 4.8. Comparison of the Estimated Profiles Using Industrial Data 2.5 2 1.5 0.5 0 —0.5 —1.5 —2 Figure 4.31: The raw scan without noise 4.8 Comparison of the Estimated Profiles Using Industrial Data In this section, the DCT method is applied to process data from the same paper machine as in [15]. There are 600 data bins across the sheet and 300 scans are recorded. The original data has NAN (not a number)’s on both sides of each scan corresponding to when scanner goes off the sheet. These NAN’s are later replaced with zeros and then trimmed before they are fed in to the DCT algorithm. The DCT algorithm uses windows of length 32 and the 1-D wavelet filter discussed in 3.3. The data for the controllable CD profile estimated by the wavelet filter is saved from the implemented online CD/MD separation block in [15]. Figure 4.34 shows the raw scan and Figures 4.35, 4.36 and 4.37 show the estimated profile using Exponential filtering, wavelet method and DCT method respectively. It can seen that the profile estimated by wavelet has many fast time varying CD variations which can not correspond to the slow time varying CD profile. This is due to the limited window size (32 scans) and 3 times repetition of the last scan (section 4.1) in the online 2v!D + CD Target Profile A 0 46 4.8. Comparison of the Estimated Profiles Using Industrial Data 90 80 70 60 50 40 30 20 10 Figure 4.32: The raw scan with noise The CD estimation Error for different MD frequencies 4 3 2 0 —1 —2 —3 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Fir Figure 4.33: The CD estimation Error for different MD frequencies 0 Way 47 4.9. Conclusion implementation. The estimated profile by the exponential filter has less fastchanging CD variations. Figure 4.37 shows that the DCT estimated profileis very smooth and has much less energy in the high frequency components. 4.9 Conclusion Figure 4.34: Raw Scans In this chapter, the errors achieved by the traditional CD/MD separationmethods are compared with the errors from the DCT based technique. Tosummarize, when no fast transient variations exist in the CD profile, thesuggested method can achieve lower errors in estimating the CD profile thanboth the wavelet and exponential filters. The MD profiles are very similarfor the low MD frequencies. For higher MD frequencies, the mean of each scan can not be consideredas the true variation. The DCT method has an estimate of the MD at eachdata point. If this estimation is properly filtered, the estimated MD can beused to control MD variations more effectively. The application of the DCT method to real measurements shows visuallythat the estimated CD profile is less time varying and smoother compared tothe wavelet filter currently in use. However, the real online test, which canallow a fair comparison between the two methods, will be when the profile :1 48 4.9. Conclusion CD Profile Estiinatc (Exp) 2 300 , -. i ::: _____ _____ 0 150 . ___ 100 I —1 l : 100 200 300 400 500Sensors Figure 4.35: Estimated profile using exponential filtering CD Profile Estlnoite (\Vnv) 2 30C -.. 15 _ 1 250 -1. ‘! -:j - 100 200 300 400 500 Sensors Figure 4.36: Estimated profile using wavelet 49 4.9. Conclusion 2 1.5 0.5 S 0 -0. Figure 4.37: Estimated profile using DCT estimated by the DCT method is used inside the CD control loop. 50 Chapter 5 A Picketing Detection Algorithm In this chapter a novel picketing detection method is suggested. The idea is explained in 5.1. The detection algorithm used in the method is discussed in 5.2. The simulated results are presented in 5.3. 5.1 The Idea Picketing is an undesirable unstable pattern of variations seen in the paper profile and is often considered a form of instability. This pattern is charac terized by a profile which shows adjacent undesired large amplitude streaks resembling a “picket fence”. An instance of this is shown in Figure 5.1. This occurs because the corresponding adjacent actuators are either fully opened or fully closed. A similar pattern is present in the actuator profile corresponding to this raw profile seen in Figure 5.2. Two main reasons that can cause this type of instability are aggressive tuning of the controller and misalignment. In the first case, the control signal is applied in frequencies where the process gain is low and model uncertainty is high [27]. Misalignment can also contribute to the high un certainty at the high frequencies. In the second case, the operator has to find the correct alignment after the picketing has been visually detected. This is not desirable since the corrective action can only be taken after the problem has been detected by the operator. The goal in this chapter is to devise an online picketing detection scheme that can either alarm the operator or send a signal to the controller about an impending picketing. An intuitive way to detect this pattern in the actuator profile is to ob serve the high frequency components in the actuator profile. An abnormal increase in the power of the high frequency components can signify a growing picketing pattern on the raw profile. The method devised here uses FFT to find the spectrum of each actuator profile. From this spectrum, the high frequency components are extracted 51 Figure 5.1: A raw profile which has picketing 100 90 80 70 60 50 40 30 20 10 0 5.1. The Idea 600 500 400 3 2 0 —1 hoo C’) 200 c0 Cr2 Figure 5.2: The actuator profile corresponding to the raw scan 52 5.2. CUSUM Algorithm and recorded over time. Later, a change detection algorithm (discussed in 5.2) detects the increase in these components. There are different parts of spectrum that can be extracted as high fre quencies. Assume that the power of each scan at frequency f is denoted byPf. Also assume that actuators are equally spaced at a Xa distance from each other, the highest frequency for the spectrum is the Nyquist frequency The power of this highest frequency Pfi can be monitored for 2Xachanges and so used to detect picketing. A drawback of this method is the selectivity used in this approach. On actuators with wide responses, the in crease in power will not occur at exactly the highest frequency but rather a range of high frequencies. Thus, in the second approach, the relative power of frequencies beyond the controller bandwidth is monitored for changes. In other words P defined as equation 5.1 will be monitored. Pf 2X (5.1) z; Pf In Equation 5.1, X denotes the controller bandwidth which is typically between 2Xa and 2.SXa [28]. 5.2 CUSUM Algorithm To detect the change in Pj and Pr, they should first be normalized. To2Xa normalize these two variables, they are divided by their mean value when there is no picketing present in the profile. These mean values have to be found heuristically by recording the values P and P. over time. 2XaOnce P or Pr is normalized, it can be fed to the CUSUM algorithm2Xa which can detect an increase of the mean values online( [29]). The algorithm is as follows: S = Si_i + antiw * ( — 1) (5.2)mean(P) m = rninS (5.3)J<i (5.4) The normalized P is fed to the algorithm as the input. If the mean of P increases, then mean(P)) — 1 increases and the result will accumulate in S. 53 5.3. Simulation Results The increase in S is detected by comparing its values with its minimum so far. By thresholding the output g, which is equal to the difference of S and the minimum, a flag can be set to alarm the increase in mean of P. The antiw coefficient works similarly to the anti windup in PID controllers. The default value for antiw is 1 and it is only set to zero when the threshold is passed. This allows the output g to decrease faster after the picketing pattern fades in the profile. This method will be used for both P.j and P separately in two simu 2Xalation scenarios in the next section. 5.3 Simulation Results 5.3.1 Misalignment Scenario In this scenario, a fake misalignment is introduced to the closed-loop system in order to cause picketing in the profile. Figure 5.3 shows the raw pro file. The picketing can be visually detected starting around scan 500. The corresponding actuator profile is illustrated in 5.4. The picketing can also be noticed here from scan 500. In Figure 5.5, the normalized power of the highest frequency is displayed. The system reaches the steady state at scan 50 and the misalignment is introduced at scan 90. Thus the average from scan 50 to 90 is used to normalize the power. It can be seen that the power increases monotonically from scan 200. The average shown in this Figure is equal to 1 since the power has been normalized. In Figure 5.6, the output of the CUSUM algorithm is displayed. The threshold is set to 5. The graph reaches 5 around scan 90 and then decreases until it increase rapidly after scan 240 to pass the threshold around scan 260. This enables the operator to run an alignment test before a low quality paper is erroneously produced. Figure 5.7 shows the P over time. Here also the power rises sharply after around scan 200. In Figure 5.8, the graph for g values passes the threshold of 5 around scan 90. In both cases, the method can detect the growing pattern of picketing almost 300 scans or 100 minutes (T3 is assumed 20 sec) before any symptom can be seen on the paper sheet. 5.3.2 No Spatial Filtering Scenario In this case, the spatial filtering on the profile is removed. This will allow the controller react to high frequency components and eventually result in picketing. Figure 5.9 shows the raw measurements and Figure 5.10 displays the actuator profiles. As shown in Figure 5.10 the spatial filtering is removed 54 5.3. Simulation Results 94 93.5 93 92.5 92 91.5 91 90.5 90 89.5 89 15 10 5 Figure 5.3: The raw profile 900 Actuators Profile for Basis Weight C) Ci) Figure 5.4: The actuator profile 55 — o for P1500 / TIuc,IioId 1000 100 200 300 400 500 600 700 800 900 Scan Numr Ji 100 150 200 Scan Number Figure 5.6: The g values for P1 versus scan number 2Xa 5.3. Simulation Results 100 200 300 400 500 600 700 800 900 Scan Number Figure 5.5: The power of highest frequency component versus scan number 50 250 300 56 1500 1000 500 5.3. Simulation Results Figure 5.7: P versus scan number 100 200 300 400 500 600 700 800 900 Scan Number Ut ________ g for P,. 40 Threshokl 20 50 100 150 Scan Number 200 250 300 Figure 5.8: The g values for Pr versus scan number 57 Normalized Power of the Highest Ebequency 100 200 300 400 500 600 700 800 900 Scan Number 5.4. Conclusion at scan number 85. The P.j_ and its g values are shown in 5.11 and 5.12. The sharp rise in P shows that the high frequency component in the actuator profile is growing. With the same thresholds as in previous section, the CUSUM detects picketing at scan 89. This is 200 scans (65 minutes for a scanner with scarming time 20 seconds) before the picketing can be seen on the profile. Figure 5.9: The raw profile Figure 5.13 shows the growth of P versus time. It can be seen that P increases much faster after the spatial filtering is removed. This fact is also reflected in the g values in figure 5.14. Here the CUSUM detects the picketing at scan 88. 5.4 Conclusion Two approaches were investigated for an online picketing detection scheme. It can be seen that both approaches can detect picketing well in advance with an acceptable accuracy. The second approach uses relative power to signal picketing which can be more useful for wide response actuators. The only drawback for the proposed change detection method is that the CUSUM algorithm needs to be manually fed with the right expected value of powers in no picketing steady states. Raw Profile for Basis Weight 94 93.5 93 92.5 92 91.5 91 90.5 90 89.5 89 a Cl) 50 101) 150 200 250 Sensors 58 180 160 140 120 100 80 60 40 20 5.4. Conclusion I1510 5 0 —1 —1 Actuators Profile for Basis Weight 4 : 15 0 !50 d5 Actuators Figure 5.10: The actuator profile —P — Average of P 50 100 150 200 250 300 350 Scars Number Figure 5.11: The power of highest frequency component versus scan number 59 5.4. Conclusion Figure 5.12: The g values for P versus scan number 2Xa 7 6 5 4 3 2 Normahzed Power of the Htghcst Frequency 50 100 150 200 250 300 350 Scan Number Figure 5.13: P. versus scan number 50 100 150 200 250 300 350 Scan Number 150 Scan Number 300 — Aerage of P, 60 5.4. Conclusion 2001. _________________ —gforP 1500 Threshold 50 100 150 200 250 300 350 Scan Number o:P 50 100 150 200 250 300 Scan Number Figure 5.14: The g values for PT versus scan number 61 Chapter 6 Conclusions and Future Work 6.1 Conclusions In this thesis, a novel online CD/MD separation method is presented. It is shown that the CD variations form an even periodic signal with period 2T8 inside the scanned signal. This property is used by Discrete Cosine Transform to pick the frequency components which correspond to the CD profile. Next a 1-D wavelet ifiter, separates the controllable variations from uncontrollable variations. The MD variation is computed as the mean of the difference between the the raw profile and the estimated CD profile. The performance of this method is evaluated through different simulation scenarios. The results are compared with exponential filtering and 2-D wavelets which are conventional in industry. It is shown that DCT can achieve lower estimations errors and smoother profiles. The major difficulty is when there is a fast changing CD profile or there are MD variations with frequencies around twice the scarming rate. In the first case, the suggested method suffers because it is assumed in the model that the CD variations are time invariant or slow time varying. Also there is a tradeoff between the accuracy and tracking of any online algorithm. In the second case, it has been shown that no constant scan rate method can differentiate the CD pro file from MD variations( [16]). The simulation results are verified in section 4.8 by applying the method to real measurements from a paper machine. The results from online implementation of this method will be published in near future. In the second part of this thesis, an online picketing detection method is presented. A form of instability known by the industry as “picketing” is associated with growing high frequency components in the actuator profile. This property is used in conjunction with a change detection algorithm, 62 6.2. Future Work CUSUM, to detect the increase in the high frequency components. Two variations of this method are evaluated by different scenarios which simu late picketing. It is shown that both methods can detect and predict the picketing long before any visible effect of picketing on the paper profile. The drawback of the CUSUM method is that it needs the average value of the signal in steady state. This method also has been implemented on a real paper machine and the results will be published in the near future. 6.2 Future Work The DCT method discussed in 3.3 has been implemented on a paper ma chine. The final trials to evaluate the performance of the method in the loop plots are planned. The results from this implementation can be used for better tuning of the window size or the wavelet filter. In addition, MD variations were estimated by averaging the difference between the the raw profile and the estimated CD profile. By using a better filter to remove the measurement noise, a better estimation for MD variations can be achieved which can increase the MD control bandwidth significantly. The internal even symmetry discussed in section 2.2 can be exploited to reduce the effect of noise and MD variations on the CD profile estimate. The even periodic property of CD profiles can be exploited in a CD/MD separation scheme for a scanner with two scanning speeds presented in [18]. The variable speed scanner has the potential to gain more information in comparison to fixed speed scanners. Whether the combined solution can im prove the performance is the subject of further research. Another potentially useful method which can be applied to the CD/MD separation problem is the compressive sensing ([30]). The picketing method discussed in Chapter 5 can be used in series with an online misalignment detection scheme as in [31] and [32], to detect a misalignment and find the right alignment automatically. The correct align ment can then be forwarded to the CD controller. Currently this process is carried out manually in paper mills. Also the CUSUM algorithm discussed in 5.2 can be replaced with more complex algorithms to reduce the tuning costs. The upgrade would be at cost of more computational cost. 63 Bibliography [1] J. Ghofraniha. Cross directional response modeling, identification and control of dry weight profile on paper machines. PhD thesis, University of British Columbia:191—211, July 1997. [2] G.A. Smook. Handbook for Pulp and Paper Technologists. Angus Wilde Publications Inc., 1992. [3] E.W. Leaver. An Idea is Not Enough: triumphs and disasters during 50 years as a high tech enrepreneur. Peregrine Press, London;Ontario, first edition, 1998. [4] J.G. VanAntwerp, A.P. Featherstone, R.D. Braatz, and B.A. Ogun naike. Cross-directional control of sheet and film processes. Automatica, 43:191—211, 2007. [5] E. Dahlin. Computational methods in a dedicated computer system for measurement and control on paper machines. TAPPI, pages 62.1—62.42, Sep 1969. [6] S.C. Chen. Kalman filtering applied to sheet measurement. Proceedings of the American control conference, pages 643—647, 1988. [7] C. Lindeborg. A process model of moisture variations. Pulp and Paper Canada, pages T 142—T 147, April 1986. [8] G.A. Dumont, M.S. Davies, K.Natarajan, C. Lindeborg, F. Ordubadi, Y. Fu, K. Kristinsson, and I. Jonsson. An improved algorithm for estimating paper machine moisture profiles using scanned data. IEEE Transactions on Control Systems Technology, 1(2):101, 113 1993. [9] X.G. Wang, G.A. Dumont, and M.S. Davies. Estimation in paper ma chine control. IEEE Control Systems Magazine, 13(4) :34—43, August 1993. [10] Z. Nesic. Paper machine data analysis using wavelets. Master’s thesis, University of British Columbia, February 1996. 64 Chapter 6. Bibliography [11] Z. Nesic, M.S. Davies, and G.A. Dumont. Paper machine data com pression using wavelets. Proceedings of the 1996 IEEE International Conference on Control Applications, pages 161—166, September 1996. [12] K.K. Kan. Real time paper machine data wavelet analysis and cross directional actuator response identification. Master’s thesis, University of British Columbia, April 2004. [13] X. Jiao. Paper machine data analysis and optimization using wavelets. Master’s thesis, University of British Columbia, January 1999. [14] Z. Nesic, M.S. Davies, and G.A. Dumont. Paper machine data analysis and compression using wavelets. TAPPI, 80(10):191 — 204, October 1997. [15] G. Dumont, J. Ball, G. Emmond, and M. Guillemette. Paper machine cd control using wavelet-based md/cd separation. In Advanced Process Control Workshop, April 2003. [16] S. Aslani. Estimation of cross direction and machine direction variations using recursive waviet filtering. Master’s thesis, University of British Columbia, February 2004. [17] S. Aslani, M.S. Davies, G.A. Dumont, and G.E. Stewart. Estimation of cross and machine direction variations using recursive wavelet filtering. IEEE Advanced Process Control Workshop, April 2004. [18] S. Aslani, M.S. Davies, G.A. Dumont, and G.E. Stewart. Detecting aliasing between cross and machine directions variations by variable sampling rate. Control Systems, 2004. [19] J.C. Skelton, P.E. Welistead, and S.R. Duncan. Distortion of web pro files by scanned measurements. Pulp and Paper Canada, 104:12:81—84, 2003. [20] S. Duncan and P. Wellstead. Processing data from scanning gauges on industrial web processes. Automatica, 40:431—437, March 2004. [21] 5. Duncan and A. Taylor. Using signals from two gauges to estimate md variations for cd control systems. Proceedings of the 2007 American Control Conference, pages 1341—1346, 2007. [22] A.V. Oppenheim, R.W. Schafer, and J.R. Buck. Discrete-Time Signal Processing. Prentice-Hall mc, Upper Saddle River, New Jersey 07458, second edition, 1999. 65 Chapter 6. Bibliography [231 Sylvain Gendron. Fundamental properties of scanner signals for thecorrect separation of paper machine cd and md variations. Pulp andPaper Technical Association of Canada 91st Annual Meeting Preprints,pages A289—A293, 2005. [24] N. Ahmed, T. Natarjan, and K.R. Rao. Discrete cosine transform.IEEE Transactions on Computer, pages 90—93, January 1974. [25] K.R. Rao and P. Yip. Discrete Cosine Transform: Algorithms,Advantages, Applications. Academic Press, 1990. [26] G Strang. The discrete cosine transform. SIAM Review, 41(1):135—147,March 1999. [27] G.E. Stewart. Two dimensional loop shaping controller design for papermachine cross-directional processes. PhD thesis, University of BritishColumbia, August 2000. [28] G.E. Stewart, J.U. Backstrom, P. Baker, C. Gheorghe, and R.N. Vyse.Controllability in cross-directional processes practical rules for analysisand design. Pulp and Paper Canada, 103(8):T 208 — T 214, 2002. [29] lvi. Basseville and I.V. Nikiforov. Detection of Abrupt Changes: Theoryand Application. Prentice-Hall, mc, April 1993. [30] Emmanuel Candes. Compressive sampling. Proceedings of theInternational Congress of Mathematicians, 2006. [31] F. Farahmand, G. Dumont, P. Loewen, and M. Davies. Tensor-basedblind alignment of mimo cd processes. Journal of Process Control, 2008. [32] F. Farahmand. Alignment, modeling and iterative adaptive robustcontrol of cross directional processes. PhD thesis, University of BritishColumbia, 2009. 66 Chapter 7 Appendix A - Matlab Programs for Simulations in Chapter 4 7.1 The Main Program 0 /0 Y0 MainSim.m ‘I !0 °h This is the main program for the simulations in ‘I, my thesis. °h (c) Soroush Karimzadeh File created: Nov 12,2007 Modified : Oct 20,2008 0/ /0 0/ /0 clear all close all dc °hinitializaing °h Choose : 1 for Step Change 2 for Smooth Surfaces 3 for Grade change 4forMVM 5 for Bumps 7, Produce the raw data 7,Graph raw profile surf (CDsignal+MDsignal) 67 7.1. The Main Program view(—38,70) xlabel(’Sensors’ , ‘interprete’,’latex’) ylabel(’Scs’,’interpreter’, ‘latex’) title(’MD +CD Target Profile’ ,‘interpreter’, ‘latex’) axis(t1 120 0 60 —2 2]) colorbar(’clim’,[—1.5 2]) ‘hGraph target profile figure surf (dwt) view(—38,70) xlabel(’Sensors’ ,‘interpreter’, ‘latex’) ylabel(’Scans’ ,‘interpreter’, ‘latex’) title(’ Raw Profile’ ,‘interpreter’, ‘latex’) axis([1 120 0 60 —2 2]) colorbar(’clim’ ,[—1.5 2]) 0 I0 [D,M]=size(dwt);Yfinding the size of the data: D is the %number of scans M is the number of CD bins [Ecd,Emd,Emd2]=DCT(dwt) ; °hApply the DCT method [Ecd2 , Emmmd,Emmmd2] =DCTv5(dwt ,3) ; XApply the DCT method. °hApplying the Exponential filtering [dwtExp,Emdex,Ecdex]=expfiltl(dwt ,0 .8); °h The Wavelet dunction used here is developed by 7 Zoran Nesic ‘h Inputs: °h PPCmain(dwt, wavLength, wavLevels, wavTypeNum, °h NoiseMethod, ThreshTypeNum,Flag) °h dwt 1D or 2D input data set °h wavLength wavelet length °h wavLevels the number of decomposition levels 0h wavTypeNum 1 - “daubcqf”, 2 - “haar”, 7, 3 — “Symmiet”, 4 — “Coif let” 0/, NoiseMethod noise estimation method (1 — std, 2 - MAD, 3 - MultiMAD) 68 7.1. The Main Program hard thresholding% ThreshTypeNum 1. - “hard”, 2 - “soft” °h Flag 0 — return wavelet coefficients, 1—return °h details and approximations [dwtNewF,wcoeff 1 ,wcoeff ,wavLevels, wavThresh,THlevels]=...PPC_main(dwt, 2, 3, 4,3,1,0); ‘I. estimate the MD and CD varitions for wavelet[MD, MDexp, MD_way, CD, CD_exp, CD_way] = cdmdcalc(dwt ,dwtExp,dwtNewF); MD_wav= scanrev(MD.wav,M,D);Y. change 2D MD estimate to 1D0 !0 %Graphing the CD estimate for OCT. Wavelets7.and Exponential filters using surf figure(’Units’,’inches’,’Position’,[O 0 10 5]) set (gcf, ‘Paperarientation’ ,‘landscape’) subplot (2,2,1) surf (CDsignal) xlabel(’Sensors’ ,‘interpreter’, ‘latex’)ylabel(’Scans’,’interpreter’, ‘latex’) title(’Target CD Prof ile’,’interpreter’, ‘latex’) axis([1 120 0 60 —1.5 1.5]) subplot (2 ,2, 2) surf (Ecd) xlabel( ‘Sensors’ , ‘interpreter’, ‘latex’)ylabel(’Scans’ ,‘interpreter’, ‘latex’) title( ‘CD Profile Estimate (DCT) ‘ , ‘interpreter’ , ‘latex’)axis([1 120 0 60 —1.5 1.5]) subplot(2,2,3) surf (CD_way) xlabel( ‘Sensors’ , ‘interpreter’, ‘latex’)ylabel(’Scans’ ,‘interpreter’, ‘latex’) title(’CD Profile Estimate (Way) ‘,‘interpreter’, ‘latex’)axis({1 120 0 60 —1.5 1.5]) subplot(2,2,4) surf (CD_exp) xlabel( ‘Sensors’, ‘interpreter’, ‘latex’)ylabel(’Scans’ ,‘interpreter’, ‘latex’) title(’CD Profile Estimate (Exp)’ ,‘interpreter’, ‘latex’)axis([1 120 0 60 —1.5 1.5]) 0 JO 69 7.1. The Main Program °hGraphing the CD estimate for DCT, Wavelets Y,and Exponential filters using imagesc figure(’Units’ ,‘inches’, ‘Position’, [0 0 10 5]) set(gcf, ‘PaperOrientation’, ‘landscape’) subplot(2,2,1) imagesc([0,M] , [0 Dl ,CDsignal, [—1 1]) set (gca, ‘YDir’, ‘normal’) xlabel ( ‘Sensors’, ‘interpreter’, ‘latex’) ylabel(’Scans’, ‘intereter’,’latex’) title(’Target CD Profile’, ‘interpreter’, ‘latex’) colorbar( ‘WestOutside’) °haxis([1 120 0 60 —1.5 1.5]) subplot(2,2,2) imagesc(Ecd, [—1 1]) set (gca, ‘YDir’ ,‘normal’) xlabel (‘Sensors’ , ‘interpreter’, ‘latex’) ylabel(’Scans’,’interpreter’, ‘latex’) title(’CD Profile Estimate (DCT) ‘,‘interpreter’, ‘latex’) %axis([1 120 0 60 —1.5 1.5]) subplot (2 ,2, 3) imagesc(CD_wav, [—1 1]) set (gca, ‘YDir’ ,‘normal’) xlabel(’Sensors’ ,‘interpreter’, ‘latex’) ylabel C ‘Scans’, ‘interpreter’, ‘latex’) title(’CD Profile Estimate (Way) ‘,‘interpreter’, ‘latex’) 7,axis([1 120 0 60 —1.5 1.5]) subplot(2,2,4) imagesc(CD_exp, [—1 1]) set(gca, ‘YDir’ ,‘normal’) xlabel (‘Sensors’ , ‘interpreter’, ‘latex’) ylabel( ‘Scans’, ‘interpreter’, ‘latex’) title(’CD Profile Estimate (Exp)’, ‘interpreter’ ,‘latex’) 7,axis([1 120 0 60 —1.5 1.5]) •1 10 °hFinding the J error for different methods j...exp=sqrt(sum(sum((CD..exp(1:D,:)—CDsignal(1:D,:)).2))); j_wav=sqrt(sum(sum(( CD_wav(1 :D, :)—CDsignal(1:D, :)) . j...dct=sqrt(sum(sum(( Ecd(1:D,:)—CDsignal(1:D,:)).’2))); 70 7.1. The Main Program °hFinding the Error at each scan for different methodsfor i=1:D CDerexp(1,i)= sum(abs(CD..exp(i,:)—CDsignal(i,:)))/M;CDerwav(1,i)=’ sum(abs(CD_wav(i,:)—CDsignal(i,:)))/M;CDerdct(1,i)= sum(abs(Ecd(i,:)—CDsignal(i,:)))/M;CDerdct2(1,i)= sum(abs(Ecd2(i,:)—CDsignal(i,:)))/M;end Oj 10 °hplotting the Error at each scan for different methodsfigure plot(1:D,CDerexp,1:D,CDerwav, ‘r: ‘,l:D,CDerdct, ‘g-. ‘,l:D,...CDerdct2, ‘c——’, ‘LineWidth’ ,2); legend(’Exp’ , ‘$Wav$’, ‘$DCT_32$’, ‘$DCT_3$’,’interpreter’,... ‘latex’); set(legendO, ‘interpreter’, ‘latex’); ‘/. plot(1 :D,CDerexp, 1 :D, CDerwav, ‘r: ‘,1 :D, CDerdct, ‘g—.’,...°h’LineWidth’ ,2); % legend(’Exp’,’$Wav$’,’$DCT$’,’interpreter’,’latex’);set (legendO,’interpreter’, ‘latex’); xlabe1(’Sc’,’interpreter’, ‘latex’) ylabel(’CD error’,’interpreter’, ‘latex’) Y,Plotting the ND estimates figure xes=linspace(1 ,D,M*D); plot(xes,MDsignal, ‘b’ ,xes,MD_wav, ‘r’ ,xes,Emdex, ‘y’ ,xes,Emd2,... ‘g’ , ‘LineWidth’ ,2); legend(’$MDTarget$’, ‘$Wav$’, ‘$Exp$’, ‘$DCT$’, ‘interpreter’,... ‘latex’); set (legendO, ‘interpreter’ ,‘latex’); xlabel(’Scan’ ,‘interpreter’, ‘latex’) °hylabel(’MD error’,’interpreter’, ‘latex’) title(’MD Profile Estimates’, ‘interpreter’ ,‘latex’)axis([1 D —2 2]) grid on 0 10 function [dwt ,n,m,CDsignal ,MDsignal] = MkTstSigv2(TestData) °h MkTstSigv2.m 71 7.1. The Main Program 0 10/, This function creates different test signals forY0 my thesis. It is called by FigXXXX.m files. °h (c) Soroush Karimzdeh File created: Dec 11, 1995 by Zoran Nesic ‘fl Last modification: Jun 28, 19960/ Modified for DCT : Oct 20,20080/ /0 0/ /0 if TestData == 1 0/ 10/ First data set used for thesis 0/ /0 T = 10; % time in seconds n = 58; X number of scans m = 120; 0fl number of sensors MDamp = 1; /, MI) amplitude w = 2*pi/T; °fl MD omega CDamp = 1; ‘/, CD amplitude tall = 0.0; °h CD time constant bumpAmp = 0.5; signal = 0.2; °h noise level MDfun_name = sprintf(’a(l) * sin( a(2)*TimeVector )‘);CD±un_name = sprintf(’a(l) * exp( a(2)*TimeMatrix )‘);Profile = [—ones(1,m/2) ones(1,m/2)]; t = linspace(O, T, n*m); ‘I, TimeVector[MDsignal, TimeMatrix] = MkMDsig(MDfun_name, t,...[MDamp WI, [n,m]); CDsignal = MkCDsig(CDfun_name,TimeMatrix,[CDamp tau] ,Prof lie); CDsignal(n/2+l :n,:) = bumpAmp*CDsignal(n/2+1 :n, :); randn(’seed’ ,3133); dwt = MDsignal + CDsignal + sigmal*randn(n,m); elseif TestData == 2 0/ /0/. Second data set used for thesis 0/ /0 72 7.1. The Main Program T 1000; °fl time in seconds n = 58; °h number of scans m = 120; °h number of sensors MDamp = 1; 7. MD amplitude w = 2*pi/T; °h MD omega CDamp = 1; °h CD amplitude tan = 0.0; % CD time constant bumpAmp = 0.5; sigmal = 0.5; Y noise level MDfun_name = sprintf(’a(l) * sin( a(2)*TinieVector )‘); CDfun_name = sprintf(’a(l) * exp( a(2)*TimeMatrix )‘); Profile = sin(linspace(0,6*pi,m)); t = linspace(0, T, n*m); °/ TimeVector [MDsignal, TimeMatrix] = MkMDsig(MDfun_name, t,... [MDamp WI, [n,m]); CDsignal = MkCDsig(CDfun_name,TimeMatrix, [CDamp tau], Profile); h CDsignal(n/2+1:n,:) = bumpAmp*CDsignal(n/2-i-1:n, :); randn(’seed’ ,3133); dwt = MDsignal + CDsignal + sigmal*randn(n,m); elseif TestData 3 0 I0 °h CDcos + MDsin + MDstep 0 l0 T = 1000; % time in seconds n = 58; % number of scans m = 120; °/ number of sensors MD amp = .3; °h MD amplitude w = 2*pi/T; % MD omega CDamp = .3; °h CD amplitude tau = 0.00; °h CD time constant sigmal = 0.2; 7 noise level MDfun_name = sprintf(’a(l) * sin( a(2)*TimeVector )+ [zeros(1,a(3)) 2*ones(1,a(4))] ‘); CDf un_name = sprintf(’a(l) * exp( a(2)*TimeMatrix )‘); Profile = cos(linspace(0,2*pi,m)); t = linspace(0, T, n*m); 7. TimeVector [Mosignal, TimeMatrix] = MkMDsig(MDfun_name, t, [MDamp w (n*m/2) (n*m/2)], [n,m]); 73 7.1. The Main Program CDsignal = MkCDsig(CDfun_name,TimeMatrix, [CDarnp tau] ,Profile); randn(’seed’ ,3133); dwt = MDsignal ÷ CDsignal + sigmal*randn(n,m); elseif TestData == 4 •1 to °h Moisture data created using MVM 0/ to n = 54; °/ number of scans m = 120; °fl number of sensors a = 0.8; A = [1 0;0 a]; °fl State transition matrix. B = 0.4; = 0; °h Basis weight (B=0) = {—ones(1,m/2) ones(1,m/2)]’;70Profile. P = sin(linspace(0,6*pi,m)’); ubar = 0.5; SetP = 0; B. = 0.001; °/ Covariance of measurement noise. q = (0.1Y2; °h Variance of disturbance. Uhatacc = zeros(m*m,1); °h Accumulator for Ubar + Uk x = [0;0]; ‘Ia Initial State. I = eye(2); Yacc = zeros(m,n); °h Accumulutor for saving all Y’s CDsignal = zeros(n,m); tau = —0.03; ‘h all the equation numbers are taken from %lvar Mar Jonsson’s thesis 0 /0 °h Create the MD variation (Eq. 2.3 and 2.4) a, /0 randn(’seed’ ,3133); W = sqrt(q)*randn(n*m,1); U = ubar + filter({1],{1 —a],W); sqrtR = sqrt(R); noisel = sqrtR*randn(m,n); n_order_acc = zeros(n,m); n_order_f = 1:m; n_order_b = m:—1:1; 74 7.1. The Main Program n_order_acc(1:2:n,:) = n_order_f(ones(1,length(1:2:n)),:);n_order_acc(2:2:n,:) = n_order_b(ones(1,length(2:2:n)),:);for k=1:n n_order = n_order_acc(k,:); urange = ((k—1)*m+1):k*m; °h slow change in CD aa = exp(k*tau); CDsignal(k,:) = (P*aa)’; C = (ones(m,1) + B*P*aa)*[1 1]; ‘h Generate measurements. (Eq. 2.7) Y = CDsignal(k,n_order)’ + C(n_order,1).*U(urange)+...noisel(:,k); 7. Eq. 2.6 Yacc(:,k) = Y(n_order); end Yacc Yacc ÷ SetP; dwt Yacc’; elseif TestData == 5 0/ I, Y CDcos + MOsin + 2 x CDbump ÷ streak T = 1000; °h time in seconds n = 58; °h number of scans m = 120; number of sensors MDa.mp = .5; % MD amplitude w = 2*pi/T; °h MD omega CDamp = .2; % CD amplitude tau = 0.00; h CD time constant sigmal = 0.2; °I noise level streakPos = 60; streakWidth = 2; streakAmp = .6; bumpRaw = 30; bumpCol = 30; bumpHight = 10; bumpWidth = 5; bumpAmp = 1; MDfun_name = sprintf(’a(l) * sin( a(2)*TimeVector )‘); 75 7.1. The Main Program CDfun_naiue = sprintf(’a(l) * exp( a(2)*TimeMatrix )‘); Profile = cos(linspace(0,2*pi,m)); t = linspace(0, T, n*m); % TimeVector [MDsignal, TimeMatrix] = MkMDsig(NDfun_name,... t, [MDamp WI, [n,m]); CDsignal = MkCDsig(CDfun_name,TimeMatrix,... [CDamp tau], Profile); CDsignal (bumpRaw: bumpRaw+bumpHight-1,... bumpCol :buinpCol+bumpWidth-i)=... CDsignal (bumpRaw: buinpaaw+bumpHight- 1,... bumpCol : bumpCol+bumpWidth-1)+... bumpAmp * ones(bumpHight,bumpWidth); CDsignal(: ,streakPos:streakPos+streakWidth—1) = CDsignal(: ,streakPos:streakPos+streakWidth—1) + streakAinp*ones (n, streakWidth); bumpCol = 80; bumpRaw = 20; bunipHight 4; bumpWidht = 6; bumpAmp = -1; CDsignal (bumpRaw : buinpRaw+bunipHight-1, bumpCol bumpCol+bumpWidth-1) CDsignal (bumpR.aw : bumpRaw+bumpHight-1,... buxnpCol : bumpCol+bumpWidth- 1) +... bumpAmp * ones(buinpHight,bumpWidth); randn(’seed’ ,3133); dwt = MDsignal + CDsignal + sigmal*randn(n,m); end ‘I I0 function [MD, MD_exp, MD_way, CD, CD_exp, CD_wav]=... cdmdcalc (dwt , dwtExp, dwtNewF) 0 I0 °h Setup parameters narg = nargin; MaxArg = 3; if narg < MaxArg error ‘Not enough parameters!’ end Oj 10 76 7.2. The DCT Program ‘fl Calculation of CD and MD profiles 0 /0 [nl,ml] = size(dwt); n = 1:nl; m = 1:ml; °h Raw profiles MD = mean(dwt’); MD = MD(ones(1,ml),:)’; CD = dwt - MD; °h EXPO profiles MD_exp = MD; CD_exp = dwtExp(1:nl,1:ml) — MD_exp; °h Waviet profiles MD_wavl = mean(dwtNewF(n,m)’); MD_wavi = MD_wavi(ones(1,ml),:)’; CD_way = dwtNewF(1:nl,1:mi) - MD_wavl; MD_way = mean((dwt-CD_wav)’); MD_way MD_wav(ones(1,ml),:)’; 7.2 The DCT Program function [Ecd,Emd,Emd2]=DCT(y) Y0 Inputs y : The raw profile °h Outputs Ecd : The estimated CD scan Emd : The estimated MD scan(filtered by Wavelet) °h AUTHOR: (c) S.Karimzadeh °h DATE: 2008—09-09 I, 10 °h REVISION HISTORY X 2007—11—30: CREATION °h 2008—10—27: Modified with a filter for CD to separate °h controllable and uncontrollable variations °h M should be change to -M in line66 °h 2008—12—7: Last modified 77 7.2. The DCT Program °h Comments: °h In this code the PPc_main function developed by the ZoranX Nesic is used for wavelet filtering. The MATLAB’s wden °k function can be used alternatively. °h Finding the number of scans and number of headboxes[D,M]=size(y); %Setting the window length : ss is the number of scansY0in the window ss=32; /0NN is the number of data boxes in the windowNN=ss*M; °h Subtracting the mean of the raw profile. The values in °hy are saved for later use yy=y; for i1:D yy(i, :)=yy(i, :)—mean(yy(i,:)); end Y0initialization scdl=zeros(1, (floor(M*D/NN)+1)*NN); XChanging the 2D raw data to a 1D vector scdl(1,1:M*D)M2V(yy); rescdl=zeros(1 , length(scdl)); X=zeros(1,NN); duin=zeros(1,NN); % Processing the first ss number of scans °h finding the DCT coefficients dum(1,1:NN)=dct(scdl(1:NN)); Y0Keeping the coefficients which are integer multiple%of 2Ts X(1,ss+1:ss:NN)=dum(1,ss+1:ss:NN); h Reconstructing the CD Profile by taking %the inverse DCT rescdl(I. , 1 :NN)=idct(X); hfiltering to separate the controlable variationsY0 from uncontrolable °hvariations rescdl(1,1:NN)=PPC_main(rescdl(1,1:NN), 2, 3, 4,3,1,0);0/, dumm=zeros(1,NN); ‘h This loop processes the rest of the data by updating 78 7.2. The DOT Program °hthe window and °h estimating the CD profile for i=M:M:M*(D—ss)°fl °h resetting the dummy variables dum=zeros(1,NN); X=zeros(1,NN); XXXzeros(1,NN); dumm=zeros(1,NN); °h Finding the DCT coefficients for the current window dum(1,1:NN)=dct(scdl(1,1+i:NN-I-i),NN); °flxeeping the coefficients which are integer multiple Y0of 2Ts X(1,ss+1:ss:NN)=dum(1,ss÷1:ss:NN); XXX=X; °h Reconstructing the CD Profile by taking the Y,inverse DCT dumm=idct(XXX); °h filtering to separate the controllable variations °fl from uncontrollable °fl variations du!mn=PPC_nlain(dunun, 2, 3, 4,3,1,0); °h Copying the estimated signal in the final vector rescdl(1,NN+i—M+1:NN+i)=dumm(1,NN—M+1:NN); end ‘h Changing the 1D vector to a 2D matrix Ecd=V2M(rescdl(1,1:M*D),D,M); h Estimating the ND variation by subtracting the estimated %CD profile from the raw scan and compensating for the 7.scanning direction Emd= scanrev(y—Ecd(1:D,:),M,D); Y.filtering the MD variation by a wavelet filter Emd=PPC_main(Emd, 2, 4, 4,3,1,0);’ ‘h Estimating the MD variation by subtracting the estimated /DCD profile from °/ the raw scan ‘I. Filtering the MD variation by a averaging Emdmid=mean((y—Ecd)’); Emdmid=Emdmid(ones(1,M),:)’; 1For codes of the function “PPCmain” refer to {1O. 79 7.2. The DCT Program Emd2= scanrev(Emdmid,M,D); end 80 Chapter 8 Appendix B - Matlab Programs for Simulations in Chapter 5 °hPicketing Detection Simulation version 3.0 °h(c) Soroush Karimzadeh °flDate: 12/08/2008 dc ciear close all °hLoading the data to the workspace °hdata should be in the same folder as this mat file.Yo Loading the actuator inputs load autoslice.mat °h Loading the CD measurments load dryweight . mat °hextracting the data and converting from structure to variableuw=autoslice . setpoints_out_prof ile; s=uw’; °hfinding the dimensions of the process °h N is the number actuators °h D is tghe number of CD bins [M,D)=size(s) °hPioting the actuators input figure imagesc (s’) xlabei(Actuators’,’interpreter’, ‘latex’)ylabel(’Scans’ ,‘interpreter’ ,‘latex’) set(gca, ‘YDir’,’normal’) title(’Actuators Profile for Basis Weight’,’interpreter’... ,‘latex’) 81 Chapter 8. Appendix B - Matlab Programs for Simulations in Chapter 5 colorbar( ‘WestOutside’) figure °hPloting the measurments clim= [89 94] imagesc (dryweight . cd_bin_profile, dim) xlabel(’Sensors’, ‘intereter’,’latex’)ylabel(’Scans’ ,‘interpreter’, ‘latex’) set (gca, ‘YDir’, ‘normal’) title(’Raw Profile for Basis Weight’,’interpreter’, ‘latex’)colorbar(’WestOutside’) figure OfO) 0/0/0)0)0/0/0)0/0/0,00)0/0/0/0/0) 0/0/ 0/ 0/0/ O/0/0/0/0/0)O/0/O/O/O/ 0/0/0/0)0/0/0/0/0/0)0/0/0/0/0/0/0/0) 0l0/o/0l0lOl I, I, to to lo /otolo/oto/o/o!o 1* to lob do do to/of.!,!, 1010 (olO (ole!, do 1. 1, (0 (o (ole ( lo 1,/ole!,!, le!oboIo blob bo Io lob boTh% This loops finds the power of the fourier transform of °fl each scan for i=1:D °h Finding the fourier transform of each scan%the mean of each scan is removed first to get rid of the °h DC component. sff=fft(s(: ,i)—mean(s(: ,i))); °fl Finding the power of the fourier transform of each scan °h sf has th power for the current scan °h sout has the power of all scans sf=sff.* conj(sff) /M; sout(: ,i)=sf; °h extracting the power of the highest frequencyp(1,i)=sf(22,:); pH=sum(sf (17:22,:)); °h* (sf(8:22,:)’); pP=sum(sf(1:22,:)); Pr(1 ,i)=pH/pP; end Pr(1,1:2)=0; 0/0/0/01010/0/0)0/010/010/0/0/0/0/0)0) 0)0)0)0/0/0/0/0/0/0/0)0/0/ 00/ 0/0/O/0/O/O0/O/0/O/0/0/0)O/0/O)o) 0J0/0/O/0)0/0)0l lob I. lob to I to to to to to to bo to to lob Io /0/0/0/0 lo to to I, /0 bo /0 /0 to 1010/0 blob lob to/oh/o /0 I, blob!. (o/e (0 /o/o!oThlo/o °hplotting the Power of all scans stored in soutimagesc ([1 D] , [0 0.5] ,sout(1:23,:)) ylabel(’Normalized Spatial Frequency’) xlabel(’Scan ‘) colorbar DDD=zeros (size (sout)); DDD(17, :)80; figure 82 Chapter 8. Appendix B - Matlab Programs for Simulations in Chapter 5 clim=[0 50] imagesc ([1 D],[0 3.2808],sout(i:23,:)) hold clim=[0 50] imagesc ([1 D] , [0 3.2808] ,DDD(i:23,:)) ylabel(’Normalized Spatial Frequency cycle/meter’) xlabel(’Scan ‘) colorbar °flTo normalize the power we find the average of the power °/4n the steady state and divide the power by this factor.0flIn this case the steady state is assumed to be between scan 50 and 90.Thus the average of the highest °hfrequencie’s power between scan 50 and 90 are stored °hin “Cosm”. figure Cosm=mean(p(i,51:89)); °hplotting the normalized power for the highest frequnecyplot(p’/Cosm) 0/ xlabel(’Sensors’,’interpreter’, ‘latex’) °h ylabel(’Scans’ , ‘interpreter’,’latex’) °h set(gca, ‘YDir’ ,‘normal’) Y0 title(’Raw Profile for Basis Weight’,’interpreter’... ,‘latex’) 0,t0J000/ 0/0/0/0/0/0/ 0/0/0)0/Oj 0/0/0/0/0/0)0/0/0)0)0/Oj 0/ 010101 01 ) 0/0/0/0/0/0/0/0/Of 0/0/0/Of 0/0/Of 0/0/0/0/0/0/0/0//0/a la/a fo/ol,lofo to/a /o/a/a ,‘o (a (0 (0 10 ía (a Ia Ia blob (ala I. lab blob lola/a blob ía bo Ia/a / lob /0 Ia lola/a lola blob °hCUSUM Algorithm %initializing the S function S=zeros(1,D); °hfinding the mean of process at steady state. °flSince the process is nomarlized , mu will be equal to 1.mu=mean(p(i,51:89)/Cosm); °himplementing the CUSUM ‘h S is the S variable in CUStJM algorithm °h G is the g function m=S (1, 1); for i°2:D S(i,i)°=S(i,i—1)+(p(1 ,i—i)/Cosm)--mu; if S(1,i)<m mS(i,i); end g(i,i)=S(1,i)—m; 83 Chapter 8. Appendix B - Matlab Programs for Simulations in Chapter 5 end °kPlotting the POWER of the highest frequency%The red line shows the mean of %process in steady state(=i). plot(i:D,p/Cosm, ‘b’ ,i:D,mu, ‘r’, ‘LineWidth’ ,2)legend( ‘$P_{\frac{i}{2X_a}}$’ , ‘Average of P’); set (legendO, ‘interpreter’ ,‘latex’); title(’Normalized Power of the Highest Frequency$ ‘,... ‘interpreter’ ,‘latex’) axis({i D 0 200 ]) xlabel( ‘Scan Number’, ‘interpreter’ , ‘latex’) °hPlotting the S function °hfigure °hplot (S) ‘fltitle(’s’) °hPlotting C Function Y.Threshold is set to 5.The green line shows the thresholdfigure subplot(2,1,1) plot(1:D,g, ‘b’ ,i:D,5, ‘r’) legend( ‘g for $P_{\frac{i}{2X_a}}$’, ‘Threshold’); set (legendO ,‘interpreter’, ‘latex’); axis({1 D 0 2000]) xlabel(’Scan Number’,’interpreter’, ‘latex’) subplot(2,i,2) plot(1:D,g, ‘b’ ,i:D,5, ‘r’) legend( ‘g for $P_{\frac{i}{2X_a}}$’ , ‘Threshold’); set(legendO,’interpreter’, ‘latex’); axis([i 300 0 30]) xlabel(’Scan I’Tumber’ ,‘interpreter’, ‘latex’)0/0/0/0/0/0/0/0/0 0/0/0/0/0/0)0/0) 0t0/ 0/0)0/0/0/0/0/0/0/to to to fob to lob !olobolo to to to to to to to to to to /0 to 10 to to to °himplementing the CUSUM 7, S is the S variable in CUSUM algorithm °h C is the g function S2=zeros(i,D); Cosm2=mean(Pr(i,51:89)); mu2=mean(Pr(i ,51 :89)/Cosm2); m2S2(1,1); figure plot(i:D,Pr/Cosm2, ‘b’ ,1:D,mu, ‘r’, ‘LineWidth’ ,2) 84 Chapter 8. Appendix B - Matlab Programs for Simulations in Chapter 5 legend( ‘$P_r$’, ‘Average of $P_r$ ‘);set (legend 0,... ‘interpreter’ ,‘latex’); title(’Normalized Power of the Highest Frequency ‘,... ‘interpreter’ ,‘latex’) axis([1 D 0 8]) xlabel(’Scan Nber’,’interpreter’, ‘latex’)for i=2:D S2(1,i)=S2(1,i—1)+(Pr(1,i—1)/Cosm2)—mu2; if S2(i,i)<ni2 m2=S(1,i); end g2(1,i)=S2(1,i)—m2; end °hPlotting G Function hThreshold is set to 5.The green line shows the thresholdfigure subplot (2 , 1, 1) plot(1:D,g2, ‘b’ ,1:D,20, ‘r’) legend(’g for $P_r$’,’Threshold’); set(legendO,’interpreter’, ‘latex’); axis([i D 0 2000]) xlabel(’Scan Nber’,’interpreter’,’latex’) subplot(2,1,2) plot(i:D,g2, ‘b’ ,1:D,5, ‘r’) legend(’g for $P...r$’,’Threshold’); set(legend() ,‘interpreter’, ‘latex’); axis({1 300 0 60]) xlabel(’Scan Nber’,’interpreter’, ‘latex’) 85
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Online detection of picketing and estimation of cross...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Online detection of picketing and estimation of cross direction and machine direction variations using… Karimzadeh, Soroush 2008
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Online detection of picketing and estimation of cross direction and machine direction variations using the discrete cosine transform |
Creator |
Karimzadeh, Soroush |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | A scanning sensor at the dry end of a paper machine samples a noisy mix ture of both cross direction and machine direction variations. Although many techniques have been developed to separate these two variations and estimate a filtered profile, industrial practice has changed little over the years. In this work, a novel online Cross Direction (CD) /Machine Direction (MD) separation approach is developed. The method uses a fundamental prop erty of the CD variations to separate them from the MD variations. It is shown that the scanned CD variations build an even periodic function in the time domain and appear only as integer multiples of twice the scanning frequency in the frequency domain. Based on this property, the Discrete Cosine Transform (DCT) is utilized to separate the CD profile from noise and the MD variations. The performance of this method is verified through comprehensive simulation tests and is compared with existing wavelet and exponential filters. It is shown that the suggested method has better per formance in simulation scenarios. This is further validated by applying the method to industrial data. In the second part of this thesis, a novel picketing detection method is devel oped. The picketing pattern is associated with growing components in the spectrum of actuator profiles. A conventional change detection algorithm, CUSUM, is applied to detect these growing components. The method is verified in two simulation scenarios. It is shown that this method can detect and predict a picketing pattern before any picketing pattern is visible on the raw profile. |
Extent | 2576568 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067190 |
URI | http://hdl.handle.net/2429/7565 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_karimzadeh_soroush.pdf [ 2.46MB ]
- Metadata
- JSON: 24-1.0067190.json
- JSON-LD: 24-1.0067190-ld.json
- RDF/XML (Pretty): 24-1.0067190-rdf.xml
- RDF/JSON: 24-1.0067190-rdf.json
- Turtle: 24-1.0067190-turtle.txt
- N-Triples: 24-1.0067190-rdf-ntriples.txt
- Original Record: 24-1.0067190-source.json
- Full Text
- 24-1.0067190-fulltext.txt
- Citation
- 24-1.0067190.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067190/manifest