Cross directiona l contro l o f basi s weigh t o n pape r machine s usin g Gra m polynomial s by Kristinn Kristinsson B.Sc. University of Iceland, Reykjavik, 1986. M.A.Sc. University of British Columbia, Vancouver, 1989. A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA January 1994 © Kristinn Kristinsson, 1994 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. Department of Electrical Engineering The University of British Columbia Vancouver, Canada Date: January 1994 Abstract In this thesis a control algorithm for cross-directional basis weight control based on orthogonal polynomials is presented. Instead of controlling the cross-directional values at each point, the parameters of the orthogonal polynomials are controlled, taking into account the maximum bending moment of the slice lip and the relative stiffness of the slice with respect to the actuators. Problems with slice lip bending are dealt with in two ways. Implicitly, by modelling it as a polynomial such that the slice is not bent more than the highest degree of the polynomial, which is chosen beforehand. Explicitly, by calculating the bending moment and then use quadratic programming to limit the bending moment. By modeling the measurements and the actuators positions by polynomials the mapping between the actuators and the sensors, which is the key element of the current control systems, is effectively eliminated. Computational requirements are greatly reduced by avoiding inverses of big matrices. H Table o f Content s Abstract ii List of Tables vi List of Figures vii Acknowledgments xiv Epigraph xv 1 Introductio n 1 1.1 Literature overview 4 1.2 Motivation 8 1.3 Contribution of this work 10 1.4 Outline of the thesis 1 1 2 Discret e orthogona l polynomia l approximatio n usin g Gram polynomial s 12 2.1 Orthogonal polynomials 12 2.1.1 Creating Gram polynomials 13 2.1.2 Gram goes to Legendre 20 2.2 Approximations by Gram polynomials 22 2.2.1 Autocorrelation 26 2.2.2 Akaike's Information Criterion 27 2.3 Frequency characteristics of Gram polynomials 30 2.4 Summary 32 iii 3 Mode l fo r cross-directiona l basis weigh t response 35 3.1 Paper machine 35 3.2 Basis Weight 38 3.2.1 Square interaction matrix 44 3.2.2 Nonsquare interaction matrix 45 3.2.3 Interaction matrix parameterized with Gram polynomials 47 3.3 Slice lip model 52 3.3.1 Bending limit 59 3.4 Summary 59 4 Identificatio n 61 4.1 Identifying the interaction matrix 61 4.2 Identifying the parameterized matrix 69 4.2.1 Limited bending of the slice lip 73 4.2.2 Estimation of G x including dynamics 77 4.3 Identifying the inverse of the parameterized matrix 77 4.4 Summary 84 5 Basi s weigh t Cros s Directiona l Contro l 85 5.1 Dahlin controller and integrator antiwindup 85 5.2 Control of basis weight 89 5.2.1 Removing polynomials from a profile 91 5.2.2 Control of heavy weight paper 97 5.3 Comparison to conventional control 98 5.4 Chaos 104 5.5 Summary 109 6 Summar y an d Futur e wor k 113 6.1 Summary 113 6.2 Future work 114 iv References 115 v List o f Table s 4.1 Location of the PRBS inputs and the corresponding output measurements deduced from Figure 4.39 4.2 Number of data boxes between the maximum values of the cross correlations between the 11 responses in Figure 4.43 5.3 The points of bifurcation and the ratio of the differences between points of bifurcation. . . . VI List o f Figure s 1.1 A simplified paper machine 1 1.2 CD profile for industrial data and the average for 9 sections of the profile 4 2.3 Gram polynomials of first to sixth order according to equations 2.12-2.14 evaluated at the 100 points shown with the fine grid on the top x-axis 15 2.4 The world of Gram polynomials (first to third order) according to Equation 2.17, plotted for N=100 16 2.5 Gram polynomials according to Equation 2.22 for 100 points. Plotted for the first to sixth order 17 2.6 Plot of the first 30 Gram polynomials for 100 data points 21 2.7 Relationship between the number of data points and the order at which the maximum absolute value of the polynomial is greater than 2, 10 and 100. The circles in the lower left hand corner denote a fourth and an eight order polynomial for 20 data points 22 2.8 Gram goes to Legendre. The smooth solid line is the Legendre polynomial in both plots. . . . 23 2.9 Approximation of an industrial basis weight CD profile with Gram polynomials to the 10th order. The small inserted graph shows the 10 parameters of the Gram polynomial approximation 24 2.10 Measured data and 4, 10, 20 and 40th order approximation using Gram polynomials. . . . 25 2.11 Parameters of the Gram polynomials in Figure 2.10 and two times the standard deviation of the approximation error 26 2.12 Measured data and the measured data with 40th order approximation subtracted 27 2.13 Autocorrelation of a profile minus the Gram polynomial approximation for 4, 10, 20 and 40th order approximation. The 95% confidence interval is shown with the dotted lines. . . . 28 2.14 Alcaike's Information Criterion for the profile from Figure 2.9. The solid line shows the AIC, the other two lines show the AIC when it is modified to greater penalize model complexity. The circles in the plot denote the minimum for the curves 29 vi i 2.15 Left half of 100th order Gram polynomial for 4096 points. It is shown with the polynomial partitioned into 8 zones, i.e. the left half is partitioned into 4 zones, separated with different line types in the graph 31 2.16 Spectrum of 100th order polynomial 31 2.17 Spectra of 100th order poly, calculated for the foui sections of the left half of the polynomial shown in Figure 2.15 32 2.18 The original profile (dotted line) and the residual profile (solid line) when a 20th order Gram polynomial is removed. Power spectrum for the original profile (dotted line) and the profile when 20th order Gram polynomial is removed (solid line) 33 2.19 The original profile (dotted line) and the residual profile (solid line) when a 40th order Gram polynomial is removed. Power spectrum for the original profile (dotted line) and the profile when 40th order Gram polynomial is removed (solid line) 33 3.20 Paper Machine 36 3.21 O-frame 36 3.22 Response to a slice lip movement at one actuator position. The slice lip deflection is denoted by • with 0 denoting the deflection at a slice lip actuator. The + marks what the basis weight might look like 38 3.23 CD responses, according to Laughlin's interpretation from the literature , see also Equation 3.45 41 3.24 Reciprocal condition numbers of G, where it is created according to equation 3.43 with po = 1, pi = 1r and p2 = —r, for various sizes of the matrix. When G is badly conditioned, the reciprocal condition number is close to 0 43 3.25 Reciprocal condition number for a 105 x 105 matrix with po = 1, Pi = 2r and P2 = —r. Basis weight response for a machine with a singular interaction matrix, reciprocal condition number equal to zero 43 3.26 Dry weight measurements and the input to the corresponding slice lip location 44 3.27 Plot of the weight response to a straight line slice lip profile. The upper half is using G matrix from Equation 3.50 and the lower half is if the edges of the matrix are made smooth 45 3.28 Geometry when the number of actuators is not the same as the number of measurements. It is assumed that the first measurement is half a data box width from the edge 47 3.29 Meshed surface of the G x matrix for an interaction matrix of size 105x105 and po = 1. pi = 0.4 and p2 = 0.2 49 3.30 Every second diagonals of G x from Figure 3.29. The zeroth (main) diagonal is not shown 49 3.31 Reciprocal condition number of G x for a fixed number of Gram polynomials but varying the sizes of the interaction matrix, G, and its off diagonal elements which are pi = 2r and pi = - r 51 3.32 Reciprocal condition numbers of G x for G of size 105 x 105 with po = 1, pi = 2r, P2 — — r and varying the number of Gram polynomials 52 3.33 Plot of singular values of G x for an interaction matrix of dimension 105 x 105 with Po = 1> Pi = 2r, p2 = — r and with r = 0.3 53 3.34 Model of slice lip 54 3.35 Deflection of the beam in Figure 3.34 for different stiffnesses, when an actuator at x = 1 is moved one unit 56 3.36 Deflection of a beam when a actuator at x = 1 is moved to one from the initial position of zero (shown for the actuator as a dotted line with the endpoints denoted with small circles) 56 3.37 Slice lip 57 4.38 An experiment to identify the interaction matrix, G 62 4.39 Slice lip actuators and basis weight measurements from an industrial data set. To get a good contrast in the picture the mean value of every CD point is subtracted from the raw data 64 4.40 Slice lip set point and measured basis weight as a function of time, from an industrial data set 65 4.41 The cross correlation between the input and the output for the 11 PRBS signals from an industrial data set 65 ix 4.42 Estimates of steady state gain, pj, between the slice lip actuators and the corresponding measurements from Table 4.1 for an industrial data set 66 4.43 The estimated response across the sheet, for 11 actuators from an industrial data set. The plus signs mark the final estimates from Figure 4.42 66 4.44 The 11 estimated response shapes from Figure 4.43 shifted according to the first line in Table 4.2 such that they overlay 67 4.45 Average estimated response of the weight profile by moving a single actuator 68 4.46 Basis weight response to a flat slice lip (at 10/xm) with 105 actuators for two different models of a paper machine 70 4.47 Input, u c, and output, y c, used for the estimation of G x, when there are no limits on the slice lip bending 72 4.48 Overlay of 59 slice lip profiles, measurement profiles and bending moments corresponding to the inputs and outputs in Figure 4.47 73 4.49 Estimation of each row of G x when no limit is put on the slice lip bending. The true parameters for this simulation, are shown with dotted lines. . 74 4.50 Estimate of G x after 59 scans when no limits are put on the slice lip bending. The estimated matrix is denoted with * and dotted lines between its elements. The true matrix used in the simulation is denoted with o and solid lines 75 4.51 Input, u c, and output, y c, used to estimate G x, when there are limits on the slice lip bending 76 4.52 Overlay of 59 slice lip profiles, measurement profiles and bending moments corresponding to the inputs and outputs in Figure 4.51 77 4.53 Estimation of each row of G x when limit is put on the slice lip bending. The true parameters for this simulation, are shown with dotted lines 78 4.54 Estimate of G x after 59 scans when limits are put on the slice lip bending. The estimated matrix is denoted with * and dotted lines between its elements. The true matrix used in the simulation is denoted with o and solid lines 79 x 4.55 Input and output used to estimate the G x matrix when the dynamics is estimated as well 80 4.56 Estimation of the parameterized interaction matrix G x when the dynamics is estimated as well. The dotted lines show the true values 81 4.57 Estimation of a, for the 10 rows in the previous figure 82 4.58 Estimation of each row of KQ, the inverse of Gx. The dotted lines show the true values. . . . 83 4.59 The estimate of KQ after 59 scans. The estimated matrix is denoted with * and dotted lines between its elements. The true matrix used in the simulation is denoted with o and solid lines 84 5.60 Dahlin controller for a plant 86 5.61 Dahlin controller with control action taken every 3 scans with the closed loop pole, p = 0.05 and the open loop pole is, f p — 0.1. Dahlin controller with control action taken every scan with the closed loop pole, p = 0.37 and the same open loop pole. The * denotes the output at every third scan 87 5.62 The signals in a closed loop control with a Dahlin controller; with and without limits on the control action, and with and without integrator antiwindup. The plant pole is fp = 0.3 and the desired closed loop pole is p = 0.05 89 5.63 Controller with antiwindup 89 5.64 Block diagram of basis weight control strategy using Gram polynomials 90 5.65 Block diagram of basis weight control strategy using Gram polynomials with an integrator antiwindup. The controller, C, from Figure 5.64 is A (-i\ 90 5.66 Input, u c, and output, y r. The system is initially in open loop but the loop is closed at scan 40 92 5.67 Estimation of each row of G x. The true parameters for this simulation, are shown with dotted lines. The system is initially in open loop but the loop is closed at scan 40 93 5.68 Slice lip profile, bending moments and output at scan 100 94 5.69 Slice lip profile, bending moments and output at scan 100, when no limits are used for the slice lip 94 xi 5.70 Dry weight measurements for an industrial data set, and its polynomial approximation. . . . 95 5.71 Slice lip actuators settings 96 5.72 Final weight profile 96 5.73 Bending moments 97 5.74 Response to a slice lip movement for a heavy grade paper 98 5.75 Input and output for a heavy grade paper machine 99 5.76 Estimates of each row of G x for a heavy grade paper machine 100 5.77 Final estimates of G x for a heavy grade paper machine 101 5.78 Final basis weight profile, the desired profile, and the initial profile. The slice lip deflection. The bending moments 101 5.79 Block diagram of a conventional control 102 5.80 Conventional control of basis weight. The dotted line in the top graph denotes the initial profile, and the solid line denotes the final profile 102 5.81 Control with Gram polynomials 103 5.82 Disturbance rejections using conventional control on a well conditioned interaction matrix 103 5.83 Disturbance rejections using conventional control on a badly conditioned interaction matrix 104 5.84 Meshed surface of the G x matrix for G of size 105 x 105 with po = 1 and p\ — 0.4. . . . 106 5.85 Estimates of po and p\ 106 5.86 The desired profile and 10000 times the difference between the profile and its approximation with 10th order Gram polynomials. The smaller graphs show the parameters of the approximation for two sets of the Gram polynomials 107 5.87 Meshed surface of the G x matrix for G of size 105 x 105 with po = 1, p\ = 0.4 and uniformly distributed random numbers between -0.2 and 0.2 added to po and p\ 108 5.88 Meshed surface of the G T matrix using the Gram polynomials from Figure 2.3 109 xii 5.89 Meshed surface of the G T matrix for uniform G using the Gram polynomials from Figure 2.3 110 5.90 A bifurcation diagram and pseudo phase plane plots for different values of p I l l 5.91 Response of the system with p = 0.55, for 1% difference in initial condition on the desired value of the CD value at location 10 112 xiu Acknowledgments I would like to use this opportunity to thank my research supervisor and colleague, Prof. Guy A. Dumont, for his advice and guidance during the course of this work. I would also like to thank Prof. Michael Davies and other past and present members of the Paper Machine Group at the Pulp and Paper Centre for their valuable discussion. Thanks to Dr. Oliver Lee, Leoncio Estevez-Reyes, Navid Mohsenian, Dr. Ye Fu and Jahan Ghofraniha for reading parts of the thesis. Also thanks to the students and the staff at the Pulp and Paper Centre at the University of British Columbia for all their help and assistance, and for making my stay here very enjoyable. Acknowledgment is made to Paprican's Frank Leslie Mitchell memorial fellowship, for financial support, and to Jons Th6rarinssonar memorial fellowship. Special thanks to friends and family in Iceland that have helped in one way or another during the course of my study. Last but not the least my heartfelt thanks to Anna D6ra for her encouragement when I was working on my thesis. xiv / know why there are so many people who love chopping wood. In this activity one immediately sees the results. Albert Einstein xv Chapter 1 Introduction The three R's (reduce, reuse and recycle) have been catching on recently as a means of reducing the consumption of raw materials. Paper, which in Canada is manufactured out of wood chips, is one of the prime materials being targeted by the three R's. One way of reducing usage of raw materials is in the manufacturing of paper, where a single percent reduction in fiber usage translates to preservation of some 800 km2 of forest per year [79], which is about 2.5% of the area of Vancouver Island. Not only can we reduce the amount of paper we consume, but we can also reduce the amount of fibers needed to manufacture a sheet of paper. The weight of paper per unit area, or the basis weight (BW), is a fundamental physical property of paper [30]. Other important properties are moisture and caliper. These three properties are interdependent; changing one will affect the rest, but it has been observed that the oven dry basis weight profile is dominant [69]. When basis weight is expressed in g/m2 or mass per unit area, it is referred to as grammage. Grammage variation is synonymous with substance variation and basis weight variation [32]. Basis weight is the most commonly used term and will therefore be used throughout this thesis when referring to mass or weight per unit area of manufactured paper. MD - 4 [ Press _ . . 4 % CD / ' Wire n ' 1 o Calender ^ \ I , "- " Dryer s / Ree l Figure 1.1: A simplified paper machine. Chapter 1; Introduction The paper machine (see Figure 1.1), invented almost 200 years ago [68], has been constantly improved over the years. It has become wider and faster, and modern paper machines make hundreds of tons of paper per day. With the increase in their width and speed, the problem of controlling them has become more difficult. The goal of control is to make the paper as uniform as possible. Variations will undoubtedly exist in the manufactured paper sheet, because of variations in raw materials and the complicated pulping process prior to the paper machine. To better understand the source of these variations, and ultimately to remove them, it has become necessary to classify variations in the paper sheet based on their source or how they might be reduced. Variations are generally divided into 3 components, based on their source, following Burkhard's and Wrist's paper in 1954 [16]: — Machine Direction (MD) part. Variations occurring along the sheet travel on the paper machine and affecting the whole width of the machine. — Cross Direction (CD) component. Variations occurring perpendicular to the sheet travel on the paper machine. The CD variations generally change much slower than the MD variations and can be assumed to be nearly time-invariant or stable with time. The CD variations will therefore persist in successive profiles. — Random (R) element. Part of the variations that occurs neither along nor across the machine. This part is also referred to as the residual part, because variations that cannot be attributed to either CD or MD, are assumed residual or random. An alternative classification scheme has appeared in recent times as measurement suppliers have started categorizing basis weight variations according to how they are controlled. Variations that can be controlled by CD actuators are called CD, while those controllable by MD control loops are called MD [30]. The different definitions can be attributed to the fact that variations that appear in the machine direction can be better controlled by actuators in the cross direction — they have faster response time and shorter time lag since they are closer to the measurement point. Researchers have also recently started looking at ways to represent the paper sheet as a two-dimensional process in which the CD and MD are two mutually-dependent dynamical processes [123]. The goal of control is to increase uniformity of the manufactured product. A reduction in variability influences nearly every factor which tends to increase production and improve paper quality [16]. Many benefits can be gained from reduced cross direction fluctuation, e.g. the average basis weight can be set lower and the average moisture can be set higher, resulting in a reduction in raw material usage. A more uniform profile also produces a stronger web and a more uniform reel. 2 Chapter I: Introduction leading to a reduction in the number of sheet breaks on the paper machine and in the converting process. It has been argued that variations in CD are more important than variations in MD or random variations [30, 31]. Since the manufactured paper is wound into rolls the MD and random variations will average out for each roll. The paper roll is usually cut into few smaller rolls, the CD variations will therefore be particular to a specific roll. This can be observed from Figure 1.2 where a CD profile (dotted line) along with its average for 9 sections of the profile (solid line), are shown. Assuming that the MD and random variations average to zero for the production of a roll of paper, and that the CD profile will be present in the whole roll, the 9 rolls will all have different averages. It is stated in [26]: "it is a fact that changes across the machine are always greater than they are in the machine direction". That statement is supported in several references [17, 125], however, it may be due to the fact in most cases the MD is already controlled [50]. Care must be exercised when using variance to judge the performance of a CD controller. The reduction might, in part, be due to low profile resolution [95]. The measured profile mapped to the actuators might even show zero variance whereas the unmapped profile has some streaks [103]. It can, for example, be shown that the expected value of the standard deviation for normally distributed random numbers reduces by 1/V2 if the random numbers are averaged over two data boxes. A better technique for variability analysis makes use of profile auto-spectra [96], or power spectra. Fourier analysis and power spectra can be used [74] to show that the part of the variance with wavelength bigger than or equal to twice the size of the actuator spacing can be reduced, whereas it is not possible to reduce shorter wavelengths. However, it has been reported for a industrial application that wavelengths slightly shorter than two actuator spacings can be reduced when the controller makes use of full resolution measurements [63]. One of the major problems with current CD control systems is to correctly align actuators and measurements. Inaccurate mapping is mentioned in [4] as one source of poor controller performance. Misalignment of even a few inches can result in sawtooth profiles, since actuators are adjusted to compensate for errors outside their control zone [89]. For optimal control results, measurement zones should align with control zones. When the measured zone on the dry end, corresponds to the correct actuator zone on the headbox, performance is maximized [17]. If they do not align, the data must be transformed, with some loss in profile resolution [89]. As we shall see, the problem of misalignment between actuators and measurements can be overcome through the use of orthogonal polynomials. 3 Chapter 1: Introduction Figure 1.2: CD profile for industrial data and the average for 9 sections of the profile. 1.1 Literatur e overvie w There have been several review articles on the control of pulp and paper processes [76, 117], as well as several surveys on the use of CD control in the industry [108-110]. In [36] a survey is given on the use of advanced control methods in the process industry in general, and [35, 37] give a survey on the use of advanced control in the pulp and paper industry and adaptive control in papermaking in particular. The three factors that have governed the industrial implementation of CD control systems are: sensors, actuators and process control computers [121]. With the introduction of the microcomputer in the early 1980's, interest in CD profile control increased, indeed the 1980's have been called the decade of cross direction profile control [106]. With faster computers, the profile could be measured in smaller and smaller data boxes, which consequently called for smaller and smaller actuator spacings. Now it is common to have three to five measurements for each actuator for the control of basis weight. On the other hand, one early investigation on CD basis weight control assumed that there were more actuators than measurements when automatically controlling the basis weight profile [116], but that is clearly inadequate, as it is not possible to control variations that are not visible in the measurements. Traditionally the measured profile was mapped to the actuators before calculating the control action, but recently, use has been made of the high resolution profile when calculating the control action 4 Chapter I: Introduction [63]. The idea here is that the controller will control high frequency variations better if it knows they are present. The controller did reduce wavelengths slightly shorter than two actuators spacing. As mentioned above, Burkhard and Wrist [16] were the first to classify basis weight variations into three components. Initial efforts at reducing variations concentrated on the control of MD variations [5], mainly because CD sensors were not available and more computer power was required for control of cross directional variations. CD control first became possible with the invention of sensors to measure the paper across the sheet. Initial efforts were mainly focused on troubleshooting machine faults rather than on control. In 1957, Cuffey [28] obtained a picture of basis weight profiles with an optical basis weight scanner. He used a bank of lamps over the moving paper sheet to illuminate the sheet; a rotating mirror reflected the light to a lens and then to a photosensor. With this setup it was observed that a movement of a single slice lip actuator caused the basis weight to change not only at that location but also at the adjacent locations. It was also observed that the slice lip opening changes at adjacent actuators when a single actuator is altered. In 1967, Astrom published one of the early papers on computer control of paper machines [5]. This paper was on machine direction control of basis weight and summarized a project, dating back a few years, to implement computer control in industrial settings. A minimum variance controller was designed by use of linear stochastic control theory. The beta-ray sensor utilized was able to traverse the paper web but it stayed at a fixed position. An application of adaptive control to moisture and basis weight in the machine direction was reported in 1975 by Cegrell and Hedqvist [19]. In 1969, about 10 years after Cuffey's paper, Dahlin presented an algorithm based on exponential filtering, to extract the MD and CD component from raw data, and gave a description of the sensors used for basis weight and moisture measurements [33]. Around the same time Beecher and Bareiss reported on their efforts to use a closed-loop analog controller on the slice lip for basis weight control; the inverse of an interaction matrix was used to describe how the measurements and actuator settings were related [8]. They had 16 motorized slice actuators on a 2.3 m wide machine capable of speeds between 0.76 m/s and 3.81 m/s. Standard deviation from the mean of the basis weight profile was reduced. Better results were obtained a couple of years later when the system was made digital [44]. With installation of a titanium slice, further reduction in basis weight variability was observed [26]. In 1975, Carey and co-workers [17] at Kodak reported a reduction in fibre usage following installation of a CD control system, which used the interaction matrix approach. Two years later Boyle [13] showed that the interaction matrix approach could result in excessive bending and/or displacement of the slice lip and that the solution may well be infeasible, calling for actuators to 5 Chapter 1: Introduction be more than full open or less than fully closed [14]. Two methods emerged to limit bending and displacement of the lip. The first method, quadratic programming, minimized variations subject to hard limits on the actuators movements; however, it was found to be computationally intensive. The second method, minimization of quadratic objective function, penalized actuators moves or bending of the slice. In 1977, Boyle [13] described a quadratic programming formulation for simulation and noted that computational considerations limit the applicability of the method to only a few actuators. The following year, the same author [14] revealed a simulation study, which used a quadratic objective function for a system with unequal number of actuators and sensors. A drawback with this method is that a matrix of size number of actuators x number of measurements has to be inverted or its pseudoinverse found, although this task can be somewhat simplified by taking advantage of the symmetry of the matrix. A particular problem with stiff slice lips according to Garnett [52] is that the lip cannot make the moves necessary to control some profiles. This concern was repeated a year later by Boyle [15]; the author suggested that tests should be carried out on slice actuators, prior to purchase of a CD weight control system to evaluate the flexibility of the slice lip. Cutshall and Hamel [31] give an overview in 1988 of some factors that need to be evaluated when slice lip automation is considered. In Wilkinson and Hering [127], a linear quadratic optimal design using on line identification of the interaction matrix is used in industrial settings. This work, which was partly based on earlier work by Richards [100] and Wilhelm [124] from 1982, gave between 44—55% reduction in basis-weight variations on several installations. The same year, Wilhelm and Fjeld [126] reviewed some control algorithms for cross directional control. They fonnulated control methods for both the inverse interaction matrix approach and for quadratic optimal design. They noted that strong coupling in the interaction matrix could result in the matrix being singular and therefore not invertible. They suggested that strong coupling simply represented a poor choice of actuator spacing. The problem could be avoided by controlling groups of actuators together, or computing the control based on a wide spacing of actuators and then using interpolation to calculate the settings of the remaining actuators. This work was extended three years later by Chen and Wilhelm [24] and Chen, Snyder and Wilhelm [25]. They looked at both quadratic programming (QP) and quadratic penalty function (QPF) methods; the latter is a form of quadratic optimal design. The QP gives an optimal solution, however, time and memory requirements for the QP method prohibit real-time control. Therefore, they turned to QPF, which gives a suboptimal solution. Generalized least squares was used to identify a spatial non-causal model of the process and more than 50% reduction in standard deviation of BW 6 Chapter I: Introduction was achieved with an adaptive control algorithm based on QPF [25]. This algorithm was shown to be equally appropriate to manage the wide response model for heavy weight machines [21]. Later, it was used as a basis for controlling a single CD profile with multiple actuators and for using the CD actuators to assist in MD control [57]. Interactions in CD have been a cause for great concern in basis weight control. The interactions depend on the operating point of the paper machines, from which the need to identify the interactions follows. One such identification scheme is described in a patent granted in 1991 [22] for determining cross direction responses and mappings of measurements onto actuators used for control. The responses were determined by means of PRBS signals and by calculating the covariances between the input and the outputs. The perturbing signals were defined by PRBSs with multiple signals being selected to be statistically independent [22]. In the early 1980's researchers with STFI in Sweden developed and installed a CD control system in a mill. Karlsson et al. [72, 73] observed a 20% decrease in standard deviation for basis weight using an undisclosed control method. Some of the highlights of the control package were given by Elliot in 1984 [45, 46]. The algorithm was used in several installations to reduce CD variations [95, 96, 105]. Hansson et al. [59] have also reported on the use of the STFT control package, where a slice lip was used to control moisture on a kraft machine. The paper machine had 48 slice actuators and a width of 7.1 m, with slice opening at 50 mm. In 1987 Bergh and MacGregor [9] employed LQG (linear quadratic gaussian) control theory to jointly control MD and CD moisture variations using spatially distributed actuators. They modeled the disturbances for a sheet forming process as a multivariate time series. A multivariate LQG control with recursive identification to sequentially adjust slice lip actuators with a traversing robot is given in Halouskova et al. [58]. None of the control strategies mentioned so far have addressed the issue of error in the interaction model. An exact model will never be possible because of wave propagation on the wire. In 1988, Laughlin [79, 81] came up with a MIMO (multiple input multiple output) robust control scheme based on internal model control (IMC) [80]. He was able to guarantee both robust stability1 and robust performance2, by assuming knowledge of the parameter uncertainty. A major drawback with this scheme is that the gain matrix describing the relationship between actuators and sensors has to 1 A closed loop control system is said to be robustly stable over a set of plants if and only if it is stable for all plants in the set. 2 A closed loop control system is said to exhibit robust performance if for all frequencies the sum of weighted sensitivity and complementary sensitivity functions is less than one. 7 Chapter 1: Introduction be positive definite. In [79], a table of reported matrices from the literature, shows that five out of thirteen matrices are not positive definite [79]. Furthermore the author does not take into account the constraints of the slice lip, so the control action from this scheme may not be realizable. There have also been reports about reduction in cross directional variation, an increase in production and/or reduction in raw material usage with the installation of slice lip control. A few of them can be found in references [10, 12, 47, 53, 34, 60, 82, 91, 119, 128]. Matsuda observed in 1990 [88] that increasing the number of measurements to twice the number of actuators improved control because narrow streaks could be eliminated. A year later Yamamoto [129] reported an application of fuzzy logic for profile control of a paper machine. The author reported using only 4 rules to eliminate the sawtooth form from the CD profile. A patent describing the use of a fuzzy logic controller on a slice lip was granted in 1992 [104]. In the patent an example is given of a paper machine with 360 measurements and 37 actuators. The fuzzy rules are derived from manual manipulations by a skilled human operator. Because of the complexity in controlling the cross direction dry weight Amyot, Nuyan, and Evans designed an expert system to help operators of a paper machine troubleshoot the cross direction basis weight control system [4, 48]. 1.2 Motivatio n The objective of this thesis is to develop a strategy for control of CD variations on paper machines which can address the shortcomings of current control methods, namely problems with interactions in the cross direction, and the need to map the measurements to the actuators. Since the basis weight is a fundamental property of paper, cross-directional basis weight control has been chosen for this study. The sampling rate on the CD basis weight loop is slow (around 30 seconds), so this is essentially a steady-state control problem where the purpose is to solve the decoupling problem between the numerous actuators [36]. Most current techniques for cross-directional basis weight control manipulate the cross directional values at each point separately from each other. Therefore, a large interaction matrix from the slice lip to the measurements has to be identified and its inverse found. The inverse matrix or mapping between measurements and actuators is the key element of current control methods and, if not correctly done, a sawtooth profile will result. Filtering is used in order to protect the slice lip from excessive bending. The filtering may take the form of sending the slice lip set points to a low pass 8 Chapter I: Introduction filter; of penalizing bending when calculating the slice lip profile; or of placing some hard limits on the movement of actuators. Another way to deal with excessive bending is to compute control settings of a few actuators and the remaining actuators settings are then found by interpolation using a polynomial of a certain order [126]. The actuators can also be treated in groups to simulate a change in the spacing of the actuators, since some paper grades might not be able to take advantage of close spacings [125]. Current industrial trend is towards an increased number of CD points. Measurement data boxes as small as 1 cm are reported by Taylor in 1991 [115]. For a 10 m wide paper machine this implies 1000 data boxes. To make full use of this increased measurement resolution, the number of actuators in the cross direction has increased as well. However, this gives rise to large interactions between adjacent CD points so that an actuator movement in one CD position affects many of its neighbors. Such interactions cannot be ignored when calculating a new control move and adjacent actuators must work not to exceed the bending limits of the slice lip they are controlling. Another problem associated with high resolution profiles is that in order for a CD control system to operate effectively, it is necessary to accurately determine which positions in the CD profiles will correspond to the locations of the slice screws. If this is not done, self-oscillation will result in the profile when a slice lip actuator tries to compensate for an error outside its area [95]. Recently, work has been devoted to the extraction of the cross direction (CD) and machine direction (MD) profiles from scanned measurements [20, 23, 39, 40, 70, 92, 122]. The main reason for extracting profiles from measurements is to exploit the information obtained for control. In this thesis, knowledge of the CD profile will be used for control. In this work, the CD profile or its shape is approximated by discrete orthogonal polynomials and then the parameters of the polynomials are controlled rather than the cross directional values at each point. The slice lip has certain limits on bending so its profile is a smooth curve that is naturally described by a polynomial of a given order. By modelling the slice lip with a polynomial, problems with slice lip bending are dealt with implicitly because the slice lip will not bend more than the highest degree of the chosen polynomial. It will also be shown in this thesis how bending of the slice lip can easily be limited by modifying the parameters of the polynomials. Fewer parameters are controlled and fewer parameters need to be identified because the degree of the polynomial is much smaller than the number of CD points. Computational requirements are greatly reduced by avoiding inversion of big matrices and inverses are found for a diagonal matrix. Finally, by modelling the measurements and actuators positions with polynomials, mapping between actuators and measurements, which is the key element of current control systems, is effectively eliminated. 9 Chapter 1: Introduction Modelling profiles by polynomials has already been used in the pulp and paper industry to control wet end moisture profiles on a paper machine through application of pressure to the ends and center of the press rolls [75]. The moisture profile was approximated by a set of polynomials and then the coefficients of the polynomials were controlled. Since only three actuators were used, second order polynomials were considered and consequently not much could be gained by use of orthogonal polynomials. Modelling profiles with discrete orthogonal polynomials has previously been used in Sendzimir steel rolling mills [41-43, 56, 101, 102]. Models for these plants can be derived from first physical principles and, therefore, analytical models exist. Only about 10 actuators and up to 31 measurements were used for cross directional control, so that only low order polynomials were required. Paper machines, in contrast, have complex interactions in CD that are difficult to model. Furthermore, the paper sheet is made up of small fibers that tend to flocculate, and are thus difficult to control. The number of measurements on paper machines is about 20 times that on steel mills and the number of actuators is about 10 times. Consequently, higher order polynomials are needed to describe the shape profile of manufactured paper. 1.3 Contributio n o f thi s wor k 1. Use of Gram polynomials to describe profiles on paper machines. Several forms of the Gram polynomials are investigated and it is shown that the scaling of the polynomials is important. Some of the frequency characteristics of the Gram polynomials are evaluated. 2. Use of Gram polynomials to describe slice lip deflection and bending moments. Problems with slice lip bending are implicitly dealt with by modelling it as a polynomial series, since the control method will not bend the slice more than the highest degree of the series. Slice lip bending is related to the parameters of the orthogonal polynomials describing the slice lip, thereby allowing the bending moments to be restricted, directly from the parameters of the orthogonal polynomial of the slice lip deflection. 3. Profile control design in the Gram space. Instead of controlling the CD values at each point, the parameters of the Gram polynomials are controlled, taking into account the maximum permissible slice bending. Because the number of polynomials is much smaller than the number of CD points (10 polynomials seem sufficient to describe most profiles in our industrial data), this technique greatly reduces computational requirements by avoiding inversion of large matrices and by using the fact that the polynomials are orthogonal. Furthermore, we effectively 10 Chapter 1: Introduction eliminate the need for mapping between the actuators and the data boxes, a crucial difficulty of current control systems. 1.4 Outlin e o f th e thesi s The discrete orthogonal polynomials used in this work and how they are used to approximate a profile or vector are discussed in chapter 2. A model is developed in chapter 3 to relate i) slice lip deflection to measured basis weight profile, ii) actuator setpoints to slice lip deflection, and iii) parameters of the discrete orthogonal polynomial approximation of the measured output to parameters of the discrete orthogonal polynomials of the slice lip. In chapter 4 identification of the interaction and the parameterized matrices is described. Chapter 5 presents a simulation of the result of applying the method described in this thesis to control a paper machine via a Dahlin controller. Finally, in chapter 6 some conclusions are drawn, and a future work is discussed 11 Chapter 2 Discrete orthogonal polynomial approximation using Gram polynomials Polynomial approximation is a technique used to obtain a smooth approximation to noisy measurements [99]. Instead of including measurement errors in the approximation, we seek to find an approximation which gives the general trend of the data without being too concerned with local fluctuations [90]. The use of orthogonal polynomials achieves certain computational simplifications. There are many different orthogonal polynomials that can be employed. For discrete measurements it is natural to choose discrete orthogonal polynomials. Some discrete orthogonal polynomials are appropriate for uniformly-spaced data points. Others can be used when data points are to be weighted. 2.1 Orthogona l polynomial s Let y(x) be a sequence of measurements observed at x = {x(i); i = 0,... ,N — 1}, where N is the number of data points. Let Pj(x) be a sequence of functions defined for every x(i) and j = 0 , 1 , . . . , m, where m is the highest degree of the polynomials. The plan is to write a function, y(x), as a linear combination of Pj(x): y(x) = yrX0)P0(x) + y c(l)P1(x) + y c{m)Pm{x) (2.1) or y(x) = [P0(x) Pi(x) PnAx) ydo) yAm) (2.2) with the parameters y c(j) to be determined such that the weighted square of the difference between the true (or measured) function and the approximated one is a minimum. In an equation form this is expressed as a least-squares problem: l 2 • r.v-i min < 2, w(x{i)) y(x(z))-^P,-(a;(i))y r0') (2.3) 12 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials To minimize the above expression, its partial derivatives with respect to yc(j) can be found and the result set equal to zero. This gives: A r - 1 J2 w (x(i)) y(x(i))-J2PMi))ycU) Pk{x(i)) = 0 (2.4) i=0 or with a little rearrangement: A r - l m N-l J ) w(x(i))y(x(i))P k(x(i)) = ^ y c(j) £ ™(:r(i))P,(*(*))Pfc(:r(z)) (2.5) j'=0 j=0 »'= 0 Calculation of y0(i) from the above equation is greatly simplified if the set of functions (Py(x)} is orthogonal over the set of points x(i). By definition a set of functions {P,(x)} is orthogonal over the set of points x(i) with respect to the weight function w(x) if A ' - l 53«»(x(t))P,-(a:(i))Pfc(x(i)) = { ° . $ ] t \ for j,k = 0,l,...m (2.6) »'=o Furthermore, for equally spaced data points and the particular case when w(x(i)) = 1, the discrete orthogonal polynomials are called Gra m polynomials [99 p.259]. They are sometimes referred to as discrete Legendre orthogonal polynomials (DLOP's) because it is observed that the numeric coefficients of the DLOP's are, apart from the sign and possibly a scaling factor, precisely the same as those of the shifted Legendre polynomials on the interval [0,1] [93]. In [66 p. 290] it is mentioned that the polynomials are sometimes referred to as Chebyshev polynomials, although that name is usually reserved for continuous polynomials with w(x) = 1/Vl — x2. 2.1.1 Creatin g Gra m polynomial s A recurrence relationship to create discrete orthogonal polynomials is given in [18, 51, 67, 99] and can be derived from the Gram-Schmidt orthogonalization of the powers l , x , x 2 , . . . x m [78]. The recursive formula in [67, 99] is: P m + i (x ) = (x - a m + i )P m (x ) - /?TOPm_i(x) (2.7) where "ft w(x(*)Mi){Prn(x(t))} 2 (Xrn+l = ' ~ A f _ 1 ( 2 - 8 ) E w(x(i)){Pm(x(i))} 2 i=0 13 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials and E 1 ' w(x(l))x(i)P m-i(x(j))Pm(x(i)) A E w{x{i)){Pm(x{i))f f3m = i ^ L _ = _|d> ( 2 . 9 ) E wixiiMP^ixii))} 2 E w(x(z)){P m^(x(i))}2 and the initial conditions for the recursive equation is given by P0(x) = l and P_i(x) = 0 (2.10) For an equally spaced x on the interval [—a, a] and for the particular case w(x(i)) = 1 for all i, «m+i = 0 so the recursive relationship becomes P m + i (x) = xPm(x) - /?mPm_!(x) (2.11) or Pm(x) = xPm_i(x) - /3TO_1Pm_2(x) (2.12) If the interval the polynomials are defined over is chosen as [—1,1], the polynomials will be evaluated at the equally spaced points: x(i) = 2 ^ - y - 1 ie[0,N-l] (2.13) and the parameter 0 m-i defined as E {Pm-l(x(i))} 2 {m _ 1}2 (N2 _ (m_ 1 } 2\ «'=0 The first few polynomials are i Pi{x{i)) = 2 N-\ i i2 i 2(N - 2) P2(x{l)) = 4W^f " 4^^+ ^^ (2 -15) 3 ?2 i 6iV2-15JV + l l 2(iV - 2)(JV - 3) P3(x(i)) = 8 =- - 12 5- + 4- 9 , V W ; ( i V - 1 ) 3 ( i V - 1 ) 2 ^ - 1 5 ( i V - l ) 2 5 ( i V - l ) 2 Figure 2.3 shows the Gram polynomials calculated using Equations 2.12-2.14 for 100 points (N=100). For clarity, these polynomials are shown as continuous lines, even though they are defined only at the discrete abscissas indicated with tick marks at the top of the figure. The maximum value of 14 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials Figure 2.3: Gram polynomials of first to sixth order according to equations 2.12-2.14 evaluated at the 100 points shown with the fine grid on the top x-axis. the polynomials gets smaller and smaller as the degree of the polynomials gets higher and higher. Accordingly, the parameters for the higher order polynomials can become very large. Another recursive scheme for calculating orthogonal polynomials for N evenly spaced measure-ments is [7]: m2(N2 -m 2) Pm+1(i) = P^Pxii) - 4 ( ; m 2 _ 1 } 'Pm-i(i) i = 0 , 1 , . . . , J V - 1 and m = 1,2,. . . , N - 2 (2.16) The starting conditions are Po(i) = 1 and Pi(i) = x(i) = i — ^-j^~- With the same initial conditions, the recursion can also be written as: ( m - l ) 2 ( i V 2 - ( m - l ) 2 ) Pm(i) = Pm-WP^i) r - ^ r J -Pm-2(l) 4f4(m-l)2- lJ i = 0,l,...,N -1 and m - 2 , 3 , . . . , N - 1 The first three polynomials are then written as: A r - 1 (2.17) Pi(i) = i -P2(i) = i2 - i(N - 1) + (N - 1)(N - 2) (2.18) P,(i) = i i - i 3 :1 3 ( i V - l ) ,6iV2-15iV + l l (N - 1)(N - 2)(N - 3) + i 10 20 15 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials Figure 2.4 shows the Gram polynomials calculated using Equation 2.17 for 100 points, and plotted as a function of x. Here, the maximum value of the polynomials gets larger and larger as the degree gets higher and higher. Due to the poor scaling of these polynomials, the scaling of the parameters here becomes even worse than for the first set of polynomials (see Figure 2.3). For example P3(0) « -N 3/20 = -1003/20 = -50000, whereas Pi(0) « -N/2 = -100/2 = -50. Poor scaling of the polynomials is apparent in Figure 2.4 where the first and the second order polynomials are barely visible because the third order polynomial has much larger values. x lO* Figure 2.4: The world of Gram polynomials (first to third order) according to Equation 2.17, plotted for N=100. Yet another scheme for creating discrete orthogonal polynomials is derived in [90]. For N equally spaced points, the polynomials are written as: pm(i) = j^(-iy 7=0 m J m + j\ vi* J ) (N - 1)(J) (2.19) The product i^ = i(i — l)(i — 2) • • • (i — j + 1), where j is a positive integer, is called the factorial jth power of i [66], the jth fading factorial [93] or the backward factorial function of order j [1]. Remembering that the binomial coefficients are written as ( "' J = jUm-j)'.' Eq u a t i ° n 2.19 can also be written as: (-iy (m + jfj) &) Prr E j = 0 0"! (N-l) or (2.20) 16 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials The orthogonality of the polynomials is not affected by multiplication of the polynomials with a scalar. The general form for the polynomials can therefore be written as [66 p. 289]: (-iy(m + j){2j) fi) Pm(i) = CmN ^ j=0 (iO2 (2.21) (N - l) (j) The N points, i, are denned over the interval [0, JV —1], and m is the polynomial order. The polynomials, with CmN = 1, from degrees one to four, inclusive, are: i Pi(i) = l - 2 Pz{i) = 1 - 1 2 P 2 (0 = 1 - 6 + 30 -+ 6 J V - 1 N-l (N-l)(N-2) »'(» ~ 1) on »(* - 1)(* - 2) (2.22) 20-N-l (N-l)(N-2) (N - 1)(N - 2)(N - 3) 140 — + 70-P4(i) = 1 - 2 0 - 4 — + 9 0 *(Z 1 } ' i V - 1 ' ~-(N-l)(N-2) ""(N-lf) (N-1)M Those polynomials are plotted in Figure 2.5 for 100 points and up to the sixth order, as a function of x. It can be seen that the six polynomials shown in the graph, all have the same maximum value and all the odd polynomials have the same minimum value. These polynomials are much better scaled than the previous two. Figure 2.5: Gram polynomials according to Equation 2.22 for 100 points. Plotted for the first to sixth order. 17 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials It can be shown that Equations 2.12-2.14 (the first form of the polynomials given) can be written as Equation 2.21 with (-l)mm\(N - 2 ) ( m - 1 >2 m - 1 CmN = • — ; 7—T— ( 2 . 2 3 ) ( A T - l ) r a - 1 ( 2 m - l ) ( m ) and Equation 2.17 (the second form of the polynomials) can be transformed into Equation 2.21 by means of _ (-l)mm!(iV - 1)('»> CmN = ; -7--N ( 2 . 2 4 ) 2(2m - l ) ( m ) with its respective initial conditions. A recursive form for the creation of Gram polynomials is given in [99] for N points, k, in the interval -(N - l ) /2 to (N - l ) /2 as3: 1 -Pm +i(*) = - P m W - — P m - l W £ m + l £m £m—1 m {N2-m2) Pm = -777-2 7\ (2 '25) 4(4m/ — 1) _ (2m)! 1 &rr> - — (m!)2 (N - l ) ( m ) This can be transformed into the form [54] P / M - Q V - W 2 ™ - 1 ) 2k P m (m-l)(N-l + m)n Fm{k) = T— r — -P m-i(k) T—: r P m-2{k) (2.26) m(N — m) N — 1 m(N — m) for m > 1 , with Po(k) = 1 and P-\(k) = 0. The nonrecursive form is given in [66, 99] as: - ( _iy+- ( j + mp)(iYfa + A .)0) The first few polynomials, calculated with the recursive or the nonrecursive equations, are: A r - 1 P^-1-6-^r+6l 2 ( i v - i ) ( i v - 2 ) (2-2*) ^ + k (HL=l+k)(^ + k-l) ( ^ + * ) ( 3 ) 3 It is sometimes easier to assume that N — 1 is an even number, such that (Ar — l ) /2 will be an integer, but it is not required in any of the analysis that follows. 18 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials Using the same index as usual, i.e. let i run from 0 to N — 1, the recursive Equation 2.26 can be written as m(jv — m) V N — 1J m(N — m) which gives the same polynomials as the general form (Equation 2.22) with CmN = 1. A number of properties of DLOP's are listed in [93] and it also lists an algorithm for constructing the polynomials with a minimum number of operations. The algorithm is as follows: Initialization: n = N-l co = 0 6i=c\+c\ c\ = n - i - i 6 2=n C2 = n P 0 = 1 (2.30) n Iterations for m = 2,3, Pi = ci/n (2.31) 6$ = 6Q + 2 62 = 62 — 2 CQ = Co + 60 C\ — C\ + 61 C2 = C2 + 62 Pm = ( c i P m - 1 - CoP m-2)/c2 Since m, N, and i are integers, the coefficients CQ,CI,C2 and <5o,^ 1,<?2 will also be integers. The total number of operations required for this scheme to calculate a single point, i, is: m real divisions 2m — 2 real multiplications m — 1 real subtractions 4m — 3 integer additions m + 1 integer subtractions 19 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials The other Ar — 1 points would need the same number of real operations, however if the integer additions and subtractions are stored, the number of integer subtractions would only be two for each point and the number of integer additions would be m. This scheme does not require any integer multiplications like Equation 2.29, and does not require rn2 operations as in Equation 2.19. There is also an algorithm to compute the multiple of the binomial coefficients in Equation 2.19, [ "' J f mt-? j , by forming a simple matrix using addition only, similar to the formation of Pascal's triangle, and then multiplying two elements of the matrix to compute the coefficients [1]. As a last resort, tables of orthogonal polynomials can be used. Fisher and Yates [49] give tables of integer orthogonal polynomials for equally spaced points for order 1 to 5 and number of data points from 3 to 75. They mention that more extensive tables are available for up to 104 data points and also for order up to 9 but only for 52 data points. With today's computers it is quite easy to generate all necessary integer polynomials, or scaled such that the first or the last value of every polynomial is equal to 1. Care must be taken when dealing with limited word lengths to represent large numbers since round-off errors might become a problem. Unless stated otherwise, in this thesis, Equation 2.29 is used in order to create the Gram polynomials, since it will produce polynomials that are all of comparable magnitude as long as the number of data points is large enough compared to the polynomial order. Figure 2.6 shows a plot of the first 30 polynomials for 100 data points. It can be seen that the higher order polynomials become larger and larger as the polynomial order gets close to the number of discrete data points. Figure 2.7 shows at what order the maximum of the absolute value of the polynomials is greater than 2, 10, and 100 for number of data points between two and 1000 points. It will be shown in chapter 5 how this large absolute value for a polynomial will adversely affect the controller output. 2.1.2 Gram goes to Legendre It has already been mentioned that a certain form of the Gram polynomials has exactly the same parameters as the shifted Legendre polynomials. The Gram polynomials with c, n^ — 1 will also converge to the continuous Legendre polynomials for infinite number of data points [66]. It is 20 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials 40-20-0--20--40: 0 20 40 60 80 10 15 Data points 100 0 Order 20 25 30 Figure 2.6: Plot of the first 30 Gram polynomials for 100 data points. observed that the Gram polynomials (Equation (2.22)) can be written as: Pi(i) =x(i) = 1 - 2 N-l i f ..,,3 (N-l) 2 „„ ,_ . „ N 2-l p^^-^p^y-^-^-ww (N-2) (2) 1 P4(i) = -i35(P 1(t)) 4 (N-l) 30(Pi(»)) 2(N^1)(N2 - f ) ( o(N + 3)(N + l) (N - 4)(N - 2) + 3 (2.32) v (N-2)(3) (W-2) ( 3 ) In the limit when N goes to infinity the fractions involving N in the above equations go to one and the above equations can be written as: Pi(x) = x P2(x) = -{3x2 - 1} P3(x) = ±{5x3 - 3x} P4(x) = - { 3 5 a ; 4 - 3 0 x 2 + 3 } 8 (2.33) 21 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials 1201 1 1 1 1 1 1 1 1 r 0 100 200 300 400 500 600 700 800 900 1000 Number of data points Figure 2.7: Relationship between the number of data points and the order at which the maximum absolute value of the polynomial is greater than 2, 10 and 100. The circles in the lower left hand corner denote a fourth and an eight order polynomial for 20 data points. These equations are exactly the same as the Legendre equations on the interval [—1,1] [11]. Figure 2.8 shows the fourth and eighth order Gram polynomial for 20, 30, 50, and 100 points as well as the continuous Legendre polynomial. It can be seen that, as the number of data points is increased, the Gram polynomial approaches the continuous Legendre polynomial. Furthermore, the 20 point Gram polynomial is much closer to the fourth order Legendre polynomial than it is to the eighth order Legendre polynomial, because fourth order is small in comparison to 20 points whereas eighth order is not, as can be see from the circles in the lower left hand corner in Figure 2.7. The requirement for Equation 2.32 to be written as Equation 2.33, i.e. for Gram to go to Legendre, is then N 2> m. 2.2 Approximation s b y Gra m polynomial s By approximating the CD profile with polynomials, the shape of the profile can be controlled, and the number of parameters to control can be greatly reduced. This amounts to controlling the low frequency components of the profile. This idea has been used in Sendzimir steel rolling mill during research done at the University of Strathclyde [101, 102] and the University of Sheffield [41, 56]. Now define a matrix, xn, consisting of the first m orthogonal polynomials (it is assumed that the mean value of the profile to estimate is zero, i.e. PQ does not need to be included) evaluated 22 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials 0.5 0 0.5 i i i i i « \ ^s/" 1 1 ! Fourth order polynomial i i i ^s*^ 20 points ^ W — — 30 points • — • 50 points 100 points , legendre ( I I I / A -l 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Oi Eigth order polynomial -0.8 -0.6 -0.4 Figure 2.8: Gram goes to Legendre. The smooth solid line is the Legendre polynomial in both plots. at N points as (see Equation 2.2): Xo = [Pi P 2 ••• P„ (2.34) The relationship between the profile, y, and the parameters of the orthogonal polynomials, y c, is given by y(t) = x$y c(t) + e(t) (2.35) where e consists of terms of order m 4- 1 and higher (more specifically, order m + 1 to N — 1, because a polynomial of order N — 1 will exactly match the measurements at the N data points). The least squares estimate for the y c vector is: T\-K yc(t) = (XOXQ ) xoy(t) (2.36) The matrix XQXQ is a diagonal matrix because it consists of orthogonal vectors, so its inverse is easily calculated. In the case of orthonormal polynomials, the matrix simply becomes the identity 23 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials matrix. The least squares estimate (Equation 2.36) can be written as: iic(t) = p f P i o 0 0 pfpl 0 0 1 P I P n (2.37) " pfpTP^ " pfpTP^ 1 'pTy P T P m m " It is clear from the above equation that the elements of the yc vector can be calculated independently of each other. If more parameters are needed the next parameter can be calculated from: 1 yc(t,m + l) = vT p Pm+iy (2.38) where the second variable in yc(t, m + 1) denotes the m + 1 element in the parameter vector. Hence, more parameters can be added without affecting the ones already calculated. Figure 2.9: Approximation of an industrial basis weight CD profile with Gram polynomials to the 10th order. The small inserted graph shows the 10 parameters of the Gram polynomial approximation. 24 Clmpter 2: Discrete orthogonal polynomial approximation using Cram polynomials 1.5 < 60 r o 60 '5 -0.5 -1 -1.5 -2 b Measured data 4th order approximation 10th order approximation 20th order approximation 40th order approximation 50 100 150 200 250 CD 300 350 400 450 Figure 2.10: Measured data and 4, 10, 20 and 40th order approximation using Gram polynomials. An example of approximating an industrial basis weight CD profile, with the polynomials defined by Equation 2.21 and Cm^ = 1, is shown in Figure 2.9. The measured profile, y, is drawn with a solid line connecting the measured points. The inset bar graph in the upper left corner shows the values of the estimated parameters, y c. The estimated profile, y, is drawn in the big graph with a dashed line. Here, the order of polynomial approximation is chosen arbitrarily, but there exist methods that could be employed to estimate the optimal polynomial order for a given profile. Some of these methods will be discussed below. The particular profile in Figure 2.9 has been approximated with different polynomial orders and the results are shown in Figure 2.10. The parameters needed for approximation are shown in Figure 2.11. From this plot, it can be seen that little improvement in approximation is obtained in using polynomials from the 22nd order to 40th order because most of these parameters are small in comparison with those for the 21st and 22nd order. Furthermore, the plot of 2a, twice the standard deviation of the approximation error, versus order, shows that the largest drops are indeed where the parameters are the biggest. The figure shows that the standard deviation gets smaller and smaller with increasing polynomial order, and it will in fact be zero when the order of the polynomials is equal to one less than the number of data points. In practice, the 25 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials Parameters of the Gram polynomials 0.2 0 -0.2 -0.4 I - l n _L -i U i - i i i i i i J] fi n (I n fi n „ ^ -u u i i i i uu L u u ^ u u -1 1 1 1 0 10 15 20 Order 25 30 35 40 0.7 e-eee o o o e °^H> 10 15 20 Order 25 30 ooo(t>f»eeoo 35 40 Figure 2.11: Parameters of the Gram polynomials in Figure 2.10 and two times the standard deviation of the approximation error. measured data has noise and it is not desirable to approximate noise. In the words of Hildebrand [66]: "it is foolish to attempt to determine a polynomial of high degree which fits the vagaries of such data exactly and hence, in all probability, is represented by a curve which oscillates violently about the curve which represents the true function. " The measured profile, together with the approximation error for a 40th order approximation, is shown in Figure 2.12. It can be seen that the low frequency components of the measured profile have disappeared and the peak to peak value is lower. Therefore, if the lower order polynomials are removed from the profile, the profile becomes more uniform, and the problem with different averages for different sections of the profile — a major concern in CD control — is eliminated. 2.2.1 Autocorrelatio n Assuming the true profile is an orthogonal polynomial and the measured profile is the true profile mixed with white noise, the "whiteness" of the residues could be used to calculate the optimum order of polynomials to use. If the residues are white the autocorrelation function would have a value of one at zero lag and a value of zero for all other lags. For finite length sequences the value of the 26 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials 150 20 0 350 400 450 y-ye(40) 1 0 -1 -2 50 100 150 200 250 CD 300 350 400 450 Figure 2.12: Measured data and the measured data with 40th order approximation subtracted. autocorrelation sequence for non zero lag will not be exactly zero. Therefore, a confidence interval or error bound is calculated which the autocorrelation should fall within. The confidence interval depends on the sample size used. A widely used measure for autocorrelation functions is to use 95% confidence interval, for number of data points greater than 60 it is 2/y/N, where N is the number of data points. If the autocorrelation function is white, it can be said that one out of 20 points in the autocorrelation function can lie outside the 95% confidence interval by pure chance. A part of the autocorrelation function of the four approximations from Figure 2.10 is shown in Figure 2.1?. The percentage of the autocorrelation of the residues for the 4th, 10th, 20th, and 40th order approximation, that is inside the 95% confidence bounds, are 88%, 91%, 90%, and 89% respectively. None of the approximation errors can be considered white and the whiteness does not seems to be related to the number of polynomials used. The whiteness of the residues is therefore not a good measure of fit unless something is known about the noise. 2.2.2 Akaike' s Informatio n Criterio n A perfect fit can be obtained if a polynomial order is chosen as one less than the number of data points. Such a model would be overly complex and would approximate the noise. A balance should 27 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials J I 1 I I 1 I L 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 10 0 Lags Figure 2.13: Autocorrelation of a profile minus the Gram polynomial approximation for 4, 10, 20 and 40th order approximation. The 95% confidence interval is shown with the dotted lines. therefore be found between the complexity (parsimony) of the model and the quality of the fit. A model should be as simple as possible, yet sufficiently flexible to cover all aspects of the process. The most natural way to search for a suitable model structure is simply to test a number of different ones and compute the resulting models. One formalized way has been pursued by Akaike (see for example [3]). For model errors that are Gaussian with unknown variance, the Akaike's information theoretic criterion (AIC) is written as AIC(m) = log (1 £ (y - y{m)f J + ^ (2.39) where N is the number of data points, m is the number of parameters and y{m) is the profile approximated with an rath order polynomial. The first term in the above equation measures approximation error and will decrease with increasing model complexity; the second term accounts for model complexity and increases sequentially. The idea is that as the approximation error becomes smaller with increasing complexity, a minimum is reached where reduction in prediction error is more than offset with the increased model complexity. Figure 2.14 (the solid line) shows a plot 28 Chapter 2: Discrete ortliogonal polynomial approximation using Gram polynomials T 1 1 1 1 1 1 1 r 2*Order 4*Order 6*Order 0 10 20 30 40 50 60 70 80 90 100 Order Figure 2.14: Akaike's Information Criterion for the profile from Figure 2.9. The solid line shows the AIC, the other two lines show the AIC when it is modified to greater penalize model complexity. The circles in the plot denote the minimum for the curves. of the AIC for the profile from Figure 2.9. It is found that the optimum polynomial order for this particular profile is 86. If the order is increased further the increased complexity does not warrant the reduction in the approximation error. A deficiency in the AIC is that it only penalizes the addition of parameters in calculating the approximation error. For a general ARMA type model an increase in the order of the A and/or the B polynomial is usually required when increasing model complexity. This case only requires the input and output to be stored for one extra time step. On the other hand, when a polynomial order is added to a model using orthogonal polynomials, it is also necessary to calculate the new polynomial using a recursive equation. If we look at the complexity of a model as a cost associated with computing, it is necessary to factor this complexity into the AIC. If it is assumed that the cost of calculating the new polynomial is the same as the cost of increasing model parameter vector, the dashed line in Figure 2.14 is obtained. This criterion is used in [83] to reduce the probability of selecting a model with one more parameter than the true parameter vector. It can be seen that the optimum polynomial order is now 22. If the cost of creating an extra polynomial is assumed to be double that of increasing the parameter vector, the dotted line is Figure 2.14 is obtained, which has a minimum at polynomial order of 10. To guard against approximating the noise in the profile, a polynomial of high degree should -1.1 -- 1 . 2 -- 1 . 3 -29 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials not be used [66]. It will be shown in the next chapter how the condition number of a transfer matrix worsens when the order of the polynomial approximation is increased. Clearly the AIC in its original form gives too many parameters. The form suggested in [83] gives a good balance between complexity and fit. 2.3 Frequenc y characteristic s o f Gra m polynomial s The frequency characteristics of the Gram polynomials were investigated using Fourier analysis. It was found that, unlike the amplitude of the polynomials, the frequency characteristics do not depend on the number of data points used because the number of zero-crossings remains the same. For this analysis, a 100th order polynomial is investigated and 212 = 4096 data points were chosen because it is convenient to calculate fast Fourier transform (FFT) for sequences that have numbers of elements that are powers of two. If the zeros of the polynomial are equidistant, the polynomial would be expected to have its major frequency around (100-l)/2=49.5 cycles per profile. Figure 2.15 shows the left half of the 100th order Gram polynomial defined for 4096 points (the right half is identical because even order polynomials are symmetric). The left half is then further divided into 4 quarters and a frequency distribution for each part was calculated. Figure 2.16 shows the spectrum of the 100th order Gram polynomial defined over 4096 points in the interval [—1,1]. A 100th order polynomial has 1°2~1 = 49.5 periods of oscillations. It is seen in Figure 2.16 that the polynomial has most of its frequencies distributed slightly below 49.5 periods over its entire profile. The reason that it is lower than 49.5 periods is that the "edges" of the polynomials have slightly higher frequencies, i.e. more cycles per sheet width. This is evident in the leftmost quarter in Figure 2.15. Therefore a large portion of the 49.5 cycles are used for the first quarter. The remaining three quarters will have fewer cycles as can be seen in Figure 2.17, where the spectra is calculated for the four different 512 points sections of the left half of the profile. The edge quarter from Figure 2.15 is shown to have frequency component corresponding to about 65 periods over the sheet width, whereas the other three quarters have frequencies corresponding to between 30 and 40 cycles per sheet width. To show how the Gram polynomials can be used to remove some frequencies, a 1000-point profile is created which consists of 5 sinusoidal waves with different frequencies (1,4, 10, 30, and 100 cycles), added together. The dotted line in Figure 2.18 shows the profile created. The solid line shows the profile when the first 20 Gram polynomials have been removed from the profile. The bottom figure shows the frequencies in the original profile (the dotted line) and the profile with the 30 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials Figure 2.15 : Left half of 100th order Gram polynomial for 4096 points. It is shown with the polynomial partitioned into 8 zones, i.e. the left half is partitioned into 4 zones, separated with different line types in the graph. 101 10" -. 10" 10" 10" 1(T 10" : : : : : : : : I : ; : : : : : I : : : : : : : I ; : : : : : : a : : : : : : : a : : : : : : : : i : : : ; : : : : l : : : : : r : : i : : : : : : ; : i : : : : : : : s 0 20 40 60 80 100 120 140 160 180 200 Cycles per sheet width Figure 2.16: Spectrum of 100th order polynomial. 31 Chapter 2: Discrete orthogonal polynomial approximation using Gram polynomials 1 0 t : , : : : : : ! : : : , ; ' : : ! : : : : : : : ^ : : , : : : J : : : : : :: : l : , : : : :: : l ; : : : : : : : l Cycles per sheet width Figure 2.17 : Spectra of 100th order poly, calculated for the four sections of the left half of the polynomial shown in Figure 2.15. 20th order removed (the solid line). As expected, the low frequencies, i.e. few cycles per sheet width, are removed, and the 10 cycles per sheet width component is slightly lowered. There are a few frequency components added around 10 cycles, but their power is low compared to the original spectrum. The 10 period component is fully eliminated if the profile is approximated with the first 40 Gram polynomials as shown in Figure 2.19. The only frequencies left are 30 cycles and 100 cycles. To remove or reduce these remaining frequencies, at least 61st and 201st order polynomials, respectively, would be required. 2.4 Summar y In this chapter discrete orthogonal polynomials with equal weight given to every data point, Gram polynomials, have been described. It is shown how a set of data, a profile for example, can be approximated with Gram polynomials. Several methods are given to create the Gram polynomials, using both sequentially and recursive algorithms. A unified sequential formula is given for all the methods. It is observed that the Gram polynomials will approximate the Legendre polynomials when the order of the polynomials is much smaller than the number of data points which gives the Gram polynomials another name, discrete Legendre orthogonal polynomials. Some methods are discussed 32 Clmpter 2: Discrete orthogonal polynomial approximation using Gram polynomials 20th order approximation 100 200 300 400 500 600 700 Databoxes across the sheet 800 900 1000 60 80 Cycles per sheet width 120 Figure 2.18 : The original profile (dotted line) and the residual profile (solid line) when a 20th order Gram polynomial is removed. Power spectrum for the original profile (dotted line) and the profile when 20th order Gram polynomial is removed (solid line). 40th order approximation 100 200 300 400 500 600 700 Databoxes across the sheet 800 900 1000 60 80 Cycles per sheet width 100 120 Figure 2.19 : The original profile (dotted line) and the residual profile (solid line) when a 40th order Gram polynomial is removed. Power spectrum for the original profile (dotted line) and the profile when 40th order Gram polynomial is removed (solid line). 33 Cliapter 2: Discrete orthogonal polynomial approximation using Gram polynomials for finding the optimum number of polynomials to approximate a profile. It is found that a modified AIC gives a good balance between complexity and fit. Finally the frequency characteristics of the Gram polynomials are investigated and they are found to have their main frequency component slightly lower than half their order. 34 Chapter 3 Model fo r cross-directiona l basi s weigh t respons e A basis weight control system is made up of three main parts: i) the sensor, a gauge which scans the paper sheet, ii) the process control computer, and iii) the slice lip actuators [106]. The process control computer uses the information obtained by the scanning sensors to calculate the new slice lip actuator settings. 3.1 Pape r machin e A paper machine, on which stock is drained and dried, is used to manufacture a sheet of paper of infinite length and up to 10 m wide. Thick stock, a water/fibre mixture with a fibre concentration around 2%, comes from the machine chest. The thick stock is diluted with white stock and sent to the headbox where the concentration is about 0.5%. The purpose of the headbox is to take the inlet stock flow and uniformly distribute the proper amount of fiber suspension across the paper machine moving forming wire at the right speed and in a uniform way, through an opening which is between 50 — 100 mm [98]. The product coming out of the paper machine will almost reverse the ratio of water to fibre, since typical newsprint paper contains 5% moisture. There are two ways of forming the sheet on the moving wire, either by a single wire (fourdrinier) or by a twin-wire former, which has wires on both sides of the forming sheet. The sketch of a simplified paper machine in Figure 3.20 shows a Fourdrinier machine. On the wire the fibres form into a sheet after water drainage using gravity and suction, and in the case of twin wire, also by pressing. On the fourdrinier wire there is a visible line of demarcation corresponding to the point where a glossy layer of water on the top of the stock suddenly fades away [112]; this line is called the dry line. The sections of the profile with high basis weight will require a longer time to dry than parts with low basis weight. Therefore, if the dry line is measured it could be used to predict the basis weight profile. In the dryer section the fibre bonds develop and more of the remaining water is evaporated by heating. In the calender section the sheet is pressed between metal rolls to reduce thickness and smooth the surface. Finally the paper is wound onto a reel [111]. Moisture, defined as percent water by weight, and basis weight (BW), defined as weight per unit area, are two of the most important characteristics of paper. BW includes both the weight of the 35 Chapter 3: Model for cross-directional basis weight response Stock valve =tfc= Calender Figure 3.20: Paper Machine. Figure 3.21: O-frame. water and of the fibres, therefore, in order to measure the amount of fibres in the paper the weight of the water has to be subtracted to obtain the dry weight. Burkhard and Wrist [16] proposed, a now widely accepted, classification of basis weight variations in three components: 1. The machine direction (MD) component which is time varying and independent of the location across the machine 2. The cross direction (CD) component which describes the profile and changes much slower than the MD. 3. The remaining signal is the random (R) component and varies in both time and CD position. Lately, some researchers have also started looking at variations in the vertical or thickness direction (ZD) [29]. The variations in ZD are primarily caused by movement of fines, the small fibres, and fillers, used for creating a smooth surface. The paper characteristics are measured with sensors mounted on an O-frame, i.e. a guiding rail on which they travel back and forth across the machine (see Figure 3.21). Due to the MD movement of the paper and the CD movement of the sensor, the measurements form a zigzag pattern on the paper 36 Chapter 3: Model for cross-directional basis weight response sheet as shown in Figure 3.21. The measurements are thus a composite of the CD and MD variations, from which the CD and MD profiles must be extracted. A simple way is through an exponential filter. Let x n(t) denote the raw measurement obtained at the measurement box n during scan t, then yn(t) = (1 - a)y n(t - 1) + ax n(t) (3.40) where y n(t) is the estimated profile for box n and scan t, and 0 < a < 1 is the exponential filter weighting factor [33, 40]. The MD profile can be calculated as the average of y n{t) over one scan and the CD profile would then be the deviation of yn(t) from the average. This method of separating the MD and CD is very slow because a small a is used to separate the MD and CD. Consequently, control action is taken every 3-5 scans so that the result of the last control action is clearly seen before a decision is made to do another control. Over the last few years, methods have been proposed in order to estimate the moisture CD profile [20, 23, 39, 40, 70, 92] and the basis weight CD profile [122] more quickly. Vendors have increased the speed of the scanner and sample the profile in smaller and smaller data boxes in their quest to reduce the variability of the paper and to better separate the CD and MD [115]. Using analysis of variance, some have come up with a method to alternatively decrease or increase the speed of the scanner in order to get the best measurements [2]. Others have used analysis of variance to quantify what can and cannot be done to reduce the variations [114]. Recently, MD variations have been estimated within a scan, and hence faster MD variations are reduced, resulting in overall reduction of MD variations [118, 120]. The algorithm in [70] consists of a modified least-squares parameter identifier for estimating cross direction profile deviations, a Kalman filter for estimating machine direction disturbances, and takes into account the nonlinear relationship between the MD and CD observed by Lindeborg [84, 85]. It is observed for basis weight [122] that there is no nonlinear relationship between the MD and CD, and that the model for MD variations is a second order model instead of a first order model, as is the case for the moisture model. Basis weight on a paper machine is routinely controlled in the cross direction by bending the slice lip of the headbox. Because of coupling, movement of a single actuator will cause the slice to be deflected between several actuators, even though the actuators settings stay the same for the neighbouring actuators (see Figure 3.22). Hence, moving one actuator of the slice lip will in general affect more than one measurement point in the cross direction. Another source of coupling is waves on the fourdrinier caused by the slice lip opening. Disturbances can also be caused by the wake effect within the headbox [61, 62]. Therefore, two steps are needed to model the response of an actuator movement on the basis weight. The first stages involves modelling the slice lip and its deflection to 37 Chapter 3: Model for cross-directional basis weight response a movement of a single actuator. The slice lip deflection is then calculated to the same resolution as the basis weight is measured. Slice lip deflection will cause a wave to form on the forming fabric and this wave will travel in the CD direction and spread further (Figure 3.22). At present, no models exist to describe that wave and consequently, identification of the actuator response to the measured basis weight is needed. Therefore, the paper machine is usually represented by a block diagonal matrix with number of rows equal to the number of basis weight measurements, and number of columns equal to the number of slice actuators. The block diagonal has to be identified. Usually, several responses across the machine are identified. Then it is assumed that the response to a slice lip movement is the average of those several responses, and is therefore the same, regardless of where it is on the machine (with the possible exceptions of the edges). I t will be shown in this chapter that it is not possible to have the same response unless the number of measurements is an integer multiple of the number of actuators. The interaction matrix can have eigenvalues close to zero; therefore, its inverse, which could be used in decoupling the actuators, is not easily calculated. 1 0.8 0.6 0.4 0.2 -0.2 I 4 + + + * * + o + Actuators Slice lip Basis weight + + + + + + a>++++®+4-++®+4-+|®a- + ® ® ++®++4-+®++-*-+®++++ ^ _1 1 I L -4 0 1 CD actuators Figure 3.22: Response to a slice lip movement at one actuator position. The slice lip deflection is denoted by • with 0 denoting the deflection at a slice lip actuator. The + marks what the basis weight might look like. 3.2 Basi s Weigh t It is claimed in [96] that the dominant constituent of the response is the result of the water waves created on the fourdrinier wire originating from the slice lip deflection, and setting somewhere before the dry line. As said before, no models currently exist in modelling the waves that appear on the 38 Chapter 3: Model for cross-directional basis weight response wire in the paper machine. For that reason the response model from the slice lip deflection to the basis weight has to be identified from the input-output data. It was noted early on by researchers looking into CD control problems that the process response is not the same in the center of the paper machine as it is at the quarter points or at the edges [106]. The main reason is the asymmetric deflection of the slice lip at the edges. Part of the edge phenomena also relates to waves being bounced back from the deckle board, the board that is used to contain the stock on the wire. A method to reduce or eliminate the deckle waves is reported in [98]. One distinctive characteristic of all paper machine cross directional models is their large dimen-sion. To measure basis weight, a beta ray gauge traverses the paper, measuring at N points across the sheet. Typically N can be between 50 and 500 although recent systems can provide as much as 2000-3000 points. Those measurements are then used to calculate settings for up to 100 actuators located upstream at the slice lip of the headbox (see Figure 3.20). The calculation of the actuators settings is based on CD response models, which have three common elements: actuators dynamics, time-delay, and interactions. It is commonly assumed that all actuators can be modeled by the same first-order lag dynamics 1 - P a 9 _ 1 (3.41) where ka is the gain and pa is the time constant. The sampling method currently employed on paper machines uses a sampling period that is much longer than this time constant. Hence, the actuator dynamics can be ignored in most cases, simply because it is not possible to observe it. Since the basis weight measurements are taken at the finishing stage of the paper machine, and the slice lip actuators are located on the head box, some time will pass before changes made at the slice will affect the measurements. Given that the transportation lag is written as q~ d, the transfer function from an actuator to its corresponding measurement would be: 9a{q'l) = , ^ ^q- d (3.42) ; 1 -paq ] Control action by one actuator generates water waves on the wire which influences the basis-weight at several locations around the actuator. The number of locations affected depends on the type of paper being manufactured — heavier paper produces wider interactions [72]. Assuming the same 39 Cltapter 3: Model for cross-directional basis weight response spatial response for all actuators, the response matrix, G, will therefore be banded diagonal: "Po Pi ••• Pk 0 • • • 0 " Pi Po Pi Pi '•• Pk '•• '•• 0 '• • '• • G 0 0 Pk Pi Pi Po Pi 0 Pk pi PO (3.43) Here, it is also assumed that the response is symmetrical such that the elements above and below the diagonal are the same, and the width of the CD response is 2k + 1 measurement data boxes. It is desirable to emphasize the fact that the matrix need not be symmetrical. The elements below the diagonal would then be subscripted with a minus. The input output relationship is: y(t) = G fca l-paq-iq -d u(t) (3.44) Moving an actuator flexes the slice lip and puts waves on the wire. The slice lip deflection can be modeled and will contribute to part of the structure of the interaction matrix. The deflection can be measured directly on the slice or observed in the measurements if there are more measurements than the number of actuators. The measurements will, however, contain more than just the deflection of the slice. They will also contain the waves produced on the wire by slice lip deflection. Work is under way at the Pulp and Paper Centre to model the wave process. It should be emphasized that the interaction matrix does not only describe the slice lip deflection but also the waveform that travels in the cross direction as a result of the slice lip deflection. Several interaction matrices have been published or can be read from low resolution graphs in the literature [14, 24, 21, 28, 64, 72, 100, 116, 119, 126, 127]. Thirteen such matrices are reported in [79] for a wide range of paper machines, from light newsprint to heavy cardboard. Po Pi Pi P3 Pi P5 P6 P7 PS P9 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.4 0.2 0.5 0.1 0.45 0.4 - 0 . 1 5 0.2 0.4 1.2 1.3 0.9 - 0 . 5 - 0 . 3 -0 .55 - 0 . 5 0.03 - 0 . 1 - 0 . 2 0.6 0.8 0.7 0.05 - 0 . 0 1 - 0 . 1 - 0 . 4 - 0 . 4 - 0 . 6 0.8 - 0 . 2 - 0 . 9 - 0 . 3 1.0 - 0 . 2 0 0.6 - 0 . 2 - 0 . 1 - 0 . 5 (3.45) -0 .4 -0 . 2 -0.2 40 Chapter 3: Model for cross-directional basis weight response -5 0 5 Relative CD location 10 -5 0 5 Relative CD location Figure 3.23: CD responses, according to Laughlin's interpretation from the literature [79], see also Equation 3.45. The corresponding responses to a slice lip actuator movement are shown in Figure 3.23 for the first nine responses (left figure), and the last three responses (right figure). It can be seen that the lighter grades have a negative side lobe in most cases but the response is quite narrow compared to the heavy weight paper where the center of response is not necessarily the largest response [21]. One of the major problems with banded diagonal matrices is that they can easily approach singularity — their condition number can grow indefinitely. The condition number of a matrix is the ratio between its largest and smallest singular value. Singular values are calculated for real valued matrices using: ^g(GTG) (3.46) where eig(-) is the eigenvalue function. For example, the determinant of a 5 x 5 matrix with po = 1, Pi = 2r, and j>2 — —r, where r is a parameter used to vary the magnitude of the side lobe, is: 1 2r —r 0 2r 1 2r —r —r 2r 1 2r 0 —r 2r 1 0 0 —r 2r det 0 0 - r 2 r 1 56r5 + 42r4 - 24r3 - 19r2 + 1 = (3.47) 56(r - 0.6404)(r + 0.7985)(r + 0.4163)(r + 0.3909)(r - 0.2149) 41 Chapter 3: Model for cross-directional basis weight response If r is equal to any of the roots of the determinant then this matrix will be singular and its condition number is infinite. Reciprocal condition numbers (condition numbers between zero and one, where zero indicates singularity) can be plotted for a matrix with po = 1, p\ = 2r, and pi = —r where r goes from 0 to 0.5 for different sizes of the G matrix (see Figure 3.24). It can be seen when there is any significant cross coupling it is very likely that the matrix will be ill conditioned, i.e. reciprocal condition number close to zero. For a similar plot of condition numbers of G matrix with Po,i,2 = [1, r, —r] see [79]. Processes with high condition number can be difficult to control, especially in the presence of uncertainty in the interaction matrix. Because of the uncertainty, the direction of a large input may not correspond exactly to the low plant-gain direction, and the amplification of these large input signals may be much larger than expected from the model. It can be shown that disturbances entering the process in the direction of the minimum singular vector for a process with bad condition number can lead to a poor performance, because the control action can be taken in the wrong direction [107]. In the limit as the G matrix becomes singular, the plant is uncontrollable, as the actuator settings cannot be determined based on the measurements [79]. Singularity is not restricted to the two forms of the interaction previously mentioned, it is also of concern in matrices of the form po.i = [1J r], Po,i = [1, —r], and po,i,2 = [1, 2r, r]. In later chapters, a paper machine with 105 actuators will be employed, it is of interest therefore to look at a 105 x 105 matrix with po = 1, pi = 2r, and p2 = —r. Its reciprocal condition number for r between 0.15 and 1.0 is plotted in the upper half of Figure 3.25. The corresponding CD response for cases where the matrix is singular is plotted in the lower half of Figure 3.25. These responses are typical of those reported in the literature for paper machines. Therefore it can be concluded that singularity of the interaction matrix can quite easily happen on paper machines. It can be seen from industrial data that the model of the slice lip movement to the dry weight measurements is simply a gain and a delay. This is shown in Figure 3.26 where the mean values are removed from the input and the output. The solid line represents the output or the measured basis weight in one CD box, and the dotted line represents the input or the slice lip adjustment for one actuator. The model for this paper machine at its operating conditions can therefore be written as: y = p Qq-d u (3.48) If the interactions of adjacent CD locations are included in the model it becomes: y = Gq~ du (3.49) 42 Chapter 3: Model for cross-directional basis weight response Size of G Figure 3.24 : Reciprocal condition numbers of G, where it is created according to equation 3.43 with po = 1, p\ = 2r and p2 = —r, for various sizes of the matrix. When G is badly conditioned, the reciprocal condition number is close to 0. 1 0 -1 -4 1 1 1 1 1 1 r' i -2 0 CD Figure 3.25: Reciprocal condition number for a 105 x 105 matrix with po = 1, p\ = 2r and p2 = ~ r - Basis weight response for a machine with a singular interaction matrix, reciprocal condition number equal to zero. 43 Chapter 3: Model for cross-directional basis weight response 2 -0 --1 --3 r — — — — — — ~\ r - — — — — — , — I I ' ' I ' I I I ' I ' Output in g/mA2 *-' ' 1/ I — — Input in lCpm ' ' . \ ' ' I ' | 10 20 30 MD 40 50 60 Figure 3.26: Dry weight measurements and the input to the corresponding slice lip location. where G is a block diagonal matrix. The G matrix will be estimated in chapter 4. Now we will look at some of the forms the G matrix will take. 3.2.1 Squar e interactio n matri x When the number of actuators is equal to the number of measurements, the G matrix will be square. There are two main problems associated with this matrix form that are also prevalent for a nonsquare G matrix. The most serious problem is mapping, that is, determining which actuator belongs to which data box. Another problem is how to model the edges. If the edges are not accounted for, the model might show some edge effects that are not present in the real process. As an example, assume that G, defined by Equation 3.43, is a 20 x 20 matrix with k — 3 or. G = P0 Pi P2 P3 0 0 Pi Po Pi P2 P3 0 P2 Pi Po Pi Pi P3 P2 Pi Po Pz u 0 P3 P2 Pi P2 P3 U U Pi Po Pi P2 P3 u P3 P2 Pi Po Pi P2 P3 0 P3 P2 Pi Po Pi P2 0 P3 P2 Pi Po Pi 0 0 P3 P2 Pi Po (3.50) where the edge elements of the interaction matrix that might need adjustments are bold faced. If the actuators are set such that the slice is a straight line from -1 to 1, the resulting basis weight profile, 44 Chapter 3: Model for cross-directional basis weight response i t _ i i i i i i i i I 2 4 6 8 10 12 14 16 18 20 2 4 6 8 10 12 14 16 18 20 CD Figure 3.27: Plot of the weight response to a straight line slice lip profile. The upper half is using G matrix from Equation 3.50 and the lower half is if the edges of the matrix are made smooth. with po = 1, pi = 0.6, pi = 0.2, and pz = —0.2, chosen arbitrarily, will have ripples at the ends (upper half of Figure 3.27). The ripples are artifacts which arise because the sum of the middle rows of the G matrix is equal to po + 2(pi + P2 + Ps) = 2.2 whereas the first three rows have column sum equal to 1.6, 2.2 and 2.4, respectively. If the G matrix is fixed such that the first three elements in the first column (bold faced in Equation 3.50) are set to po +Pi + P2 + P3, Pi + P2 +P3. and p2 +P3 (and the same thing is done for the other edge), the effects of the edges of the matrix are reduced somewhat (lower half of Figure 3.27). However, the ripples are not completely removed from the output profile. To do so, it would require tuning of the G matrix based on the input signal. 3.2.2 Nonsquar e interactio n matri x When the number of actuators does not match the number of measurements the interaction matrix will be nonsquare for the unmapped measurements. Again, the edges of the interaction matrix might have to be identified separately otherwise, the model might create some phantom ripples at the ends of the BW profile. For control purposes it is very important that an accurate mapping of the slice screw positions to the measurements data boxes of the weight profile is established [96]. If the number of data boxes is an odd integer multiple of the number of actuators, and the actuators are 45 Chapter 3: Model for cross-directional basis weight response perfectly aligned to some measurements, the interaction matrix is of the form: P-l P-2 P-3 P-4 P-5 P-6 P-l P3 P2 P\ Po P-l P-2 P-3 P-4 P-5 P-6 P-7 Pi Pa Ps PA P3 P2 Pi Po p-l P-2 P-3 P-4 P-5 P-6 Pi P6 P5 P4 P3 P2 Pi PO P-l P-2 Pi P6 P5 P4 P3 P2 Pi Pa where it is assumed that the number of measurements is four times the number of actuators. Note that in this matrix there are four kinds of rows: [p_4 PO P4 ] [P-5 P-l P3 Pi] [P-6 P-2 P2 Pd] [p_7 p_3 pi p5] If the positive and negative elements are the same, and po for the nonsquare matrix is the same as po for the square matrix, and p4 for the nonsquare matrix is the same as pi for the square matrix, then the outputs of both systems are going to be equal for every fourth point of the nonsquare matrix and every point of the square matrix. The general shape of the output profile should therefore be the same. The nonsquare matrix gives a higher resolution of the profile than the square matrix. Figure 3.28 shows the geometrical relationship between the actuators and the measurements. It is rare (number of measurements are odd multiple of the number of actuators) that an actuator will have a measurement box located at the same CD location. In some cases (number of measurements 46 Chapter 3: Model for cross-directional basis weight response 16 15 14 13 12 11 10 9 I 0 o o 0 o o o o o o o : o o « o : o o o o o o o \ o « o : o o o : o o 6 » : o I I o o o o o o o o o e o o o o o * • » • I o o o o o o o o o O 9 o O 0 o I o o o o o o * o o o o o o o o o o o ® o o <b o o I • o o o o o o o o I Figure 3.28: Geometry when the number of actuators is not the same as the number of measurements. It is assumed that the first measurement is half a data box width from the edge. are even multiple of the number of actuators), the interaction matrix is best written as: P-5 P-6 P-l P-l P-2 P-3 P-4 P-5 P-6 P-l P4 P3 P2 Pi P-l P-2 P-3 P-4 P-h P-6 P-l Pi P6 P5 Pi P3 Pi Pi P-l P-2 P-3 P-4 P-h P-6 Pi P6 P5 P4 P3 P2 Pi P-l P-2 Pi P6 Pi P3 P2 Pi P6 (3.52) In general, when the number of measurements is not an integer multiple of the number of actuators, each column will have different entries and each would have to be identified separately or be found by some weighted mean from an average response. 3.2.3 Interactio n matri x parameterize d wit h Gra m polynomial s It was shown in chapter 2 how a profile can be approximated using Gram polynomials. The actuator settings can be written in terms of Gram polynomials as well. The number of actuators does 47 Chapter 3: Model for cross-directional basis weight response not have to be the same as the number of measurements so another set of orthogonal polynomials will be used. For N a actuators and up to degree m a the actuator settings at time t can be written as: act(<) = XQ auc(t) (3.53) If the above equation and Equation 3.44 are used in the least-squares equation from last chapter (Equation 2.36), the least-squares equation can be written as yc(t) = (xoxjf) xoGg a{s)xoauc(t) (3.54) where the matrix XQ 0 is defined as [Pia P2a • • • Pmaa] — the Gram polynomials with number of data points equal to the number of actuators. For a linear ga(s) an m x ma matrix can be defined as: r \ - i Gx = (XOXQ ) x0G'x0a Using Equation 2.37 the G x matrix, for m a actuators and m measurements, can be written as (3.55) Gx = P'J'PI PJGPla PS'P* m ±2. L PIP m PfGP 2 o P'i'Pl P ^ P 2 . . p'i'Pj PJ;GP2„ P^Pm PTGP P?Pl PTGP P2JP2 P T GP P^Pm (3.56) Having established a form for the G matrix from the previous section, the form or the shape of the G x matrix can be investigated. The form of the G x matrix for a square interaction matrix, G, with dimension 105 x 105 and with the diagonal elements given by ^0,1,2 = [1-0, 0.4, 0.2], i.e. k = 2, is plotted in Figure 3.29 and Figure 3.30. From Figure 3.29 it can be seen that G x is almost a diagonal matrix. Figure 3.30 is a plot of every second diagonal of G x except the main diagonal which is positive and much larger in magnitude than the off diagonals. The odd diagonals which are not shown in the figure are all zero for this particular choice of G. What that signifies is that even order polynomials in the input do not contribute anything to the odd order polynomials in the output, and the odd order polynomials in the inputs do not contribute anything to the even order polynomials in the output. That property does of course depend on the structure of the G matrix. In general it can be expected of a real paper machine to contain some nonlinearities, resulting in odd order polynomials in the input affecting even order polynomials in the output, and vice versa. From equation 3.43 (the matrix does not necessarily have to be a square matrix, if it is not a square matrix only the matrices G,'s will change), it can be seen that the G matrix can be split 48 Chapter 3: Model for cross-directional basis weight response 2-1.5-1-0.5-0--0.5: 0 20 ma Figure 3.29 : Meshed surface of the G x matrix for an interaction matrix of size 105x105 and po = 1, pi = 0.4 and p2 = 0.2. 0 0.1 0.2 1 -i i - * • i . * - • ' ' % • • • ' J K . 1 • • • * • • ' * • ' . I Lower diagonals of Gx i i . . . i .> ; - ;$: y . • •6 .8 10 * " " * ' • * ' ' w 1 1 1 • . * • • • • * -i •m '. * • 1-2 ' *• •I: • . * : .*'' •?'• -iff : ' • * • • '. • * • 14 • • * • . . * • : ' • • * • • / . < " 16 .' • If- • • -x • • • It m 10 12 14 Upper diagonals of Gx 16 18 20 Figure 3.30: Every second diagonals of G x from Figure 3.29. The zeroth (main) diagonal is not shown. 49 Chapter 3: Model for cross-directional basis weight response into several matrices, as G ~p 0I + pi 0 1 1 '•. 0 ' • • 1 1 0 H + Pk 0 ••• 0 1 0 1 '•. LO 1 0 (3.57) or more generally G = p0G0+p1G1 + ---+ PkGk (3.58) where Gk describes the effects moving the slice actuators has on their kth neighbour. The Gx matrix can then be written as Gx = (xox^) xoGx^ = (xox^)_1xo(PoGo +P1G1 + • • •+PkGk)x£a (3.59) = poG Xll +piGXl-\ + pkG Xk where the matrices G Xk are defined as GXk = (xox$)~1x0Gkx$a (3.60) The parameters of the slice lip approximation can then be written in terms of the parameters of the Gram polynomials for the slice as: yc =• G xuc = (PoGr„ + PiGXj + h PkG Xk)uc P0 Pi Ipk (3.61) = {G x„uc G Xluc ••• G Xkuc] The parameters pk in the above equation can be estimated in the following fashion : {A TAylATyc Po Pi Pk-i where GXkuc] (3.62) (3.63) A = [G Xltuc G Xluc • This method of finding the Gx matrix could be used to find initial estimates of the Gx matrix, although its use is limited to the case when the mapping from the actuators to the measurements is well known. 50 Chapter 3: Model for cross-directional basis weight response Also, the structure of G is assumed to be well known and all its diagonal elements identical. Having said that, some results using the equations above to estimate G x will be presented in later chapters. The condition number for the G x matrix can also become quite large. Avoiding G x matrices with large condition numbers can become one of the designing criteria when deciding what order of polynomial to use. The condition number of the G x matrix is generally not a problem as long as the polynomial order is much smaller than the number of data points, or there is small cross coupling between adjacent measurements. Reciprocal condition numbers for a fixed dimension of G x matrix but varying G matrix are shown in Figure 3.31. The G matrix is created as before with the elements on the main diagonals as Po,i,2 = [1> 2r, — r\. If 10 polynomials are used, as in this example, there is no danger of a bad condition number if the G matrix is at least 26 x 26, no matter how severe the cross coupling is between adjacent measurements. If there is very little cross coupling (r < 0.15) there is also no danger of a badly conditioned interaction matrix. Reciprocal condition numbers of Gx (dimension 10x10) 0 . 1 - , : 0 . 0 8 - • • • • : • " " 0.06- • • • • : ' " " 0 .04- - • • • : • " " 0 .02- • • • • : • " " o > ^ ••• •"•''." 30 ^ - ^ _ 25^""""^ Size Figure 3.31: Reciprocal condition number of Gx for a fixed number of Gram polynomials but varying the sizes of the interaction matrix, G, and its off diagonal elements which are p\ = 2r and p2 — —v. Similar results can be obtained by fixing the size of the interaction matrix but varying its elements and the size of the G x matrix as shown in Figure 3.32. The G x matrix will not become badly conditioned as long as its size is less than about one fifth of the size of the G matrix, or if cross coupling is small. The cause of the bad condition number for the G x matrix can be seen in Figure f G 51 Chapter 3: Model for cross-directional basis weight response 3.33, which shows a plot of the singular values for the G x matrix for three different dimensions of Gx. It can be seen that the lowest singular value approaches zero for a bigger G x matrix, and as a result the condition number grows. Reciprocal condition numbers of Gx (G is 105x105) 0 . 1 - . • • • • : • " " 0.08 - • • • • • • ' " " 0 .06 - • • • • : ' • ' " 0 .04 - . - • • : • • ' " 0 .02 - • - • • : • " " oJ^-• •'•"•'•'." 10 ^ ~ ^ \ ^ 31 Size of Figure 3.32: Reciprocal condition numbers of G x for G of size 105 x 105 with Po = 1. Pi = 2r, p2 = ~~r and varying the number of Gram polynomials. It will be shown in chapter 5 how the G x matrix can be used for control. It involves finding an ma x m matrix, Ko, such that G TKQ = I. The rank of G x is the minimum of m and m a, therefore, the rank of G XG^ is also the minimum of m and m a. The inverse of G x, which can be written as K0 = GT(GXGT)~1 (3.64) will only exist if ma > m [94], i.e. the inverse only exists if the m x m matrix GXGX has a full rank. The inverse simplifies to KQ — G j 1 if G x is a square matrix (m — ma). Note that the dimension of G x is rn x ma and will be much smaller than the dimension of G (NxN a) so calculating Ko requires less computations than the plant inverse. 3.3 Slic e lip mode l The actuators are assumed to only apply force but no bending moment. Hence, the moment can be completely specified by the deflection of the slice lip at the points where the actuators are 52 Chapter 3: Model for cross-directional basis weight response x O 3 3 > 3 „ w> 2 X o + Size of Gx 10x10 20x20 30x30 X x O ° 4 - + ' .q>. fy.+ .+.+• + •+•+•+ + +• o o--+-t-10 15 20 Number of singular values 25 30 Figure 3.33: Plot of singular values of G x for an interaction matrix of dimension 105 x 105 with po = 1. Pi = 2r, p2 = —f and with r — 0.3. located. The intermediate values can be obtained by a straight line fit between the nearest neighbors because the moment is defined as a force at a particular point times the distance from the force. This is characteristic of beams with point loading. The slice lip is assumed to have free ends, which implies zero moments and zero shear force at the ends. Furthermore, the slice lip is assumed to have a constant flexural rigidity, meaning that the bending modulus, EI, is constant for the whole slice. Assuming small deflection theory [27], i.e. no horizontal displacement at the ends, the bending moments can be written in terms of the second derivative of the slice lip deflection as M-EI (fu(x dx2 (3.65) where EI can be treated as a constant. Only about 1 % error is introduced if the deflection is up to one twentieth of the length of the beam. For a paper machine the deflection of the slice lip is less than 20 mm for a slice lip of 10 m so that the horizontal displacement can be neglected. 53 Chapter 3: Model for cross-directional basis weight response i I I 1 I I I I I Figure 3.34: Model of slice lip. The moments for a beam with forces acting on the ends (p(0) and p(L)) and at a distance a from one end (p(a)) is (see Figure 3.34) d2u El—^r = M(x) = p(0)x+p(a)(x - a) - (a)b ( 3 6 6 ) = ^ x + P(a)(x ~ a ) Since the bending modulus EI is constant along the beam, integration of Equation 3.66 gives Now assuming the forces acting on the beam are springs with a certain spring constant, the force created by an actuator located at ia is p(ia) = k s(z(ia) — u{ia)) (3.68) where u(ia) is the deflection of the slice at the point ia on the beam, and z(ia) is the actuator setpoint at the same point, i.e. the positive difference z(ia) — u(ia) is the expansion of a spring with a spring constant k s. The above equation can also be written as: u(ia) = z{ia) - ^~^- (3.69) For the beam in Figure 3.34 the force acting at x = 0 is —p(a)b/L so from Equation 3.67 and the above equation „(0) = c2 = 2(0) + ^ (3.70) The C2 coefficient from Equation 3.67 can therefore be found if the force at point a, the length of the slice lip, the spring coefficient of the actuator, and the actuator setting at point 0, are all known. Evaluating the other boundary condition, the c\ coefficient can be found as: 54 Chapter 3: Model for cross-directional basis weight response leading to z(L)-z(0) . , Cl = v ' + p ( a ) Equation 3.67 can now be written as ^r^V^" (3.72) t t (x ) = y t ( £ ' -^ -^ + i ( , _.)»}+» 6£ iX . . z(L) - z(0) + z(0) + _ ; z a — b b+ ——x p(a) (3.73) If the above equation is evaluated at x = a and u(a) is written as u(a) = 2(a) — p(a)/ks the force at point a can be written as: z(a) - 2(0) - *(L)-'(°>a Therefore the slice lip deflection can be written as: (3.74) u(x) =• ^{t^-^-^x+^-a^i+^^+^y z(a) - z(0) - y ' w a + 2(0) + v y , w x (3.75) The slice lip deflection can thus be written in terms of the actuators locations (z(0), z(a) and z(L)). As an example assume a = 1, L — 18, k s = 1/18, and EI is varied from around 0.001 to 0.85. The resulting deflections are shown in Figure 3.35. Figure 3.36 shows the detailed deflection and the three actuators working on the beam for the stiffest lip from Figure 3.35. The dotted lines show the actuators at rest when they all have setpoint equal to zero, and the beam has zero deflection. The actuator at x = 1 is then moved one unit and the solid lines show the resulting deflection as well as the actuators compression (at x = 0 and x — 18), and the actuators expansion (at x = 1). If it is assumed for the slice in Figure 3.37 that superposition holds, the deflection at the points where the forces are applied can be calculated as (note that if all forces except p(la) are equal to zero, the equation below is equal to Equation 3.73): u(la u(2a u((Na-2)a)_ or in vector form Q p(la p(2a j>((Na - 2)o). + z(0) + L la 2a .(Na-2)a. (3.76) , , z(L) - z(0) u = Qp+z(0) + v ! w x (3.77) 55 Ciiapter 3: Model for cross-directional basis weight response Figure 3.35 : Deflection of the beam in Figure 3.34 for different stiffnesses, when an actuator at x = 1 is moved one unit. Figure 3.36: Deflection of a beam when a actuator at x = 1 is moved to one from the initial position of zero (shown for the actuator as a dotted line with the endpoints denoted with small circles). 56 Chapter 3: Model for cross-directional basis weight response W *t f C O s i € € L = (N.-I X Figure 3.37: Slice lip. where matrix Q € ^ Na~2)x(N^~2) is symmetric with elements given by 1 / k saz i(N a - 1 - j)(2j(N a - 1) - j 2 - i 2) Q(i,j) = ks { 6EI (N a - 1) (Na-1- j)(N a - ! - » ) + ij (3.78) (Na - l ) 2 for i < j , and u = [«(lo) u(2a) • • • u((JVa - 2)a)]T, p = \p(la) p(2a) • • • p((N a - 2)a)] T, and x = [la 2a • • • (N a — 2)a] . By defining a vector of actuator setpoints z = [z(la) z(2a) • • - z((N a — 2)a)] we can write the forces created by the springs as: p = k s(z - u) (3.79) Then the deflection of the slice can be written as u = ( I + ksQ)' 1 (k sQz + z(0) + ziyL) ~ Z (°K) (3.80) and the forces as p =k s(I + ksQ)-1 ( z - z(0) - ~v~y ~v"yx ) (3.81) A-1 / _ . ./^ *C&) - *(<>) with the values at the endpoints given by Na-2 N _ j _ . JV0-2 p(°) = - Ep ^ a Na-il ajxd p{L) = ~ Ep ( i a ) i v ~ i (3 -82) 1 = 1 != 1 u(0) = z ( 0 ) - ^ and u(L) = z(L) - ^ (3.83) If it is desired to find the slice lip deflection between the actuators, Equation 3.81 can be used to calculate the forces exerted on the slice by the actuators. Equation 3.77 can then be used together with the matrix Q defined for noninteger values of i in Equation 3.78 such that it becomes an iVfine x (Na - 2) matrix where N^ne is the desired number of points for which the slice lip deflection 57 Chapter 3: Model for cross-directional basis weight response is required. The Equations above are useful if the actuator setpoints are known, and it is desired to find the slice lip deflection. On the other hand the slice lip deflection might be known, and it is required to find the actuator setpoints needed for the particular slice lip deflection. Equation 3.73 can be written in terms of the slice lip deflections at the endpoints as: u(x) 1 QEIL b(L2- b2)x + L(x- a) p(a) + tt(0) + u{L) - w(0) If superposition applies as before, the above equation can be written in a vector form as: ^ . , u(L) - «(0) u = Q„ p + u(0) + v ' —— x LJ With the symmetric matrix, Q„ € K0V..-2WA'--2) defined as o c -A a3 t(N a-l-j)(2j(Na-l)-f-i2) 6EI Nn using p = k s(z — u), this can be written as 1 k. z = j-Qu 1{(I + k sQu)u-u(0) u{L) - u(0) L x with the forces given by p = Q - 1 u - u(0) - u(L) - u(0) = QZ - i A-„-2 "A '„ - l ' l L 'Na-l 1 A ^ - l " (0) U = Q ^ A U The forces for the whole slice, including the endpoints can be written as _ _ _ l V _ 1 ' M p = LP&) 'A' - 1 I The bending moments are then written as: 'M(O) M M{L) AV-T ' 0 1 Na'- 1 Nq-2 ' N a-1 0 P = A r p p(0) p(la) p((Na - 2)o) • P(L) (3.84) (3.85) (3.86) (3.87) (3.88) (3.89) = S P (3.90) S A r p = SATQ,; 1AU = SA rQ,71Ax^K f . Accordingly the bending moments can be calculated directly from the parameters of the Gram polynomials of the slice, which allows limiting the bending moment by modifying the parameters. Note that all the matrices in Equation 3.90 can be precomputed, as they only depend on the slice lip configuration. 5 8 Chapter 3: Model for cross-directional basis weight response 3.3.1 Bendin g limi t In general, the Gram polynomial approximation of a specific profile can be found using least-squares estimates, according to Equation 2.36. The least-squares solution can result in a shape profile that is not possible to achieve with the slice lip since the stroke is outside the actuator's range or the bending limit of the slice will be exceeded. To satisfy those restrictions, the square of the approximation error (e(t) from Equation 2.35) can be minimized subject to restrictions on the maximum and the minimum value of the slice lip position as well as the minimum and the maximum bending moments. The minimization of e2(t) (i.e. least-squares) can be transformed into a quadratic programming problem according to min \\u - x^auc | | = min (u - x^ auc) (u - xjfaur) "c (3.91) = min (u Tu - u^x 0au - (xQ au)Tuc + «Jxoaxjfaucj The value of u Tu does not have any effect on the minimization with respect to u c, therefore, it can be omitted. Also the minimization of uJxoaU is the same as the minimization of (xo nu) u c which can be seen from the fact that (uJxoatt) = (xo au) u c. The above equation can then be written as min \\u - x^-uJI = min (ttJxoax^,«c + 2(-x 0au)Tuc) Uc ' U c \ / = min I -uJx0aXQauc + {~XQ nu)Tuc subject to Ruc < B (3.93) which can be solved using Matlab's optimization toolbox [55]. The restrictions can be written as maxM —SA?Q~lAxQaUc < - m i n M (3.94) x0a u r < max u —x0auc < — mmu 3.4 Summar y In this chapter the paper machine has been described and a model that relates slice lip movement to basis weight response has been developed. A model has also been given that can be used to 59 Chapter 3: Model for cross-directional basis weight response describe the bending of the slice lip. The main feature of the model is an interaction matrix that describes how the movement of a single actuator will effect a large section of the measured profile. The effect is transmitted through bending of the slice lip which sets up waves on the forming wire and spreads out on the wire. The interaction matrix can, for many paper machines, become singular causing problems when controlling the basis weight. Some of the problems in creating the interaction matrix are discussed and it is found that it is difficult to create a matrix that will not create ripples, particular to the model, in the model output. Gram polynomials are then used to parameterized the interaction matrix, resulting in a much smaller matrix which can be identified directly and will be discussed in the next chapter. The form of the parameterized interaction matrix is discussed and it is found that for a general interaction matrix the parameterized matrix will have nonzero entries mainly on the diagonal. A model of the slice lip as a flexible beam with multiple forces is given and the slice lip deflection, which is written as Gram polynomials, is related to the bending moment. The parameters of the Gram polynomials can therefore be manipulated to limit the bending moments for the whole slice. 60 Chapter 4 Identification Basis weight, which is measured at the finishing section of the paper machine, is controlled by the slice lip actuators located on the headbox. For control purposes it is necessary to know the transportation delay from the slice lip to the measurements. The delay can be estimated by carrying out bump tests on the slice lip, or by knowing the distance the paper travels through the machine, and the speed of the machine, it can be calculated. A knowledge of the interaction matrix, G, which describes how each actuator affects the basis weight, is also required. The interaction matrix of a paper machine depends on the speed of the machine, the target basis weight, the consistency of the pulp in the headbox and other variables. Until such time as a model is available that will calculate the interaction matrix given all the miscellaneous variables, the interaction matrix has to be estimated. If Gram polynomials are used, the parameterized interaction matrix, G x (see previous chapter) could as well be estimated. Currently no paper machine data is available from which the Gx matrix could be identified. Therefore, a paper machine was simulated using an interaction matrix identified from industrial data. The identified matrix was then used as a model test-bed to identify the G x matrix. There are two ways to identify the Gx matrix from those simulations. One way is to use the estimated G matrix and then calculate G x using Equation 3.55, i.e. G x = (XOXQ1)- XOGXQ U. Another way is to simulate a paper machine using a particular interaction matrix and then identify directly the G x matrix, from the parameters of the orthogonal approximation of the measured profile, y c, and the parameters of the orthogonal polynomials of the slice lip, uc. In this chapter the latter method, which identifies the G x matrix directly, will be used. 4.1 Identifying the interaction matrix There are many ways to estimate the interaction matrix [22, 64]. Basically it involves perturbing several actuators and using least-squares estimation or correlation studies to find the steady-state gain from the actuators to the measurements over an area of the measured profile corresponding to the perturbed actuators (see Figure 4.38). For a single actuator, several estimations might have to be run, one for each measurement, and then from those estimations the response shape could be obtained. The experiments are then repeated for several actuators across the machine. The response shapes 61 Chapter 4: Identification from those several actuators are generally averaged to get a single response shape to represent the response to any actuator across the machine. As was mentioned in the previous chapter, there might be some problems with mapping the measurements to the actuators. Figure 4.38: An experiment to identify the interaction matrix, G. Experiments were done on a 8 m wide linerboard paper machine, with target weight at 61.1g/m2 [64, 65, 77]. The 105 slice actuators are located 76 mm apart and there are 445 measurements across the machine with data box size around 18 mm. The data were collected for 60 scans of the profile, and moisture and basis weight profiles were recorded, as well as the slice lip setpoints and basis weight profiles mapped to the 105 actuators. It was observed for this particular paper machine on this particular day that an actuator could affect an area of the paper sheet up to nine actuators wide. Therefore, in order to minimize the effect of cross-coupling on the identification, every ninth actuator is perturbed with a PRBS signal, starting away from the edge, at actuator number eight. Hence 11 actuators are perturbed with PRBS signals, at the CD positions given in the first column of Table 4.1. Figure 4.39 shows a black and white image of the actuators and the measured basis weight, with the mean value of every CD point subtracted in order to better see the response. From the black and white image, the area of the measured profile, corresponding to the 11 actuators can be found. The model describing the effect of each of the inputs on its corresponding area can then be estimated. The centers of response as seen from the input-output data in Figure 4.39 are assumed to be the measurement points corresponding to the center of response. Those measurement data boxes are shown in Table 4.1 and will only be used to limit the area where the steady state gain from 62 Chapter 4: Identification Actuator number 8 17 26 35 44 53 62 71 80 89 98 Measurement number 27 64 105 144 182 220 259 298 339 377 416 Table 4.1 Location of the PRBS inputs and the corresponding output measurements deduced from Figure 4.39. the actuators to the measurements is found. The exact location of the center of response will be calculated later from the estimated responses, by calculating the cross-correlations between the 11 response shapes. Figure 4.40 shows input and output data for one CD location. It can be seen that the dynamics are very fast, indeed the system reaches steady state within one scan, after a pure delay of two scans. The delay can also be established by looking at the cross correlation between the 11 PRBS signals and their corresponding outputs. Figure 4.41 shows the cross correlation between the 11 PRBS signals and their presumed center of response. The dynamics of the system are negligible therefore, the maximum of the correlation function indicates the delay of the system. It is evident that the maximum correlation is reached for a delay equal to two, for all the 11 sequences. The same delay is obtained when finding the maximum of the absolute value of the correlation between the 11 inputs and their closest measurements. Hence the knowledge of the delay will be assumed to be known in the estimation that follows. From the knowledge of the process delay and dynamics, the process can be approximately modeled by y(t,x)=pou(t-2,x) (4.95) 63 Chapter 4: Identification 60 50 40 3 3 0 20 10 • 1 1 1 1 • 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 • 1 1 1 1 • 1 • 1 • 1 I 1 1 1 1 20 30 40 50 60 70 80 90 100 Actuator 60 f*„9* ? w > v ' . »-» saw f i '- t »•*•*-* "jar •• m, •—'-—« • » - - - r —~ »> ••<•••»*• , . •< 20. v . ' ' i " S T . . T , I _ * «• !* ' - • •*"• • -* « - . ' < M" . % 1 A * * i * ^* ., , , < J - i t* - - " fc» , ^ , , . ** ' ' ' , < < • i « i * * * A m ' * - Sr 50 100 150 200 250 300 350 400 Measurement Figure 4.39: Slice lip actuators and basis weight measurements from an industrial data set. To get a good contrast in the picture the mean value of every CD point is subtracted from the raw data. where y(t,x) and u(t,x) are the basis weight and actuator set point respectively, at time t and location x across the sheet. Each input will affect several outputs, which can be written as \y(t,x-k) y(t,x-k + l) ••• y(t.x + k)} = [p_* p_A .+ 1 ••• Pl ,]u(t-2,x) (4.96) Since the basis weight and actuator setpoints have a nonzero DC level, it is necessary to look at the incremental inputs and outputs, to estimate all the Pj's: At[y{t,x-k) y(t,x-k + l) ••• y{t, x + k)} = [p_ t p_A-+] ••• p k.]Atu(t - 2,x) (4.97) 64 Chapter 4; Identification 0 10 20 30 40 50 60 Scan Figure 4.40: Slice lip set point and measured basis weight as a function of time, from an industrial data set. Figure 4.41: The cross correlation between the input and the output for the 11 PRBS signals from an industrial data set. All the 2k + 1 estimations for each input, can be done separately with A (y(i ,x - k) = y(t,x - k) - y(t - l,x - k) = p-k{u(t - 2,x) - u(t - 3,x)) (4.98) which can be identified using batch least squares, recursive least-squares or recursive instrumental variables [87]. Recursive instrumental variables are used for these data to get unbiased estimates. Only the gain of the system is estimated and the system is assumed to be time-invariant during the 65 Chapter 4: Identification 0.6 0.4 S3 0. 2 > 4J 0 6 w -0. 2 -0.4 0 -0.6- | / I ' I \ / 10 20 30 Scans 40 50 Figure 4.42 : Estimates of steady state gain, pj, between the slice lip actuators and the corresponding measurements from Table 4.1 for an industrial data set. 60 150 200 250 CD Measurements Figure 4.43 : The estimated response across the sheet, for 11 actuators from an industrial data set. The plus signs mark the final estimates from Figure 4.42. trial, therefore a forgetting factor of one is used. Figure 4.42 shows the 11 estimates using the slice lip actuators and measurement data boxes in Table 4.1. Those are 11 independent estimates, shown here to see if all the 11 inputs are sufficiently rich to give good estimates. It can be seen that the estimates converge within the 60 scans, but they do not all converge to the same value, because the outputs in Table 4.1 do not necessarily correspond exactly to the inputs. The final estimates are shown in Figure 4.43 with the plus signs. To obtain a cross directional model, the responses for moving an actuator is estimated at 20 measurement locations on both sides of the corresponding measurement location from Table 4.1, using recursive instrumental variable as before. Figure 4.43 shows the 11 response shapes from the estimation. In order to average the 11 responses it is necessary to find the number of data boxes 66 Chapter 4: Identification X> S 3 C c 'So BQ & 0* / 2 3 4 5 6 7 8 9 10 11 PRBS signal number i 1 0 -38 -78 -117 -156 -194 -233 -272 -312 -351 -390 2 38 0 -39 -79 -118 -156 -195 -234 -274 -313 -352 3 78 39 0 -39 -78 -117 -156 -195 -234 -274 -312 4 117 79 39 0 -39 -77 -116 -155 -195 -234 -273 5 156 118 78 39 0 -38 -78 -117 -156 -195 -234 6 194 156 117 77 38 0 -39 -78 -118 -157 -196 7 233 195 156 116 78 39 0 -39 -79 -118 -157 8 272 234 195 155 117 78 39 0 -40 -79 -118 9 312 274 234 195 156 118 79 40 0 -39 -78 10 351 313 274 234 195 157 118 79 39 0 -39 11 390 352 312 273 234 196 157 118 78 39 0 Table 4. 2 Number of data boxes between the maximum values of the cross correlations between the 11 responses in Figure 4.43. 0 .41 1 1 •— i 1 1 1 1 r CD measurement positions relative to the peak of response Figure 4.44 : The 11 estimated response shapes from Figure 4.43 shifted according to the first line in Table 4.2 such that they overlay. 67 Chapter 4: Identification -15 -10 -5 0 5 10 15 20 CD measurement positions relative to the peak of response Figure 4.45: Average estimated response of the weight profile by moving a single actuator. between the responses. The cross-correlation between all the responses is calculated and the number of data boxes between the maximum of the cross-correlation is shown in Table 4.2. It can be seen that the number of data boxes between the two outermost responses, the PRBS signals number one and 11, is always 390 data boxes, no matter which PRBS signal is used as a reference. The distance between the second response and the third response is between 39 and 40 data boxes, depending on which response shape is used to calculate the difference. The first line in the table, i.e. the correlation between the first response shape to all the other responses, gives the difference between the peak of response that is very close to the mean value, therefore it is used to align the responses on top of each other to find the average response. A plot of the results are shown in Figure 4.44. The response shapes are averaged, giving a typical response shape for a paper machine as can be seen in Figure 4.45. The figure shows the mean and the standard deviation for the 11 responses from the previous figure, as well as the actuator spacing. It can be seen that beyond 16 measurement zones from the center of response, there is almost no response to the actuator bump. Therefore it will be assumed that the width of response is 33 measurement zones. To make the response symmetric the left and right hand side in Figure 4.45 are averaged, hence 17 parameters will be used in the interaction matrix. The block diagonal interaction matrix, G, is constructed by putting the averaged response in every column of the matrix, centered in the rows corresponding to those columns. As mentioned in Chapter 4: Identification Chapter 3, it is difficult to decide in which row to place the center of the response, especially if the number of measurements is not an integer multiple of the number of actuators or if the sheet wanders. 4.2 Identifyin g th e parameterize d matri x As previously mentioned, once the G matrix has been estimated the Gx matrix can be calculated using its definition (see Equation 3.55). Another way, which does not depend on the knowledge of the G matrix, is to estimate the matrix directly from the knowledge of the parameters of the orthogonal approximation of the measured profile, y c, and the parameters of the orthogonal polynomial of the slice lip, u c. To simplify the estimation, the structure of the G matrix could be used. But doing so would defy one of the reasons for using orthogonal polynomials, i.e. not having to establish the mapping from an actuator to the measurements. For that reason every element of G x has to be estimated. The problem is therefore MTMO (multiple input multiple output) estimation. A MIMO system can be written as [86 p.80] Iyc{t) + Aiyc(t - 1) + • • • + An ayc(t - n a) = B xuc(t - d) + - • • + B Ubuc(t - d - n h + 1) (4.99) where y c(t) is a vector of length m, Ai is a matrix of size m x m, u c(t) is a vector of length m a and Bi is a matrix of size m x m a. Keeping the analogy with SISO systems, a parameter matrix 9 can be defined as e = [A x A 2 ••• A n. Bj ••• B nb (4.100) with the size as [n am + ribma] x m. The regression vector ip is a vector of vectors and is written as: ¥>(<) = -Vc(t - 1 -Vc(t - 2 -yc(t- n a) uc(t — d) (4.101) uc(t — d — ni, + 1). and will have length of [n am + nbma}- The output vector y c can then be written as yc(t) = 0Ttp(t) (4.102) It has been shown that a paper machine response is governed by a delay and a gain. The model for the basis weight response is not found to have any zeros, so re& = 1. The time constant for the 69 Chapter 4: Identification system is very small and can be ignored in most cases, so «„ < 1. Furthermore, the matrix A\ can be assumed to be diagonal, i.e. the outputs do not affect each others. Every row of the Gx matrix can then be estimated separately, so the estimation has been reduced to m MISO (multiple input single output) estimation problems, each estimating m a + na parameters. The m systems can be written as: yc(t,i) - [ai Bi(i,:) -yc(t-l,i) uc(t - 2) (4.103) where a,; is the pole for the ith input, Bi(i,:) is the ith row of the B\ matrix and y c(t,i) is the ith element of y c(t). If there is no dynamics, for each MISO system, i, where i € [ 1 , . . . ,m], the system is written as: yc(t,i) = [G x(i,l) G x(i,2) ••• G x{i,ma)] * ( * )>(<)• uc(t — d, 1 uc(t — d, 2 uc(t - d,m a)_ (4.104) The G x matrix can be estimated using standard least squares or recursive least squares [86, 87, 113]. If the system has some dynamics (for simplicity the dynamics is assumed to be the same for all the polynomial orders), for example Vc(t) 1 + a —^Gxuc(t - d) 1 + aq the G x matrix can be estimated from the estimates of B\ and a such that Gx = 1 + a (4.105) (4.106) 0.4 0.3 S0.2 0.1 0 445 points 420 points 0 50 100 150 200 250 CD measurements 300 350 ' I 400 Figure 4.46: Basis weight response to a flat slice lip (at 10/zm) with 105 actuators for two different models of a paper machine. 450 70 Chapter 4: Identification For the simulation studies the plant is assumed to be described by the interaction matrix from the previously identified G matrix except here it is assumed to be 420 x 105, i.e. the number of measurements is assumed to be an integer multiple of actuators. The reason for doing so can best be seen in Figure 4.46 where a flat input gives (apart from the edges) a flat output if the number of measurements is an integer multiple of the number of actuators. On the other hand if the number of measurements is assumed not to be an integer multiple of the number of actuators, the model output is very jagged as is evident from the solid line in Figure 4.46. This jagged profile is particular to the model chosen, where all responses are assumed to be the same and the center of every actuator is assumed to be aligned with a particular measurement data box, which is impossible for a noninteger ratio of actuators to measurements. Therefore there will be multitudes of rows and their sum will take on varying values, if all the response profiles to each actuator are assumed to be the same. It is therefore very difficult to fine tune a model to accurately describe a paper machine. It can be somewhat simplified by use of a number of measurements that is an integer multiple of the number of actuators. The slice lip profile used for the simulation is controlled by 10th order Gram polynomials and 10th order polynomials are used to approximate the measured profile. Independent PRBS signals of magnitude ±0.1 are added to the parameters of the polynomials of the slice lip, the paper machine is simulated, the basis weight profile is measured and the Gram polynomial parameters of the simulated output profile are calculated. Figure 4.47 shows the parameters of the polynomials of the slice lip, the input u c, and the parameters of the Gram polynomial approximation of the simulated output, y c, used for the estimation of G x. The PRBS signals added to the parameters of the slice lip results in the slice lip deflection of at most ±0.6 • 10//m = ±6/zm (see the top graph in Figure 4.48) from its initial settings. Therefore the peak to peak deflection of the slice lip required by the excitation signal is 12/x m which is only about a quarter of the peak to peak value of the slice lip before applying the PRBS signals. The slice lip bending, that results from moving the slice lip for the estimation trials is shown at the bottom of Figure 4.48 with length of the slice lip, i = 8m, and the bending modulus, EI = 4500 Nm2. The dotted line in the graph shows the bending limit that will be imposed upon the slice in the next simulation. The measured simulated basis weight profile will then have a deflection of at most ±0.24 g/m2 (Figure 4.48) from the initial measured profile, which amounts to about 50% of the peak to peak value of the initial basis weight profile before applying the PRBS signals. This input-output data is used to estimate the Gx matrix using 10 independent recursive least squares (RLS) on Equation 4.104 with i, the row number, between one and 10 inclusively. Each of 71 Chapter 4: Identification 1.5 0.5 o 3 o--0.5 -1.5 . • • • 1 -i- . 1 . • • • 1 1 - . • • • • % • • • '• •' : ' • • ' - . • ' • • -2- \ / \ _'~\_ 3 - " \ / ^ "\ 7 - - - \ / " \ -9-^ ^ \ _ / W * v _ ' " ~ X _ / \ / ~ \ fry---. •-/•- • v- - • 1 ' 1 / \ , \ -_ . _ / " x - . . " X - /" " -_^={ / _FX r ^ p^ — ^- -* \ — / \ _ / \ _ _ • / \ / • • • • \ 1 1 1 10 20 30 40 50 60 Figure 4.47 : Input, u c, and output, y c, used for the estimation of G x, when there are no limits on the slice lip bending. 72 Chapter 4: Identification o u . 3 IS o a ^ T3 i r 3 C t> 6 o 6 60 C c « 100 150 200 250 CD measurements 300 350 40 50 60 70 Actuator number 80 90 10 20 30 40 50 60 70 80 90 100 Actuator number 400 Figure 4.48: Overlay of 59 slice lip profiles, measurement profiles and bending moments corresponding to the inputs and outputs in Figure 4.47. the 10 independent estimators, estimates 10 parameters, one for each input. Figure 4.49 shows the estimated parameters from the 10 estimators, one for each row. It can be seen that the parameters converge in about 40 scans. The estimates at scan 59 are used to make the G x matrix. It can be seen in Figure 4.50 that the estimated matrix is the same as the true matrix. At least there is no visible difference between them. 4.2.1 Limite d bendin g o f th e slic e li p As said earlier, the estimation of the parameterized matrix can be performed as long as the regressors, ip(t) in Equation 4.104 are linearly independent When the bending is limited, the regressors will become correlated, because if bending limits are exceeded, all the parameters of the Gram polynomials will have to be adjusted and will therefore be correlated. For that reason the estimates might not converge to the true values. Also a limit on the slice bending will result in the higher order parameters being limited the most, because higher order polynomials have higher bending moments. Therefore the excitation for the higher order polynomial will be severely limited, resulting in difficulty in estimating the corresponding row of the G x matrix. 73 Chapter 4: Identification Rowl Row 2 0.3 0.2 0.1 0 0.4 0.2 0 -0.2 0.2 0.1 n A 1 ^ Y '• H/Li^- f l \k- '• : ^* "4 i l l! A FL ^ 0 20 40 60 Row 7 0 20 40 60 Row 9 • • • j r r ^ 1 , • ^ • ^ t , w , M , ^ r , | - | - h * M ' - ' - , ~ — 0 20 40 60 Scans Figure 4.49 : Estimation of each row of G x when no limit is put on the slice lip bending. The true parameters for this simulation, are shown with dotted lines. 74 Chapter 4: Identification Figure 4.50: Estimate of G x after 59 scans when no limits are put on the slice lip bending. The estimated matrix is denoted with * and dotted lines between its elements. The true matrix used in the simulation is denoted with o and solid lines. Figure 4.51 shows the input and output used for the estimation, when the bending moments are limited at ±0.23 Nm. It can be seen that the input parameters are correlated and they also have smaller values than when there was no limit on the bending moments. Because of smaller input parameters the output parameters will also be smaller as can be seen in the same figure. Figure 4.52 shows the slice lip profiles, measured basis weight profile from a simulated paper machine and the bending moments of the slice lip. The slice lip is now smoother than before (compare to Figure 4.48), especially at the edges. It shows clearly in the plot of the bending moments, where they are successfully limited at ±0.23 Nm. Figure 4.53 shows the estimated parameters of the G x matrix. It can be seen that the convergence is slower than before, and for the higher orders, the parameters of the corresponding row do not converge within the 59 scans. Figure 4.54 better illustrates the estimates of G x after the last scan. Now the difference between the true matrix and the estimated one is clearly visible, mainly for the higher order terms. 75 Chapter 4: Identification 1.5-0.5 0--0.5 -1 -1.5 -- * • < _ 3 - r 7 y 6" N _ /-~J~ / ' i ^ •s ~ \ t " / "" / / ~ =/ZT. /* — / \ ~ ' \ ^ V \ - ~ - s » ' • • \ 1 • •/. -' • V _ . -/ -^IX^ \ 1 . ./ - . — . — ^ r - / - - — / i -/--.-. \ / t - . -J ^ ^ = ^ V 1 . r J . -^~_ \ i "\ . — . - "~ 7 i / r\^ f. 7 • *v» - . _ _ \ • \ --"V ^ _ . -— 10 20 30 40 50 60 Figure 4.51: Input, u c, and output, y c, used to estimate G x, when there are limits on the slice lip bending. 76 Chapter 4: Identification 50 100 150 200 250 CD measurements 300 350 10 20 30 40 50 60 70 80 Actuator number 400 90 100 Figure 4.52: Overlay of 59 slice lip profiles, measurement profiles and bending moments corresponding to the inputs and outputs in Figure 4.51. 4.2.2 Estimatio n o f G x includin g dynamic s The paper machine might have some dynamics other than the delay. For example there might be that only the filtered profile, i.e. the profile out of the CD estimator, is available. Then the model given in Equation 4.105 can be used and one additional parameter, a will be estimated for each row, and the pole averaged to give the pole of the system. For this simulation the filter pole is set at 0.2 or a = —0.2. The same input as previously will be used and no limit will be used for the bending moments. The inputs and outputs are shown in Figure 4.55. The dynamics are clearly visible in the output. The estimates of G x do not seem to be affected by estimating one more parameter for each row as can be seen by comparing Figure 4.56, the estimates of G x when the filter pole is estimated as well, to Figure 4.49, the estimates of G x without the filter. The estimates of the filter pole converge quickly to the correct value of —0.2 as can be seen in Figure 4.57. 4.3 Identifyin g th e invers e o f th e parameterize d matri x In the controller design (see Chapter 5), use is made of a matrix KQ, which is calculated from 77 Cliapter 4: Identification 0.2 0.1 o Row 1 [~ / " " I 1' ' Gi in' '"' ' ^'i~'~n''' iqj* Row 2 -0.1 Row 9 -0.1 0.4 0.2 0 20 40 Scans 60 20 40 60 Scans Figure 4.53 : Estimation of each row of G x when limit is put on the slice lip bending. The true parameters for this simulation, are shown with dotted lines. 78 Chapter 4: Identification Gx Figure 4.54: Estimate of G x after 59 scans when limits are put on the slice lip bending. The estimated matrix is denoted with * and dotted lines between its elements. The true matrix used in the simulation is denoted with o and solid lines. GXKQ — I, where / is a m x m matrix. When G x is a square matrix, the inverse of G x is equal to KQ. Equation 4.105 can then be written as uc(t - d) = ——yc(t) + — —yc(t - 1) (4.107) 1 + a y ' 1 + a' Instead of having to calculate the inverse of the estimate of G x, which might be close to singularity, the matrix KQ could be identified directly. The identification of KQ can only be done in open loop, because in closed loop the controller would be identified as well as the KQ matrix. From Equation 4.104, for i 6 [ 1 , . . . ,m a] we can write: ic(t-d,i) = [B 1(i,:) B 2(i,:)] = 8(if<p(t). Vc{i) Vc(t - 1) (4.108) From the above equation the dynamics and the inverse of G x can be calculated as KQ = B\ + i?2 a = B 2/B1 (4.109) 79 Chapter 4: Identification 1.5 1 0.5 0 -0.5 -1 -1.5 2- N 3 - " \ _ / \ ^ j v / • w . / "v. . / 9 ^ " Q^IH^K^ / \ _ / J7M A / \ ' / \ •/^v-. • / \ r~\-10 20 30 40 50 60 Figure 4.55: Input and output used to estimate the G x matrix when the dynamics is estimated as well. 80 Chapter 4: Identification Row Row 2 0.2 0.1 0 : j^ ~^33feP^ Figure 4.56 : Estimation of the parameterized interaction matrix G x when the dynamics is estimated as well. The dotted lines show the true values. 81 Clmpter 4: Identification Figure 4.57: Estimation of a, for the 10 rows in the previous figure, where the division is an element by element division. Another way to estimate KQ is to assume knowledge of a. The input u c can then be sent through the filter, with pole at o, and the filtered input used to identify the inverse matrix. The filtered input is therefore or , . 1 + a . yc(t) = G xu'c(t) (4.110) (4.111) uc(t) = K 0yc(t) (4.112) Therefore the KQ matrix can be identified using output and filtered input, provided the filter is known, or Equations 4.108 and 4.109 used. Note that the elements of the vector y c are correlated as is evident from Equation 4.104, where a single input, u c(t,i), will in general affect all the outputs, y c(t). Therefore it can be expected that the estimates of KQ will be slower to converge than for the estimate of G x. Figure 4.58 shows the estimates of the 10 rows of KQ using the same input and output data as before (see Figure 4.47). Comparing those results with Figure 4.49 it can be seen that the convergence is about 10 scans slower for the estimation of KQ. Nevertheless the estimates do converge to the true values and there is no visible difference after 59 scans between the estimated matrix and the true matrix used in the simulation as can be seen in Figure 4.59. 82 Chapter 4: Identification Row 1 Row 2 20 4 0 6 0 Row 10 Figure 4.58: Estimation of each row of KQ, the inverse of G x. The dotted lines show the true values. 83 Chapter 4: Identification KO Figure 4.59: The estimate of KQ after 59 scans. The estimated matrix is denoted with * and dotted lines between its elements. The true matrix used in the simulation is denoted with o and solid lines. 4.4 Summar y In this chapter the interaction matrix, G and the parameterized interaction matrix, G x have been identified. It is shown that even though a response shape can be identified that is used in creating the interaction matrix, it is not trivial to construct the interaction from the response shape. The identified G matrix, with number of measurements as integer multiple of the number of actuators, is used to simulate a paper machine to use for identification of the parameterized matrix. Few examples are given on the identification of G x, with or without bending limits. An example is also given on how to identify the inverse of the interaction matrix, which might be helpful in the control design. 84 Chapter 5 Basis weight Cross Directional Control The main problem with controlling basis weight on paper machines is the cross coupling in the cross direction. It will be dealt with by the use of Gram polynomials introduced in chapter 2. A significant dead-time is present in many pulp and paper processes. The cross directional control of basis weight has been shown in chapter 3 to have transportation delay. The Dahlin controller, which is a dead-time compensator, is used because of its simple design and ease of use. 5.1 Dahli n controlle r an d integrato r antiwindu p The design of a Dahlin controller states that the closed loop response of the feedback system should be first-order plus dead time, with unit gain and a closed loop dead time equal to the open loop one [38]. The closed loop time constant is the only tuning parameter. Let a discrete plant, HOL(<1~ 1), be described by H°L{q ] ~-A(^)q - i + a i r i + . . . + a„ a ( r». q (5-113) With the use of a feedback, it is possible to design a controller, C = A c)q-il ( s e e Fi§u r e 5.60), such that the closed loop system is a first-order process plus a delay with steady state gain equal to one. It is written as: Hct{q-l) = - y —q-d (5.114) 1 — pq i where p, the tuning parameter, is the closed loop pole. The closed loop system in Figure 5.60 is then , n HOLC BB C _ d 1-p _ d / c 1 1 C s HCL{* ) = TTH^C = AAc + BBcq-*q ^Y=W q (5 '115) It is then straightforward to show that the controller becomes: C(Q-*) = ^ -1 = !-* Ai"-1) (5 il 6 ) [Q } A c(q-1) l-pq- 1-(l-p)q-dB(q-i) P ' j This is the Dahlin controller and it inverts the process dynamics, i.e. it cancels out all the process poles and zeros. The inverse of the process has to stable which means that all the process zeros have 85 Clmpter 5: Basis weight Cross Directional Control Be A, « * ! H OL Figure 5.60: Dahlin controller for a plant, to be stable. In practice it is not possible to know the poles and zeros with infinite accuracy. The unstable process poles can therefore not be cancelled with unstable zeros, hence the process has to have stable process poles as well. The process zeros have to be well damped so as not to cause ringing in the control signal. The model used for basis weight control is a stable first order model with gain and a delay so problems with poorly damped or unstable zeros and unstable poles is not an issue. In the slice lip controller used in the industry today, control action is typically taken every three scans. The Dahlin controller can then be designed based on the steady state value of the process, with the implicit assumption that the process reaches steady state within the three scans. Before we make that assumption, lets look at the Dahlin controller with control taken every three scans. Assume the model is written as a first order process with delay of two, as The model of the plant that the controller sees (every three scans) is then H0L{q7l) = ^ i + Zp + yfrr1) - i q7Mk) (5.117) (5.118) 1 - ffc1 where q~ 2 is used to denote sampling every three scans. The controller, assuming the closed loop system is ^~*lig~1, is then -pg3 (5.119) If the dynamics of the process has also died out within the control interval, i.e. fp is close to zero, the controller is then basically a steady state controller that takes into account the delay of the system. The controller is then written as C(i7') = 1 - 1 1 - P 1 - / P 7 - 1 (5.120) 1 - 9 s -i 1 - q~i b and is calculated every three scans. The control signal, which is calculated at k = 1,4,7,..., is then: u(k) = fr-^-^eW + uflb-S) (5.121) 86 Chapter 5: Basis weight Cross Directional Control Its performance can be seen in the upper graph in Figure 5.61 for f p = l — b = 0.1 and p — 0.05. It can be seen that the system does not quite reach steady state within the three scans and consequently there is a small overshoot in the output. To compare its performance to the performance of a control taken every scan, the same plant is used, but now the closed loop pole is set equal to the cubic root of 0.05 or p = 0.37, which gives the same time constant as 0.05 for control action every third scan. The response for the control taken every scan is shown in the lower half of Figure 5.61. It can be seen that the output is almost the same as for control action taken every third scan, except there is no overshoot. Much more frequent control action is taken in this latter case and not much is gained on performance for a system with such a fast dynamics. On the other hand much faster response is possible when control action is calculated every scan and hence it will be used in what follows. -1 * r ~-M - - - * -ydes Jt " _ - * - - J * - * - * -- } ( 10 15 20 25 -1 * • ydes y u L - * - - J K --•__ j- - ?K * -* -10 15 Samples 20 25 Figure 5.61: Dahlin controller with control action taken every 3 scans with the closed loop pole, p = 0.05 and the open loop pole is, f p — 0.1. Dahlin controller with control action taken every scan with the closed loop pole, p = 0.37 and the same open loop pole. The * denotes the output at every third scan. An example of a Dahlin controller for a plant, HQL — i 7 -i q~ 2, is according to Equation 5.116 c = (5.122) i -pq 1 - ( l -p)q 2 The plant has a steady state gain of one so the control signal will have to be at the same level as the output, for the output of the plant to follow the reference signal. Figure 5.62 shows the reference 87 Chapter 5: Basis weight Cross Directional Control signal as well as the output and the control signal for a plant with a pole at f p = 0.3 and with the closed loop pole chosen as p = 0.05. It can be seen that the control signal required to achieve the closed loop response is always greater than or equal to one. Now if the control signal is bounded at ±0.5 it is observed that the controller performance deteriorates (see middle plot in Figure 5.62), especially after long periods of constant input it takes a long time for the controller to respond. The reason is the integral action of the controller, it integrates an ever present error, which results from limited control action. When a new step is applied to the system the controller has to integrate (wind back), out of the saturation. It can be seen in the middle plot in Figure 5.62 that because of the integral action it takes considerable time for the control signal to respond to changes in the input. For example the step change applied at sample 13 is not seen in the output until at sample 19 instead of sample 15. If there would be a change in the reference signal (t/des) of short duration during the integrator wind back, it is possible that its effect will not be seen in the output. It is possible to guard against this windup [6]. The controller, C, which is written in terms of its nominator and denominator as ^ or Acu = B ce (5.123) The bounded control signal can be observed, with the aid of an observer A>- Adding the observed signal to both side of the previous equation results in A0u + Acu — B ce + A0u (5.124) or A0u = B ce + (A0 - Ac)u (5.125) A controller with antiwindup compensation (see Figure 5.63) is then given by [6] AoV — B ce + (A 0 — Ac)u (5.126) u = sat{t)} where sat{} is the saturation limit on the control signal. When there is no saturation this controller is equivalent to Equation 5.123, when it saturates it can be interpreted as an observer with dynamics given by A 0. Figure 5.62 shows clearly that the performance of the controller with integrator antiwindup in place is superior to the case when no integrator antiwindup is in place. With antiwindup in place the controller responds immediately after the time delay of the process, and the output follows the input. The output will not be able to reach the reference signal because of the limited control action. 88 Chapter 5: Basis weight Cross Directional Control Desired input-output response 1 1 L _ 1 l ydes y 1 1 1 1 1 1 1 1 1 1 1 _ — • 10 15 Without integrator antiwindup 20 25 10 15 With integrator antiwindup 25 -" — -• -1 - -, _ f -_^- . 1 1 V u 1 -1 1 _ 1 1 1 1 -— 10 15 Samples 20 25 Figure 5.62: The signals in a closed loop control with a Dahlin controller; with and without limits on the control action, and with and without integrator antiwindup. The plant pole is f v = 0.3 and the desired closed loop pole is p — 0.05. <> + I Efe ' A„ A-A 0 c F - u t\ Figure 5.63: Controller with antiwindup. 5.2 Contro l o f basi s weigh t The basic control strategy for basis weight control using Gram polynomials is shown in block diagram in Figure 5.64. With the use of the pre-compensator, KQ, the controller can be designed for a SISO feedback system with controller C and plant dynamics which is the same as the dynamics from u c to yc. There is a delay in the model of the plant, therefore a Dahlin controller is used. For a 89 Chapter 5: Basis weight Cross Directional Control plant written as B{g -d _ —g d the Dahlin controller is written as Equation 5.116. Limits will have to be placed on the control signal in order to protect the slice lip from excessive bending and the stroke of the slice lip is also limited. When the slice lip is in saturation the controller has to know, else the signal out of the controller might windup to a bigger and bigger value. To protect against the windup [6] the controller can be designed using Equation 5.126. For B{q~ l) = b and A(q~l) = 1 — f vq~l the controller can then be written as Acvc(t) = K Q-j£[e(t) - f Pe(t - 1)] + \ A O - ( l - pq~ l - (1 - p)q- d)]uc{t) (5.127) For d = 1 and a dead beat observer, A„ = 1 the controller can be written as vc[t) = K 0^-^[e(t) - f pe(t - 1)] + u c(t - 1) (5.128) and for a delay equal to two, the controller is vc(t) = K 0^Y?-[e(t) - f Pe(t - 1)]+ Puc(t - 1) + (1 - p)u c(t - 2) (5.129) Figure 5.65 shows the block diagram of the controller with antiwindup in place. vVty *o + * 1 , ' 0 i ! ~\ i Est . C — 1 ? "! yc -*• q 2B - (XoXj)'x 0 - -Figure 5.64: Block diagram of basis weight control strategy using Gram polynomials. A - A H . 0 c Jdes . . -jy! H x x r x o r - ^ - ^ B C K. h-.* < bending •- x „ ,-2B - i Est . p : L_. (X^jX^) x 0 Figure 5.65: Block diagram of basis weight control strategy using Gram polynomials with an integrator antiwindup. The controller, C, from Figure 5.64 is 4 (^-i) • 90 Chapter 5: Basis weight Cross Directional Control Estimation results from the previous chapter show that G x is close to a constant times the identity matrix. Here we use this information to initialize the estimated matrix. First the system is run in open loop to get a better estimate of the G T matrix. At scan 40 the system is put into closed loop with, Ko, the inverse of G x as the pre-compensator. The limit on the bending moment of the slice was arbitrarily set to be ±0.23 Nm. The dynamics are simply a time delay of one scan and the closed loop pole, p, is set at 0.5. Figure 5.66 shows the inputs and outputs of the simulation both of the open loop using PRBS inputs and the closed loop with some step changes. It is observed that there are small overshoots for the parameters of 8th-10th order. Those overshoots arise when the controller needs to limit the slice lip bending. For this particular slice lip profile the above mentioned parameters need to be modified in order to limit the bending. Figure 5.67 shows the recursive estimate of the rows of the G x matrix, for the same period as the last figure. As mentioned above the estimated matrix is initially taken to be the identify matrix. Some of the parameters converge within 40 scans, when the system is in open loop. Other parameters do converge when the step changes are applied to the system at scans 60 and 90. Figure 5.68 shows the slice lip profile at scan 110 and the calculated bending moment. The figure also show the output and the desired output at the same scan. It can be seen in the figure that there is a difference between the two. The reason for the difference is that the bending moment of the slice lip is limited at ±0.23 Nm. If those limits are removed the difference between the output profile and the desired profile disappears, at the cost of bending moments increasing substantially as can be seen in Figure 5.69. It was mentioned in Chapter 3 that a conventional control might become unstable in the presence of uncertainty for a plant with a badly conditioned interaction matrix. The condition of the interaction matrix is of no concern if Gram polynomials are used, only the condition number of the parameterized interaction matrix. It was shown in Chapter 3 that the condition number of G x can be used as one of the selection criteria for the number of orthogonal polynomials used in the approximation of the profiles. As long as the condition number of the parameterized interaction matrix is low, the error term G xAGxl in the loop transfer matrix, c(l + GxAG~l) will be small. 5.2.1 Removin g polynomial s fro m a profil e Another simulation similar to the one above was performed, with the dry weight profile assumed to be one of the profiles from the industrial data set (see the solid line in Figure 5.70). The Dahlin controller with closed loop pole at 0.5, is then to remove all the low order polynomials from the profile (the dashed line in Figure 5.70) without exceeding the bending limits of the slice. The initial 91 Cliapter 5: Basis weight Cross Directional Control « - ' ' V ., > , / > • n^vi^_r^ 20 40 60 80 100 120 60 8 0 Scans 120 Figure 5.66: Input, u c, and output, y c. The system is initially in open loop but the loop is closed at scan 40. 92 Chapter 5: Basis weight Cross Directional Control Row 1 Row 2 0.5. 1 .5 0 * \ • • •^^TTttSif 'feE I 50 100 Scans Figure 5.67: Estimation of each row of G x. The true parameters for this simulation, are shown with dotted lines. The system is initially in open loop but the loop is closed at scan 40. 93 Chapter 5: Basis weight Cross Directional Control "i 1 1 r o a. c o S o e 60 g -3 e pa o 3 150 200 250 CD measurements 300 350 Figure 5.68 : Slice lip profile, bending moments and output at scan 100. o q> V3 Cu a o *—, 2 0 c E o e 00 c T5 c ID 03 I 0. 5 40 50 60 70 80 90 100 Actuator number o & "3 D. o 150 200 250 CD measurements Figure 5.69: Slice lip profile, bending moments and output at scan 100, when no limits are used for the slice lip. 94 Chapter 5: Basis weight Cross Directional Control 0 50 100 150 200 250 300 350 400 450 Measurement Figure 5.70: Dry weight measurements for an industrial data set, and its polynomial approximation. position of the slice lip is assumed to be the same as one of the slice lip profile from the same industrial data set (see the dot-dashed line in Figure 5.71). The system is initially kept in open-loop to estimate the parameterized interaction matrix, and the loop is closed at time sample 40. The initial and the final setting for the slice lip deflection, as well as their Gram polynomial approximations, are shown in Figure 5.71. It can be seen that only minor adjustments are needed for the slice lip position to remove the first 10 Gram polynomials from the profile as can be seen in Figure 5.72. Here the final output profile is shown with a solid line and its 10th order Gram polynomial approximation with a dashed line. That can be compared with Figure 5.70 to see how the low order polynomials have disappeared. The 2 sigma value of the final profile is 0.80 g/m2 down from 0.94 g/m2 for the initial profile. It is worth noting that a very little change is observed in the bending moments when the slice lip is manipulated to remove the 10 first polynomials from the output. Figure 5.73 shows the bending moments of the slice lip. The figure has actually two lines, one before control, and one after control, but the changes are so small that they cannot be distinguished on the graph. The bending moment changes at most 5Nm from the initial bending moment, or only about 2% of the maximum value of the bending moments. 95 Chapter 5: Basis weight Cross Directional Control T 1 — , r J i L _ 30 40 50 60 70 Actuator 90 100 Figure 5.71 : Slice lip actuators settings. 200 250 Measurement 450 Figure 5.72 : Final weight profile. 96 Chapter 5: Basis weight Cross Directional Control 50 60 Actuator Figure 5.73: Bending moments. 5.2.2 Contro l o f heav y weigh t pape r A paper machine with 132 actuators and 132 measurements is simulated assuming its response to the center actuator moving one unit is as shown in Figure 5.74. It can be seen that the response to a single actuator covers one fourth of the width of the machine. This is a similar response as suggested in [21] and [71]. This fictional paper machine is assumed to have an initial output profile that consists of a single sine wave across the width of the machine. This sinusoidal wave is supposed to be removed from the profile, but first the machine is run in open-loop to estimate the parameterized interaction matrix. The inputs are PRBS signals as can be seen in Figure 5.75. The outputs are also shown in the same figure. The estimates of the rows of G x are shown in Figure 5.76. It can be seen that there is a small bias in most of the parameters. The reason is the small amplitude in the input signal, which needs to be small in order to limit the stroke of the actuators. If the amplitude of the input signal is increased a better convergence is achieved. The final estimates are shown in Figure 5.77. Figure 5.78 shows the final basis weight output, the desired basis weight, and the initial basis weight profile. The slice lip deflection required to reduce the initial basis weight profile to the final basis weight profile, is shown in the same figure, as well as its bending moment. 97 Chapter 5: Basis weight Cross Directional Control 21 1 •— | 1 1 — r r 1.5-1 -i S 0.5-< "5b 0 --0.5 -] I i i i i i i I 0 20 40 60 80 100 120 Actuator Figure 5.74 : Response to a slice lip movement for a heavy grade paper. 5.3 Compariso n t o conventiona l contro l Figure 5.79 shows a block diagram of a very simple conventional controller utilizing the inverse of the interaction matrix to decouple the actuators for control. More elaborate scheme exists that take into account the control action as well as minimizing the output error [126]. For simplicity only the basic control scheme will be considered here, and the number of actuators are set to be equal to 250 and the number of measurements are four times the number of actuators, or 1000. The interaction matrix, G, has the same p m<s parameters as before (see Figure 4.45). Figure 5.80 shows a control of a profile consisting of 5 sinusoidal waves when a conventional control is used. The five sinusoidal waves have wavelengths that are 100, 30, 10, 4, and 1 cycles per sheet width. It can be seen that the five sinusoidal wave are completely eliminated from the profile except the high frequency one which has still got a small magnitude. Notice that the slice lip deflection is very jagged and consequently the bending moments are very high. When the process is controlled by utilizing Gram polynomials up to 40th order, it is observed that the slice lip is much smoother as can be seen in Figure 5.81. With the use of Gram polynomials in the control design, the lower frequencies have been eliminated but the higher frequencies, which would require high bending moments on the slice lip, are not altered. The lack of high frequencies 98 Chapter 5: Basis weight Cross Directional Control 0.04-r 0 . 0 2 -o 3 o--0.02- -? y y /• ^ ^ ^ ^ -? -?*-10 9 8 7 6 5 4 3 2 1 100 Scan Scan Order Figure 5.75 : Input and output for a heavy grade paper machine. 99 Cliapter 5: Basis weight Cross Directional Control Row 1 Row 2 15 10 5 0 -5 ^ - U • • '•' . 20 F Row 4 Row 6 Row 10 Figure 5.76 : Estimates of each row of G x for a heavy grade paper machine. 100 Chapter 5: Basis weight Cross Directional Control Gx Figure 5.77 : Final estimates of G x for a heavy grade paper machine. Output, desired output, and initial output S 20 \ 2 6 o £ -2 0 20 40 60 80 Actuator 100 120 Figure 5.78 : Final basis weight profile, the desired profile, and the initial profile. The slice lip deflection. The bending moments. 101 Chapter 5: Basis weight Cross Directional Control "j) — ^ i c u G"1 ; — i y - q G i Figure 5.79: Block diagram of a conventional control. £ 2 CO S 0 •- 9 -*MNV**V»^ vJftf *v^i*v<rt^^'W-^^»vv^-''*>M*wv.*^i-jy«vwAM<y^ w < t f " V ^ * 0 100 200 300 400 500 600 CD measurements 700 800 900 1000 T - 4 c o o (1) TD 10 0 -10 1 1 . , . 1 i i , , I , I . . i A A ^ V V V M A A * A A / \ A A I , ,H/iA/A i i i A U I " , I I i,i.i .,i , uKAm i LmAM'^WMW^AiAi^mmrlm^ iLLiM nifMAhi^hi t iMM0kmW^ ' ivMfH wye f-iWMWMlMn V'VMKtMhftk iilA/VV'i'i'f ' 'WW ! ' i 1 ' ' u p u | ' ' T I ' vWrWrnl ' I ' l l 50 100 150 Actuator 200 250 100 150 200 250 Actuator Figure 5.80: Conventional control of basis weight. The dotted line in the top graph denotes the initial profile, and the solid line denotes the final profile. deflections of the slice lip results in the inability of the controller to eliminate the high frequencies in the output. What is gained, is that the maximum value of the bending moment is only 156Nm, compared to 5750Nm for the conventional control. The condition number of the interaction matrix can cause problems for the control of the paper machine for certain type of disturbances, as was discussed in Chapter 3. Two examples are shown here, one for a badly condition matrix and the other for a better conditioned matrix. In the first a 105 x 105 interaction matrix with pn = hpi — 0.82366, and p2 — —0.41183 is used. Those nominal values are assumed to be known and its reciprocal condition number is 0.013. The diagonal of the interaction matrix has uncertainty of the form of uniform random numbers between —5 - 10~9 and 102 Chapter 5: Basis weight Cross Directional Control c o o <D "0 O O 00 1 10 0 -10 i £ 200 0 a 0 | -2000 £ -4000 100 200 300 400 500 600 700 800 900 1000 CD measurements 50 100 150 Actuator 100 150 Actuator 200 250 200 250 Figure 5.81 : Control with Gram polynomials. 10 Time Figure 5.82 : Disturbance rejections using conventional control on a well conditioned interaction matrix. 103 Chapter 5: Basis weight Cross Directional Control Figure 5.83: Disturbance rejections using conventional control on a badly conditioned interaction matrix. 5 • 10 - 9 . The Dahlin controller is then to remove a disturbance, or an initial profile, in the output, which is the same as the minimum singular vector of the interaction matrix. The results are shown in Figure 5.82. It is observed that the disturbance is eliminated with the controller. The second example uses an interaction matrix of the same size, but with the elements on first and the second diagonals slightly different. The values are: p1 = 0.784935028, and p2 = -0.392467514. The reciprocal condition number is now 1.43 • 10~10, or much closer to zero (much worse) than in the first example. The same uncertainties are used and the control results are shown in Figure 5.83. It is seen that the controller goes in the wrong direction initially but is able to remove most of the disturbances from the profile. It is worth noting in order to remove the disturbances the slice lip has to make deflection that are 100 million times bigger than in the first example, or in other words not realizable. For both examples, controller using Gram polynomials will not reduce the high frequency initial profiles. The slice lip will therefore not be altered and consequently the bending moments are zero. 5.4 Chaos From the diagonal shape of the parameterized interaction matrix, G x, which is apparent in this and the previous chapter, it is tempting to assume the matrix is a diagonal matrix with zero 104 Chapter 5: Basis weight Cross Directional Control entries outside its main diagonal. Then, only the elements on the diagonal of the Gx matrix need to be identified. Furthermore it is assumed that the interaction matrix, G, is described by only two parameters; a parameter on the main diagonal, and another parameter on the diagonal above and below the main diagonal. The parameterized interaction matrix, G x, can then be described by two parameters and estimated with k = 1 using Equation 3.62. Having obtained the G x matrix the pre-compensator matrix, KQ, is then used to calculate current control signal. To simplify the calculation of the pre-compensator matrix, it is further assumed that the G Xl matrix (see Equation 3.59) is also a diagonal matrix so the matrix KQ is calculated by Ko(i,i)= ^ \ ,. ., (5.130) Po+PiGXl(i,i) For simulation purposes and interaction matrix, G, is created for 105 actuators and 105 measurements with po = 1 and p\ = 0.4. The form of the Gx matrix is plotted in Figure 5.84 for 10th order polynomials used to approximate the measured basis weight profile and also for the control of the slice lip. The controller is chosen as a Dahlin controller, with the model of the process as ^ r Z J = 1, i.e. all the gain will be attributed to the interaction matrix, G. The delay of the simulated process will be assumed to be known and equal to one. The Dahlin controller is then -n i - P Mi' 1) cfo-1) = [l-pq-1-(l-p)q-d]B(q-i) ( 5 m) 1-p 1-q-1 and p, the closed loop pole, is used as the design parameter. The block diagram shown in Figure 5.64 is used. The paper machine is modeled by an interaction matrix that is block diagonal with uniformly distributed random numbers on the 3 main diagonals. It is maintained constant throughout. The desired profile is a sine wave and is approximated with 10 Gram polynomials of order one to 10 inclusively. The Gram polynomial approximates the sine wave perfectly as can be seen in the upper half in Figure 5.86, where the desired profile is shown as well as the difference, magnified 10000 times, between the profile and its approximation. The lower left graph in the figure shows the parameters of the orthogonal approximation of the desired profile. When this system is simulated with the aforementioned control and estimation method, the output profile is found to be the same as the desired profile and the parameters converge to their true values. Noise is then added to the entries on the diagonals in the interaction matrix. The noise is uniformly distributed random numbers between -0.2 and 0.2. The range of values for/>o are between 105 Chapter 5: Basis weight Cross Directional Control Figure 5.84: Meshed surface of the G x matrix for G of size 105 x 105 with po = 1 and p\ = 0.4. 2 1 1 1 1 1 1 1 1 1 r A I i i i i i i , i i i L _ 0 10 20 30 40 50 60 70 80 90 100 Experiment Figure 5.85 : Estimates of po and p\. 106 Chapter 5: Basis weigh! Cross Directional Control « -0.5 -2 n c ^ ^Equation 2.2: U 5 Order 10 Figure 5.86: The desired profile and 10000 times the difference between the profile and its approximation with 10th order Gram polynomials. The smaller graphs show the parameters of the approximation for two sets of the Gram polynomials. 0.8 and 1.2 and the range for pi is between 0.2 and 0.6. The Gx matrix for this particular interaction matrix is shown in Figure 5.87. It can be compared to Figure 5.84 and it is seen that their difference is very small. The process is simulated with this particular structure of the interaction matrix. The estimated po and pi are shown in Figure 5.85 for 100 interaction matrices with different noises. It can be seen that the estimated parameters are around the true values even though there is a small bias for each experiment. The output profile is very close to the desired profile. The same tests are done using different form of the Gram polynomials, i.e. the polynomials described by Equations 2.12 - 2.14. The parameters needed to approximate the sinusoidal profile are shown in the lower right graph in Figure 5.86. The G T matrix for G with dimension 105 x 105v pi = 1.0 and p2 = 0.4, is plotted in Figure 5.88. It can be noted that the first order polynomial has quite a large effect on the ninth order polynomial, because the ninth order polynomial is much smaller in amplitude than the first order. The same goes for the effect the second order polynomial has on the 10th order polynomial. For an interaction matrix with no noise the output profile is found to be the same as the desired profile, just like the case when the other form of the Gram polynomials was 107 Chapter 5: Basis weight Cross Directional Control Figure 5.87: Meshed surface of the Gx matrix for G of size 105 x 105 with pQ = 1, p\ = 0.4 and uniformly distributed random numbers between -0.2 and 0.2 added to po and p\. used. When the uniformly distributed noise is added to the interaction matrix the system is observed to go into chaos with on-line identification of the plant matrix. Figure 5.90 shows a bifurcation diagram (the output for CD point 10 as a function of the previous output for that point) for certain values of the closed loop pole, p. In the figure, pseudo phase plane plots for different values of p are also plotted. The fractal dimension can be estimated by calculating the correlation dimension for the seven phase plane plots in Figure 5.90, using an algorithm in [97]. It is found that the dimension is zero for periodic orbit, for the limit cycle it is equal to one and for p — 0.580 and p = 0.550 the dimension is 1.5, which suggests the system is chaotic. By zooming in on Figure 5.90 it can be observed that the output goes from being at a fixed value, when p — 0.64, to a limit cycle, when p = 0.60200. Then at p = 0.60199 the system has period three cycle, and when p is decreased further, the system has period 6, 12, 24 etc. cycle oscillation and then around p = 0.584501, chaos occurs. The values of p where the bifurcation occurs is shown in Table 5.3. The table also lists the ratio between the difference of the points of bifurcation. The ratio is very close to Feigenbaum's number (approximately 4.669). The small difference is mainly because of errors in finding the exact point of bifurcation. For example the point where the system goes from period 96 to period 192 appears around p = 0.58451000, but according to previous two values of point of bifurcation, and 108 Chapter 5: Basis weight Cross Directional Control the Feigenbaum's number, this should appear at p = 0.58450992 which is within the error bounds for determining the points of bifurcation. Figure 5.88: Meshed surface of the Gx matrix using the Gram polynomials from Figure 2.3. A 1% difference in initial condition on the desired value of the output profile does affect the response in such a way that it becomes uncorrelated to the nominal response after about 25 samples (see Figure 5.91), which further implies chaotic motion. It should be emphasized here that the chaotic behaviour only show up in the system when the badly scaled polynomials are used. 5.5 Summary In this chapter the Dahlin controller has been described, with and without integrator antiwindup. The Dahlin controller is then used together with Gram polynomial approximation of the measured profile to control the basis weight on a paper machine by means of adding Gram polynomials to the slice lip of the headbox. The control method is compared to a conventional method, which involves the inverse of the interaction matrix. The method using Gram polynomials is always found to require 109 Chapter 5: Basis weight Cross Directional Control Figure 5.89: Meshed surface of the Gx matrix for uniform G using the Gram polynomials from Figure 2.3. much less bending on the slice lip, as it only controls the low frequencies in the output profile. It is also not sensitive to the condition number of the interaction matrix. It is observed that when some approximations are made and a form of the Gram polynomials are used that has bad scaling, the control of the process can become chaotic. The chaotic behaviour has only been observed for the polynomials that are badly scaled. No instances of chaos are observed when Gram polynomials are used, which all have a comparable amplitude. 110 Chapter 5: Basis weight Cross Directional Control p = 0.64000 0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 P 1.5 ^ 1 + >^0.5 0 0 1 y(k) p = 0.60200 p = 0.60199 1.5 ^ 1 + ^ 0 . 5 0 0 1 y(k) p = 0.55000 p = 0.58000 0 0 1 y(k) p = 0.58500 l .J c: l + M ^0.5 \ + \ l.D P 1 + M >;o.5 p = 0.59000 1.5 + >^0.5 0 0 1 y(k) Figure 5.90 : A bifurcation diagram and pseudo phase plane plots for different values of p. I l l Chapter 5: Basis weight Cross Directional Control period 6 12 24 48 96 192 P 0.6013 0.5882 0.585360 0.5846875 0.58454125 0.58451000 Pk-Pk—i Pk+i—Pk 4.613 4.223 4.598 4.68 Table 5.3 The points of bifurcation and the ratio of the differences between points of bifurcation. 0 5 10 15 20 25 30 35 40 45 50 Samples Figure 5.91 : Response of the system with p = 0.55, for 1% difference in initial condition on the desired value of the CD value at location 10. 112 Chapter 6 Summary and Future work 6.1 Summary In this thesis it has been shown how the basis weight on a paper machine can be controlled by manipulating the slice lip on the headbox, taking into account limits on the bending moments of the slice lip. A model that relates slice lip movement to basis weight response has been developed. The main feature of the model is an interaction matrix that describes how a movement of a single actuator will effect a large section of the measured profile through bending of the slice lip which sets up spreading waves on the forming wire, that further spread out on the wire. Some of the problems in creating the interaction matrix are discussed and it is found that it is difficult to create a matrix that will not create ripples in the model output. Another more serious problem is that the interaction matrix can, for many paper machines, become singular which might create problems when controlling the basis weight. Orthogonal polynomials with equal weight given to every data point, called Gram polynomials, are then introduced as a means of approximating the profile and then the parameters of the polynomials are fed back to the controller. Several methods are given to create the Gram polynomials, both sequentially and recursive algorithms. A unified sequential formula is given which all the other methods can be written in terms of. It is observed that a bad selection of the scaling of the polynomials can lead to problems in control. For a particular selection of scaling factors, it was shown that the Gram polynomials approach the Legendre polynomials for number of data points much larger than the order of the polynomial. The Gram polynomials are therefore sometimes referred to as discrete Legendre orthogonal polynomials. Some methods are discussed for finding the optimum number of polynomials to approximate a profile. It is found that a modified AIC gives a good balance between complexity and fit. When choosing a polynomial order, it is also important to consider the condition number of the parameterized interaction matrix, if a robust performance is desired. The parameters of the Gram polynomial approximation of the slice lip can be related to the bending limit of a slice lip by modelling the slice lip as a flexible beam with multiple forces. The parameters can therefore be manipulated to limit the bending moments for the whole slice. 113 Chapter 6: Summary and Future work Compared to conventional control methods, the method using Gram polynomials is always found to require much less bending on the slice lip, as it only controls the low frequencies in the output profile. The method is also not sensitive to the condition number of the interaction matrix, which cannot be chosen at will. It does however depend on the condition number of the parameterized interaction matrix, which can be manipulated by adding or removing polynomial order. It is observed that when some approximations are made and a form of the Gram polynomials are used that has bad scaling, the control of the process will become chaotic. No instances of chaos are observed when Gram polynomials are used, which all have a comparable amplitude. 6.2 Futur e wor k As has been said in this thesis the mapping, between the actuators and the measurements, which is the key element of the current control systems, is effectively eliminated with the control system proposed in this thesis. The problem with shrinkage of the paper sheet has not been addressed in this work. If there is for example different shrinkage at the edges than the middle of the cross direction a different set of polynomials might be required. For a paper machine with multiple actuators controlling the same CD profile orthogonal poly-nomials might be used in classifying the measured profile such that one set of actuators will control certain set of polynomials and another set of actuators control another set of polynomials. 114 References [I] M. F. Aburdene. On the computation of discrete Legendre polynomial coefficients. Multidi-mensional systems and signal processing, 4(2):181—186, Apr. 1993. [2] W. L. Adams. New measurements provide new control and information analysis for high-speed machines producing newsprint and groundwood papers. In 76th Annual Meeting technical section, CPPA, volume 1/2, pages A289-A294, Montreal, Jan. 30-31 1990. [3] H. Akaike. Modern development of statistical methods. In P. Eykhoff, editor, Trends and progress in system identification, chapter 6, pages 169-184. Pergamon Press, 1981. [4] J. R. Amyot and S. Nuyan. A troubleshooting expert-system consultant for cross-papermachine dryweight control systems. In 75th Annual Meeting technical section, CPPA, pages A97-A100, Montreal, Jan. 31 - Feb 01 1989. [5] K. J. Astrom. Computer control of a paper machine - an application of linear stochastic control theory. IBM Journal, pages 389-405, July 1967. [6] K. J. Astrom and B. Wittenmark. Computer controlled systems. Prentice-Hall, second edition, 1990. [7] J. V. Beck and K. J. Arnold. Parameter estimation in engineering and science. John Wiley & Sons, 1977. [8] A. E. Beecher and R. A. Bareiss. Theory and practice of automatic control of basis weight profiles. Tappi, 53(5):847-852, May 1970. [9] L. G. Bergh and J. F. MacGregor. Spatial control of sheet and film forming processes. The Canadian journal of chemical engineering, 65, Feb. 1987. [10] C. Bermond and J. Veyre. Le systeme ctp de regulation du profil de grammage. Revue A.T.I.P., 43(5):223-226, May 1989. [II] W. H. Beyer, editor. CRC Standard mahtematical Tables. CRC Press, 27th edition, 1984. p. 372. [12] D. Biros. Koehler keeps control of quality. PPI, 31(12):54-55, Dec. 1989. [13] T. J. Boyle. Control of cross-direction variations in web forming machines. The Canadian journal of chemical engineering, 55:457-461, Aug. 1977. 115 [14] T. J. Boyle. Practical algorithms for cross-directional control. TAPP1, 61(l):77-80, Jan. 1978. [15] T. J. Boyle and H. Hamby. Cross-direction weight control: pro and con; A tale of two machines. Tappi, 68(6):56-59, June 1985. [16] G. Burkhard and P. E. Wrist. The evaluation of paper machine stock systems by basis weight analysis. Tappi, 37(12), Dec. 1954. [17] E. W. Carey, C. R. Bietry, and H. W. Stoll. Performance factors associated with profile control basis weight on a paper machine. Tappi, 58(6):75-78, June 1975. [18] B. Carnahan, H. A. Luther, and J. O. Wilkes. Applied numerical methods. John Wiley & Sons, 1969. [19] T. Cegrell and T. Hedqvist. Successful adaptive control of paper machines. Automatica, 11:53-59, 1975. [20] S.-C. Chen. Kalman filtering applied to sheet measurement. In American Control Conference, volume 1 of 3, pages 643-647, Atlanta, Georgia, June 15-17 1988. [21] S.-C. Chen. Cross machine profile control for heavy weight paper. In EUCEPA, pages 77-94, Stockholm, May 8-11 1990. [22] S.-C. Chen. Acuation cell response and mapping determinations of web forming mahcines, 1991. Patent. [23] S.-C. Chen. Full-width sheet property estimation from scanning measurements. In Control Systems '92, pages 123-130, Whistler, B.C., Sept. 29 - Oct. 1 1992. [24] S.-C. Chen, R. M. Snyder, and R. G. Wilhelm. Adaptive profile control for sheetmaking processes. In PRP6, Oct. 1986. [25] S.-C. Chen and R. G. Wilhelm. Optimal control of cross-machine direction web profile with constraints on the control effort. In American Control Conference, volume 3, pages 1409-1415, June 18-20 1986. [26] F. Church. Transverse profile control of basis weight/moisture. Pulp & Paper, pages 48-50, Apr. 1973. [27] S. H. Crandall, N. C. Dahl, and T. J. Lardner. An introduction to the Mechanics of Solids. McGraw-Hill, second edition, 1972. [28] W. H. Cuffey. Some factors involved in basis weight uniformity. TAPPI, 40(6):190A-196A, June 1957. 116 [29] K. Cutshall. The nature of paper variation. Tappi Journal, pages 81-90, June 1990. [30] K. Cutshall. Cross-direction control. In B. A. Thorp and M. J. Kocurek, editors, Paper Machine Operations, volume 7 of Pulp and Paper Manufacture, chapter XVIJI, pages 472-506. Tappi and CPPA, 1991. [31] K. A. Cutshall and R. G. Hamel. Slice lip automation. Pulp & Paper Canada, 89(l):85-92 (T1-T7), Jan. 1988. [32] K. A. Cutshall, G. E. Ilott, and J. H. Rogers. Grammage variation-Measurement and analysis. Grammage Variation subcommittee, process control committee, technical section, CPPA, 1988. [33] E. B. Dahlin. Computational methods of a dedicated computer system for measurement and control on paper machines. In TAPPI 24th Engineering conference, San Francisco, California, Sept. 14-19 1969. [34] J. R. Dukes. Here are four steps to superior basis weight measurement and control. PIMA, 70(8):41-43, Aug. 1988. [35] G. A. Dumont. Application of advanced control methods in the pulp and paper industry-a survey. Automatica, 22(2): 143-153, 1986. [36] G. A. Dumont. On the use of adaptive control in the process industries. In CPC-III, pages 467-499, 1986. [37] G. A. Dumont. System identification and adaptive control in papermaking. In C. F. Baker and V. W. Punton, editors, Fundamentals of papermaking. Transactions of the ninth fundamental research symposium, volume 2, pages 1151-1181. Mechanical Engineering Publications Limited, Cambridge, Sept. 1989. [38] G. A. Dumont. Control techniques in the pulp and paper industry, volume 37 of Control and dynamic systems, pages 65-114. Academic Press, 1990. [39] G. A. Dumont, M. S. Davies, K. Natarajan, C. Lindeborg, F. Ordubadi, Y. Fu, K. Kristinsson, and I. M. J6nsson. An improved algorithm for estimating paper machine moisture profiles using scanned data. In 30th. IEEE CDC, Brighton, England, Dec. 1991. [40] G. A. Dumont, I. M. Jonsson, M. S. Davies, F. T. Ordubadi, Y. Fu, K. Natarajan, C. Lindeborg, and E. M. Heaven. Estimation of moisture variations on paper machines. IEEE Transaction on Control Systems Technology, 1(2):101-113, June 1993. 117 [41] K. Dutton. An investigation into the design and performance of an automatic shape control system for a Sendzimir cold rolling mill. PhD thesis, Sheffield City Polytechnic, May 1983. [42] K. Dutton, J. F. Barrett, and M. J. Grimble. The shape control problem for a Sendzimir cold steel rolling mill. In Symposium on application of multivariable systems theory, pages 217-227, Oct. 1982. [43] K. Dutton and R. J. Cleland. Shape control system for Sendzimir mills. In Advances in cold rolling technology, Sept. 17-19 1985. [44] Editor. New system controls CMD basis weight profile automatically. Paper Trade Journal, pages 34-36, Nov. 1971. [45] R. Elliot. Modeling the slice lip profile for basis weight control. PIMA, pages 32—36, Feb. 1984. [46] R. M. Elliot. Model-based cross directional control. In Proc. World Pulp and paper week, Apr. 10-13 1984. [47] J. Enzi. Querprofilregelungen an einer Tiefdruckpapiermachine. Wochenhlatt fur Papierfab-rikation, 116(22):958-960, 1988. [48] S. Evans, S. Nuyan, and R. Amyot. Expert system diagnosis in CD weight control. In 76th Annual Meeting technical section CPPA, volume 1/2, pages A231-A235, Montreal, Jan. 30-31 1990. [49] R. A. Fisher and F. Yates. Statistical tables for biological, agricultural and medical research. Hafner Publisher Company, 1963. p.98. [50] J. Flading. Cross-machine control. In Tappi Wet end operations, pages 319-321, 1988. [51] G. E. Forsythe. Generation and use of orthogonal polynomials for data fitting with a digital computer. SI AM Journal, 5(2):74-88, June 1957. [52] B. C. Garnett. Gaining and edge: CD profile control at the headbox. PIMA, pages 19-23. Feb. 1984. [53] K. Gibson. Profile control system at Weyco mill cuts sheet weight variations. Pulp & Paper, 58:68-71, Feb. 1984. [54] P. A. Gorry. General least-squares smoothing and differentiation by the convolution (savitzky-golay) method. Anal. Chem., 62:570-573, 1990. [55] A. Grace. Optimization toolbox for use with MATLAB™'. User's Guide. The MathWorks. Inc., Cochituate Place, 24 Prime Park Way, Natick, Mass. 01760, Nov. 1990. 118 [56] M. J. Grimble and J. Fotakis. The design of shape control systems for Sendzimir mills. IEEE Transactions on Automatic Control, 27:656-666, 1982. [57] R. Hall. Coordinated cross-machine control. Paperi ja Puu — Paper and Timber, 73(3):241-248, 1991. [58] A. Halouskova, M. Karny, and I. Nagy. Adaptive cross-direction control of paper basis weight. Automatica, 29(2):425-^129, Mar. 1993. [59] A. Hansson, L. Haglund, and I. Gladh. Experience with CD-basis weight profile control. EUCEPA, 1982. [60] D. E. Harke. Remote manual profile control on a three-ply liner machine — a mill experience. In Papermakers Conference, pages 197-200. Tappi, Apr. 6-8 1987. [61] E. Hauptmann, R. Vyse, and J. Mardon. The wake effect as applied to modern hydraulic headboxes part I:. Pulp and Paper Canada, 91(9):T357-T364, Sept. 1990. [62] E. Hauptmann, R. Vyse, and J. Mardon. The wake effect as applied to modern hydraulic headboxes part II:. Pulp and Paper Canada, 91(10):T369-T377, Oct. 1990. [63] E. M. Heaven, I. M. Jonsson, T. M. Kean, M. A. Manness, and R. N. Vyse. Recent advances in cross machine profile control using modern estimation and control technology. In Proceedings of the second IEEE conference on control applications, volume 1 of 2, pages 217-226, Vancouver, B.C., Sept. 13-16 1993. [64] E. M. Heaven, T. M. Kean, I. M. Jonsson, K. M. Vu, M. A. Manness, and R. N. Vyse. Applications of system identification to paper machine model development and controller design. In Proceedings of the second IEEE conference on control applications, volume 1 of 2, pages 227-234, Vancouver, B.C., Sept. 13-16 1993. [65] E. M. Heaven, M. A. Manness, K. M. Vu, and R. N. Vyse. Applications of system identification to paper machine model development and simulation. In Control Systems '92, pages 271-278, Whistler, B.C., Sept. 29 - Oct. 01 1992. [66] F. B. Hildebrand. Introduction to numerical analysis. McGraw-Hill, 1956. [67] D. M. Himmelblau. Process analysis by statistical methods. John Wiley & Sons, 1970. [68] D. Hunter. Papermaking. The history and technique of an ancient craft. Dover Publications, 1947. 119 [69] D. F. Jenkins. Cross-machine variability in coated publication paper production. Tappi, 56(2):89-94, Feb. 1973. [70] f. M. J6nsson. Estimation and identification of moisture content in paper. Master's thesis, University of British Columbia, 1991. [71] H. Karlsson and L. Haglund. Optimal cross-direction basis weight and moisture profile control on paper machines. In The proceedings of the third international pulp and paper process control symposium, pages 139-145, Vancouver, B.C. Canada, May 2-4 1983. [72] H. Karlsson, L. Haglund, and A. Hansson. Optimal cross-direction basis weight and moisture profile control on paper machines. Pulp and Paper Canada, 86(8):T213-T218, Aug. 1985. [73] H. Karlsson, I. Lundqvist, and T. Ostman. Principles and potentials of CD-basis weight control. EUCEPA, pages 238-243, 1982. [74] G. Kastanakis. Prediction and verifying the best quality control under particular conditions. Tappi Journal, 73(12):203-207, Dec. 1990. [75] M. A. Keyes. Cross machine moisture control system for the wet end of a paper machine. Patent, Beloit Corporation, 1971. U.S. Pat. No. 3.619.359. [76] M. A. Keyes. Computer control census. Tappi, 58(6):71-74, June 1975. [77] K. Kristinsson and G. A. Dumont. Paper machine cross directional basis weight control using gram polynomials. In IEEE Conference on Control Applications, volume 1 of 2, pages 235-240, Vancouver, B.C., Sept. 13-16 1993. [78] K. S. Kunz. Numerical Analysis. McGraw-Hill, 1957. [79] D. L. Laughlin. Control system design for robust performance despite model parameter uncertainties: Application to cross-directional response control in paper manufacturing. PhD thesis, California Institute of Technology, 1988. [80] D. L. Laughlin, K. G. Jordan, and M. Morari. Internal model control and precess uncertainty: mapping uncertainty regions for siso controller design. International Journal of Control, 44(6): 1675-1698, 1986. [81] D. L. Laughlin and M. Morari. Robust performance of cross-directional basis-weight control in paper manufacturing. In American Control Conference, volume 3, pages 2122-2128, June 21-23 1989. 120 [82] J. G. Lebrasseur. Technical advances to improve linerboard moisture. Tappi, 67(10):60-61, Oct. 1984. [83] I. J. Leontaritis and S. A. Billings. Model selection and validations methods for non-linear systems. International Journal of Control, 45(1 ):311-341, 1987. [84] C. Lindeborg. A nonlinear algorithm for the estimation of moisture variation characteristics in paper process. In Control of industrial processes, volume 3, pages 139-158, Milan, Italy, Oct. 6-7 1981. [85] C. Lindeborg. A process model of moisture variations. Pulp and Paper Canada, 87(4):T142-T147, Apr. 1986. [86] L. Ljung. System identification: Theory for the user. Prentice Hall, 1987. [87] L. Ljung and T. Soderstrom. Theory and practice of recursive identification. MIT press, 1983. [88] M. Matsuda. A new decoupling matrix method for paper basis weight profile control. In Proceedings of the 29th conference on decision and control, volume 3/6, pages 1592-1594, Honolulu, Hawaii, Dec. 5-7 1990. IEEE. [89] D. A. McFarlin. Control of cross-machine sheet properties on paper machine. In The proceedings of the third international Pulp and Paper Process control symposium, pages 49-54, Vancouver, B.C., May 2-4 1983. [90] W. E. Milne. Numerical Calculus. Princeton University Press, 1949. [91] R. Munch. Der JETCOmmander-ein fortschrittliches System zur Fernverstellun der Stoffau-flaufblende und zur Regelung des Flachengewicht-Querproflls. Wochenblatt fur Papierfab-rikation, 120(7): 259-265, 1992. [92] K. Natarajan, G. A. Dumont, and M. S. Davies. An algorithm for estimation cross and machine direction moisture profiles. Post-Graduate research laboratory report PGRL 388, Pulp and Paper Research Institute of Canada, Oct. 1987. [93] C. P. Neuman and D. I. Schonbach. Discrete (Legendre) orthogonal polynomials - A survey. International Journal for Numerical Methods in Engineering, 8:743-770, 1974. [94] B. Noble. Applied linear algebra. Prentice-Hall, 1969. [95] L.-O. Norberg. Experience from CD control of basis weight on the fine paper machines in the MoDo, Husum, mill. Pulp & Paper Canada, 90(5):83-87, 1989. 121 [96] S. Nuyan. Requirements for the cross-machine direction basis weight control of the paper web. ACC, 3:1416-1422, June 18-20 1986. [97] T. S. Parker and L. O. Chua. Practical Numerical Algorithms for Chaotic Systems. Springer-Verlag, 1989. [98] R. S. Peterson. Improving basis weight uniformity with deckle wave control. Tappi Journal, pages 121-128, July 1992. [99] A. Ralston and P. Rabinowitz. A first course in numerical analysis. International series in pure and applied mathematics. McGraw-Hill, 1978. [100] G. A. Richards. Cross direction weight control. Japan Pulp and Paper, 537:41^15,53, Nov. 1982. [101] J. V. Ringwood. The design of shape control systems for a Sendzimir mill. PhD thesis, University of Strathclyde, 1984. [102] J. V. Ringwood and M. J. Grimble. Shape control in Sendzimir mills using both crown and intermediate roll actuators. IEEE Transactions on Automatic Control, 35(4):453-459, Apr. 1990. [103] R. C. Roth. Introductory remarks CD-control of moisture profiles. In Control Systems, pages 72-74, Stockholm, May 1986. [104] T. Sasaki, T. Otani, and M. Negishi. Paper machine controller for operation slices and method of controlling the same. Patent Number: 5,170,357, Dec. 8 1992. [105] R. S. Sheel and M. G. Gichard. No. 1 p.m. newsprint modernization program at domtar red rock upgrading qualtiy on an older newsprint machine. In Proceedings 71st annual meeting, technical section CPPA, pages A237-A243, Montreal, Quebec, Jan. 27-Feb 1 1985. [106] S. J. Siler. Cross machine basis weight control machine considerations. In Engineering conference, volume 3, pages 631-636, Boston, MA, 1984. Tappi. [107] S. Skogstad, M. Morari, and J. C. Doyle. Robust control of ill-conditioned plants: High-purity distillation. IEEE Transaction on automatic control, 33(12): 1092-1105, Dec. 1988. [108] K. E. Smith. Cross-direction control is still top process automation trend. Pulp & paper, 59(2):72-76, Feb. 1985. [109] K. E. Smith. Cross-direction control sales have doubles in the past two years. Pulp & paper, 61(2):48-51, Feb. 1987. 122 [110] K. E. Smith. Survey shows more than 3000 cross-direction controls now installed. Pulp & Paper, 63(2):57-60, Feb. 1989. [ I l l ] G. A. Smook. Handbook for pulp & paper technologists. Tappi Press, 1982. [112] G. A. Smook. Handbook of pulp and paper terminology; A guidebook to industrial and technological usage. Angus Wilde Publications, 1990. [113] T. Soderstrom and P. Stoica. Instrumental variable methods for system identification. Lecture notes in control and information sciences. Springer-verlag, 1983. [114] D. A. Spitz. Control of variations in a moving web with a scanning gauge. Tappi, 73(1): 133-137, Jan. 1990. [115] B. F. Taylor. Optimum separation of machine-direction and cross-direction product variations. Tappi, pages 87-92, Feb. 1991. [116] R. M. Tong. Automatic control of grammage profile on a paper machine. In 3rd Internationl IFAC Conference on instrumentation and automation in the paper, rubber and plastics industries, pages 289-298, Brussels, May 24-26 1976. IFAC. [117] D. Wahren. Recent highlights in paper technology. Tappi, 69(3):36-45, 1986. [118] B. Wallace. The challenge of higher-frequency MD and CD control. PIMA, 72(2):56-57, Feb. 1990. [119] B. W. Wallace. Economic benefits offered by computerized profile control of weight, moisture and caliper. Tappi, 64(7):79-83, July 1981. [120] B. W. Wallace, R. Balakrishnan, and M. Rodman. Advances in on-line measurement for tighter control of paper quality. Appita, 45(2):74-77, Mar. 1992. [121] M. D. Wallace. Cross-direction control: profitability and pitfalls. Tappi, 69(6):72-75, 1986. [122] X. G. Wang, G. A. Dumont, and M. S. Davies. Estimation in paper machine control. IEEE control systems, 13(4): 34-43, Aug. 1993. [123] P. E. Wellstead and W. P. Heath. Two-dimensional control systems: application to the CD and MD control problem. In Control Systems '92, 1992. [124] R. G. Wilhelm. Self tuning control strategies: Multi-faceted solutions to paper machine control problems. In Discover the new world of instrumentation, proceedings of the 1982 joint symposium, pages 207-215, Columbus, Ohio, Apr. 13-15 1982. Instrument Society of America. 123 [125] R. G. Wilhelm. On the controllability of cross-direction variations of sheet properties. In Engineering Conference, volume 3, pages 621-629, Boston, MA, 1984. [126] R. G. Wilhelm and M. Fjeld. Control algorithms for cross directional control: the state of the art. In Proceedings of the 5th 1FAC conference on Instrumentation and Automation in the Paper, Rubber, Plastics and Polymerization Industries, pages 139-150 or 163-174. Pergamon Press, 11-13 Oct. 1983. [127] A. J. Wilkinson and A. Hering. A new control technique for cross machine control of basis weight. In Proceedings of the 5th IFAC conference on Instrumentation and Automation in the Paper, Rubber, Plastics and Polymerization Industries, pages 175-182. Pergamon Press, 11-13 Oct. 1983. [128] F. Xingren and G. Lawrence. Less variation makes all the difference. PPI, 34(6):52-53, June 1992. [129] S. Yamamoto. Industrial developments in intelligent and adaptive control. In International Symposium. IFAC, Feb. 15-17 1991. Preprints session four, intelligent tuning and adaptive control. 124
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Cross directional control of basis weight on paper...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Cross directional control of basis weight on paper machines using Gram polynomials Kristinsson, Kristinn 1994
pdf
Page Metadata
Item Metadata
Title | Cross directional control of basis weight on paper machines using Gram polynomials |
Creator |
Kristinsson, Kristinn |
Date Issued | 1994 |
Description | In this thesis a control algorithm for cross-directional basis weight control based on orthogonal polynomials is presented. Instead of controlling the cross-directional values at each point, the parameters of the orthogonal polynomials are controlled, taking into account the maximum bending moment of the slice lip and the relative stiffness of the slice with respect to the actuators. Problems with slice lip bending are dealt with in two ways. Implicitly, by modelling it as a polynomial such that the slice is not bent more than the highest degree of the polynomial, which is chosen beforehand. Explicitly, by calculating the bending moment and then use quadratic programming to limit the bending moment. By modeling the measurements and the actuators positions by polynomials the mapping between the actuators and the sensors, which is the key element of the current control systems, is effectively eliminated. Computational requirements are greatly reduced by avoiding inverses of big matrices. |
Extent | 6693958 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-08 |
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.0065251 |
URI | http://hdl.handle.net/2429/6963 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-894202.pdf [ 6.38MB ]
- Metadata
- JSON: 831-1.0065251.json
- JSON-LD: 831-1.0065251-ld.json
- RDF/XML (Pretty): 831-1.0065251-rdf.xml
- RDF/JSON: 831-1.0065251-rdf.json
- Turtle: 831-1.0065251-turtle.txt
- N-Triples: 831-1.0065251-rdf-ntriples.txt
- Original Record: 831-1.0065251-source.json
- Full Text
- 831-1.0065251-fulltext.txt
- Citation
- 831-1.0065251.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065251/manifest