Estimation and Identification for MachineDirection Control of Basis Weight and Moistureby Scott Taylor MorganB.A.Sc. University of British Columbia, Vancouver, 1988.A THESIS SUBMIT lED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF ELECTRICAL ENGINEERINGWe accept this thesis as conformingto the required standardTHE UNWERSITY OF BRiTISH COLUMBIAJune 1994© Scott Morgan, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)______________________________Department of zEZ IA/C—The University of British ColumbiaVancouver, CanadaDate 25) j9’DE-6 (2/88)AbstractFor the regulation of paper properties on a paper machine, measurements taken by a gauge tracinga zigzag pattern on the sheet must be decomposed into cross direction (CD) and machine direction(MD) signals. Researchers at the Pulp and Paper Centre, UBC, Vancouver, Canada have developeda decomposition algorithm (EIBMC) shown to give improved CD and MD estimates over previousmethods. In this thesis the potential of the improved MD estimates for MD control is examined.Traditional scan-average MD estimates provide poor anti-aliasing and a sample rate fixed by theperiod of the scanning gauge. It is shown that the algorithm can be modified to provide estimatesat any submultiple of the sensor sample rate with improved filtering flexibility. It is shown thatCD estimation errors cause MD estimation errors proportional to the gradient of the CD errors. Ascheme for reducing the CD and MD estimation errors due to poor filtering using a method involvingsensor pre-filter design and manipulation of data in the EIBMC algorithm is presented. A computersimulation of a scanning sensor has been created and is used to illustrate the performance of thealgorithm and anti-aliasing filter. A method is proposed for calculating the optimum sample rate forMD control based on the continuous-time average variance of the process output given the processstatistics and time-delay. The method differs from previous methods in that it considers the effect ofthe anti-aliasing filter in obscuring the true output variance.Table of ContentsAbstractList of TablesList of FiguresAcknowledgments1 Introduction1.1 Background1.1.1 Paper Making1.1.2 Decomposition of Paper Properties1.1.3 Scanning Sensor1.1.4 Standard Estimation Methods1.1.5 EIBMC Algorithm1.2 Description of the Problem1.3 Organization of Thesis2 Process Models2.1 Moisture2.1.1 Remarks2.1.2 Process DynamicsBasis WeightEIBMC AlgorithmMeasurement Indexing and the ScannerThe EI]3MC Algorithm for Moisture3.2.1 Estimation of B3.2.2 Modification of the Mean Deviation from Target3.2.3 EFRA33 1UviVIIix111223557• . . . 10111214141417181892.23 The3.13.23.2.4 Kalman Filter.3.3 The EIBMC Algorithm for Basis Weight3.4 Parameter Selection4 Machine Direction Estimation Using the EIBMC Algorithm4.1 Commercial Scanners4.1.1 Partial Scan Reporting4.1.2 Sensing4.1.3 The Filtering Effect of the Sensing Method4.2 Direct Measurements4.3 Scanner Simulator4.4 Data Shifting in the EIBMC Algorithm4.5 The Effect of Pre-Filtering on Spatial Smearing4.5.1 CD Distortion4.5.2 The Effects of CD Estimation Error on MD Estimation4.6 Compensation for Pre-Filter Distortions4.6.1 Maximally-Linear-Delay Filter4.6.2 Bessel Filter Design Example4.7 MD Estimation4.8 Summary5 Continuous-Time Variance Analysis of Sampled Time-Delayed Systems5.1 Criteria Selection5.2 Overview of the Variance Calculations5.3 Continuous-Time Model5.3.1 Plant and Disturbance5.3.2 Anti-Aliasing Filter 465.3.3 Combined State Space Description for Process and Anti-aliasing Filter 46iv• . 18• . 19192122222224242527272933343435373941424344455.4 Discrete-Time Model of Process and Anti-Aliasing Filter 475.5 Estimation and Prediction 485.6 Controller Design 505.6.1 Sampling the Continuous-Time Cost Function 505.6.2 Control Law 515.7 Variance Calculations 515.8 Time-Average Variance 525.9 The Calculation of Matrix Exponential Integrals 535.10 Process and Measurement Variance for the Case of No Time Delay 565.11 Some remarks on the difficulties encountered with the variance calculations 585.12 Summary 596 Conclusions 606.1 Summary of Results and Contributions6.2 Suggestions for Further Research 61References 62VList of Thbles4.1 Bessel filter order factors 374.2 Bessel filter attenuation at half sample rate with the ratio of Nyquist frequency to cut-offfrequency, in brackets 37VIList of Figures1.1 Paper machine 11.2 Scanner trajectory1.3 Nonnalized Boxcar filter Bode plot 41.4 The variable filtering of the Partial Scan Estimation method 52.5 Process indexing according to Lindeborg’s model 103.6 Indexing for scanned data 153.7 MD estimate response to a step input, showing measurement and prediction updates . . . . 164.8 Sensor-sheet orientation and measurement turn 234.9 254.10 SystemBuild scanning sensor simulator 284.11 Data shifting scheme for forward sensor traverses 294.12 Sixth order Butterworth and Bessel filter Bode Plots with time delay expressed indataboxes 304.13 Spatial shifting of raw measurements due to filter delay 314.14 Detail of the true CD and raw measurement with mean removed 314.15 CD response using sixth order Bessel filter with and without data shifting 324.16 CD response using sixth order Bessel and Butterworth filters with data shifting 324.17 MD response using sixth order Bessel filter with and without data shifting 334.18 First 50 scans of MD response using sixth order Bessel filter with data shifting 344.19 MD ifitered and re-sampled every five seconds 39vu:5.20 Controller structure.435.21 Mathematical Model of Plant and Controller 455.22 Process and measurement variance for a non-delayed plant 57VI 11AcknowledgmentsI thank my thesis supervisor Professor G.A. Dumont for his kind support and valued guidance. Ithank Professor M.S. Davies for his interest and input. I thank Fariborz Ordubadi for his friendshipand for assistance converting the EIBMC algorithm into C-code. I thank Mike Yuswack of MeasurexInc. for kindly providing information on current scanner technology. I thank Bengt Lennartson ofChalmers University of Technology for providing all of his pertinent papers on variance analysis. Ithank the friends I have made at the Pulp and Paper Centre and those who have helped along theway, in particular, Yamei Chen, Ashraf El-Naggar, Ivar Jonsson, Kristinn Kristinsson, Rick Morrison,Ye Fu, and George Wang. I especially thank Jahan Ghofraniha for his warm friendship and manyhours of engaging conversation.The highest praise and gratitude goes to my loving partner in life, Mo, for bravely weatheringmy academic storms as well as her own.IxChapter 1: IntroductionChapter 1IntroductionThe subject of this thesis is the use of a new estimation algorithm for improved control ofmachine direction variations of basis weight and moisture content in the making of paper. The nextsection provides background on papermaking and estimation to be followed by a description of theproblems to be dealt with and how they will be approached. The chapter concludes with an overviewof the thesis.1.1 Background1.1.1 Paper MakingIn the manufacture of paper, a suspension of fibres and water is deposited onto a moving wiremesh known as the web or wire and, by means of gravity, vacuum and heat, the water is removed tothe extent that a cohesive paper sheet remains. Most of the drying heat is applied by steam heatedrolls over which the sheet passes in a part of the paper machine known as the dryer section. Thepaper is then rolled onto a reel for storage and transportation (See Figure 1.1).An automatic control system is typically used to regulate paper properties such as basis weight,which refers to the mass of fibres per unit area of the sheet, and moisture content, the percentagemass of water versus water and fibre. Basis weight is commonly controlled by thick stock flow tosteamthick stockthin stockcalenderheadbox dryer scanningwire press sensorFigure 1.1: Paper machinereel1Chapter 1: Introductionthe head box which alters the rate at which fibres are delivered to the wire. Wet basis weight refersto the total mass of fibres and water per unit area, and dry basis weight refers to just the mass ofthe fibres per unit area. Unless otherwise stated, the use of basis weight in this thesis will referto wet basis weight. Moisture content is typically controlled by the dryer section steam pressure.Measurements used for control of paper properties are usually made at a point after the dryer sectionand just before the sheet is rolled.1.1.2 Decomposition of Paper PropertiesVariations in paper properties are considered to be composed of two parts: the machine directionpart (MD) and the cross direction (CD) part, where machine direction means in the direction of themovement of the paper sheet and cross direction means perpendicular to machine direction and inthe plane of the paper sheet. The MD portion of the property refers to variations that occur along thelength of the sheet and, for any point along the length of the sheet, is defined as the average valueof that property taken across the sheet. For a given position along the length of the sheet, the CDportion of the property refers to the set of values obtained by subtracting the MD value for that MDposition from the values of the property at points along the width of the sheet. For any given MDposition the set of values that constitute the CD part will collectively be known as the profile. Notethat by definition the profile values sum to zero.1.1.3 Scanning SensorMeasurements used in the regulation and control of the paper properties are made by a sensor,mounted on an 0-frame surrounding the moving sheet, that continually traverses back and forth inthe CD direction as the sheet rushes by. A traverse from one side of the sheet to the other is referredto as a scan. Samples are taken at uniform positions along each scan. At the end of each traverse,the sensor pauses for a few seconds before reversing. The number of samples taken each traverse canbe as little as 20 or as much as 3000 and scan periods can vary from ten to sixty seconds dependingon the installation. Since the paper sheet is in motion in a direction perpendicular to the path of thesensor, the sensor traces out a zigzag pattern over the sheet (See Figure 1.2). The paper travels much2Chapter 1: IniroductionMachine Direction Paper sheet4— Sensor off-sheet4Cross DirectionPath of sensorrelative to sheetFigure 1.2: Scanner trajectoryfaster than the sensor and so the trajectory of the sensor relative to the paper makes a very smallangle with the machine direction. An angle of 1° is not uncommon (see Dumont[4]). The estimatesof MD and CD quantities required for control must be derived from data available from the sensor.However, because the sensor moves in the zigzag fashion relative to the direction of the paper, thesensor data contains a mixture of MD and CD information and so the process of estimating MD andCD quantities is not straightforward.1.1.4 Standard Estimation MethodsIn most mills, the MD estimates are made by averaging the samples taken during each traverseof the sensor. This method is called scan average estimation. It is assumed that the CD profileschange slowly (see Dumont[4]). Hence, the CD variations, which by definition sum to zero for agiven MD position, will tend to cancel out in the sum of the samples in a traverse, leaving onlythe effects of MD variations. However, this estimation method imposes two important limitations.Firstly, based on the assumption that the CD component is constant in time, the method results inan averaging of the true MD during the traverse. This averaging is equivalent to using a boxcar orrectangular continuous-time filter before sampling and may not provide the required pre-filtering asdetermined by the frequency distribution of disturbances for the MD properties of interest. Figure3Chapter 1: Introduction10’frequency (Hz)Figure 1.3: Normalized Boxcar filter Bode plot1.3 shows a normalized Bode plot for boxcar filter. The sample rate is one second. The frequencybeyond which aliasing due to sampling occurs, 0.5 Hertz, is indicated with a vertical dashed line.It can be seen that pre-fikering is poor. At half the sample rate the signal has been attenuated byonly 4dB. For scans on a paper machine this frequency corresponds spatially to twice the width ofthe sheet. As much of the CD information occurs at higher spatial frequencies, there will be signalenergy beyond half that sample rate that will be aliaseci.A second limitation of the scan average estimate is that an MD estimate is available only aftereach traverse of the sensor. This limitation in the sample rate may be a factor limiting the frequencyrange of disturbances that can be rejected by the MD regulator.An enhancement of the above estimation scheme, used to increase the MD sample rate, involvessumming measurements over equal sub-portions of a scan.1 We will call this scheme Partial ScanEstimation. The MD estimate, calculated as the scanner reaches the end of one of the portions ofthe scan, is obtained by adding the average of measurements of the current scan portion to the mostrecent scan portion averages representing the remaining parts of the sheet. Figure 1.4 shows whichmeasurements are averaged for different positions of the scanner. In this example the scan has beendivided into three parts. With this method the CD variations will also tend to cancel out. However,1 Fmm conversations with Mike Heaven, Measurex-Devron Inc.,North Vancouver, Canada4Chapter 1: Introductionthe effect of the filtering is complex and is time-varying as it changes with each new scan portion. Itwill at best give the same pre-filtering as for the scan average method, and it uses older informationthan does scan average estimation.ScannerPosition .Estimation Averaging RegionEnd of scan____________________________________ScanEndoffirst IthirdofscanEnd of second I I I Ithird of scan ScanFigure 1.4: The variable filtering of the Partial Scan Estimation method1.1.5 EIBMC AlgorithmDumont et at. [5] describe an algorithm that for each moisture content sample taken along thetraverse of the sensor produces both a CD and MD estimate. In this way thousands of MD estimatescan be available during each traverse, increasing significantly the MD sample rate and potentiallyincreasing the range of frequencies for which disturbances can be rejected. The algorithm is knownas the Estimation and Identification of Moisture Content algorithm(EJMC). Wang[17] has improvedthe algorithm and also applied it to basis weight estimation. It has been renamed the Estimationand Idenqfication ofBasis Weight and Moisture Content algorithm (EIBMC). The main interest in thealgorithm so far has been in the CD estimates, for which it has proven superior to earlier methods.In this thesis the use of the EIBMC algorithm for MD control will be investigated.1.2 Description of the ProblemThe work described in this thesis helps answer the question: how can the EIBMC algorithm beused to improve MD control? Two ways in which the EIBMC will potentially aid our efforts are(a) by allowing for higher frequency MD estimates and (b) by providing flexibility in MD filtering.We show how increased sample rate and improved filtering can be realized through modification5Chapter 1: Introductionof the algorithm and provide analysis for determining when increasing the sample rate will be ofappreciable benefit.In the EIBMC algorithm MD variations are modelled as stochastic processes. Astrom[1] hasshown that MD variables such as basis weight and moisture content can be modelled as stochasticprocesses and that stochastic control methods can be applied. Bialkowski[3] has shown how thedisturbances are distributed in frequency. A stochastic model is a also a natural choice from thestandpoint of the control engineer as the variance of properties such as moisture content and basisweight is used by mill managers as a measure of paper quality.The sample rate and quality of filtering will be important factors affecting the quality of MDcontrol and so the frequency and quality of measurements made available by modern commercialscanners is a consideration. Present-day scanners report a set of the most recent measurements eachfive seconds2. Hence, MD control using the output of commercial scanners is limited to a fivesecond sample period.Because modern scanners provide samples every five seconds we will show how to modify theE1BMC algorithm to accommodate that data rate so that MD estimates are available at the samerate as measurements are reported by the measurement system. If a higher sample rate is requiredand one does not have the luxury of redesigning the measurement electronics, one may sampledirectly the raw un-filtered measurements. The EIBMC algorithm could then provide MD estimatesat the required rate. At least one manufacturer of CD controls equipment has tapped into the directmeasurement of commercial scanners in order to get faster and less filtered data3. When pre-filteringthis direct measurement, we show how to design the filter optimally and how to modify the algorithmto accommodate a spatial shift that occurs in the pre-filtering that would otherwise degrade both theCD and MD estimates. We also show that any error in the CD estimates causes errors in the MDestimates roughly proportional to the slope of the CD estimation error.Although one can expect the output variance to be reduced by increasing the sample rate, aquestion that arises is: what are the limits of this benefit? For a system with no time delay and finite2 Fmm conversatic4ls with Mike Yuswack of Measurex Inc., Ladner, BC, CanadaFor example, the Highspeed Universal Measurement Module (HUMM) of Measurex-Devron6Chapter 1: introductiondisturbance energy we can expect that the output variance can be reduced arbitrarily for sufficientlyhigh sample rate. However, a time-delay will ultimately limit the frequencies that can be controlled.Clearly there will be a point of diminishing returns, beyond which there will be little point inincreasing the sample rate. Of course, the practical limits on control effort is another limitationon variance reduction. Methods for the variance analysis follow to a great extent the work ofLennartson[8] regarding the effect of sample rate on closed ioop performance. In the analysis of thisthesis, however, the effect of the anti-aliasing filter on the dynamics of the closed loop system isalso considered. This consideration is necessary since for lower sample rates the anti-aliasing filterdynamics dominate the process dynamics.In order to determine the sample rate beyond which no further benefit is derived, the continuous-time average output variance is used as a criteria for comparing the performance of the closed ioopperformance for various sample rates. The variance analysis methodology of Lennartson has beenmodified to allow for the calculation of the actual process output variance as well as the varianceat the output of the anti-aliasing filter, since it is the output variance that is the true measure of thepaper quality. Calculations of optimal sample rates for typical paper machines have not, however,been included as numerical problems in the calculations have yet to be solved.A discussion on how the analysis methods can be applied to predict the potential product qualityimprovement due to process design changes is included. In particular it is proposed that the reductionof dead-time in the processes by moving measurements upstream may provide significant variancereduction.1.3 Organization of ThesisChapter 2 presents the paper variation models on which the EIBMC is based. Models formoisture and basis weight are given. Chapter 3 describes the EIBMC algorithm. Chapter 4 discussesstandard measurement methods and their effect on pre-filtering, describes how to filter with a Besselfilter to maximize the linearity of the phase delay, discusses how the EIBMC algorithm is modifiedto compensate for the delaying effect of filtering and provide high frequency MD estimates, and7Chapter 1: Introductiondiscusses the effect that CD estimation errors have on MD estimates. Chapter 5 contains the varianceanalysis methods and Chapter 6 concludes the thesis.8Chapter 2: Process ModthChapter 2Process ModelsThis chapter describes the process models used in the EJBMC algorithm. The algorithm wasdeveloped to measure moisture content components based on Lindeborg’s moisture variation model(MVM) (see [10]). Wang[17] has shown that the model, with some modification, can also beapplied to basis weight and, hence, the EIMBC algorithm is also suitable for measuring basis weightcomponents. We describe the process model for moisture content first and then discuss modificationof the model for use with basis weight.2.1 MoistureMoisture content in paper making has the property that as the Ml) value increases so does theamplitude of the CD profile. This effect is due to the fact that at low moisture levels the water bindsto the fibres and is less easily evaporated. Lindeborg has shown that the profile ampitude variesnon-linearly with the MD values and has modelled the process measurements as follows.y ==p+(l+Bp)uk+Vk (2.1)wheren is an index that indicates databox number and distinguishes positions inthe cross direction,k corresponds to uniform sample instants and distinguishes positions in themachine directiony is the part of the measured value at CD position n and time k that exceedsa reference value known as the target. The process variable setpointis typically used as the reference level.B is a measure of the non-linearity in the interaction between the CD andMD parts of the process, i.e., it is a measure of how the CD profileamplitude varies with changes in the MD part. B represents the9Chapter 2: Process Modelslinear term in a Taylor series expansion of the non-linear relationshipbetween the MD value and the amplitude of the profile. Lindeborg hasfound that using just the linear term provides a sufficient description.B has been found to be positive and so the profile amplitude increaseswith the MD value. B will of course vary according to the MD targetvalue around which the model is linearized.p is related to the CD value for position ii at time k. See the discussionbelow.Uk is the MD value at time k.17k contains measurement noise and model mismatch at time k and is assumedto be zero-mean normally distributed white noise with variance knownvariance R. Model mismatch refers to the difference between the truemoisture content relationship and the linearization as represented bythe parameter B.Figure 2.5 illustrates the indexing scheme for an N-position profile. Since the paper moves fromthe wet end to the reel, the MD indexes increase in the direction of the wet end relative to the paper.wet end MD reelN—•I,,• I,,•I,n+1—‘--‘-4-n — rn-i — ----• I,••k-ikk+1Figure 2.5: Process indexing according to Lindeborg’s model2.1.1 RemarksThe vector Ph = {p,j4,. . . ,p} is related to but not equal to the CD profile as defined in theopening chapter. The profile changes with MD according to the non-linearity in the moisture process10Chapter 2: Process Modthand equals Pk(1 + Buk), whereas Pk is invariant to MD changes when the changes are within theregion around the target for which linearity holds. To distinguish it from the true or instantaneousprofile Pk( 1 + Buk), k will be referred to as the pseudo-profile.The model assumes that MD changes affect only the amplitude of the profile and not its shape.This assumption is reasonable in that, due to the long time constants in the heating and cooling ofthe dry section rollers, we expect most of the moisture disturbances to be mostly due to variationsin flow and consistency that occur before the sheet is formed and, hence, would act on all parts ofthe profile equally.2.1.2 Process DynamicsAn assumption of the model is that P and B change slowly (See Lindeborg, 1986[l 1]). Pchanges slowly since it is determined mostly by the flow distribution of stock across the headboxand the wire. Moisture profiles can also be affected by the slowly varying condition of the felt thatcarries the sheet from the press through the drier section. B will vary only slowly as it is a functionof the fixed drying rate relationship of the fibres and water.The MD component is modelled as a first order process (See Natarajan[13]),(2.2)Here, ü is the mean deviation of (Uk] from the target and k is the first order process= k—1 + Wk (2.3)in which a, al < 1 is a known constant and wk is a zero-mean white noise sequence with knownvariance q. The parameters a, R, and q will of course depend on the paper machine and the samplerate.We now express the model in state-space form, a convenient representation for the EIBMCalgorithm.xl = Axk + Wk(2.4)yj: =p+Cxk+vk11Chapter 2: Process Modelswhereu 10 0Xk= , W,=0 a Wk (25)C= [(1+Bp) (1+Bp)]Note that this model does not take into account the effects of actuation. The addition of anactuation term would potentially decrease the estimation error, but would reduce the generalityof the algorithm in that more information would be required to make it work. In the pulp andpaper industries actuation and measurement suppliers are often different, making the sharing of suchinformation difficult. So, the present variation model is appropriate.2.2 Basis WeightWang[17] modified the EIBMC algorithm to include estimation of the disturbance parameters andapplied the new algorithm, called the Indentfication of Basis Weight Profile and ARMA MD modelalgorithm (IBPAM), to basis weight. It was found best to use a second order ARMA disturbancemodel for basis weight as it gave less prediction error than a first order model and no further benefitresults from increasing the order beyond second. It was also found that B tended towards zero,indicating as expected that MD basis weight variations do not affect the amplitude of the profile. Theprediction error variance did not differ significantly when B was fixed at zero. Hence, modificationof the EIBMC to accommodate basis weight estimation would involve increasing the order of thedisturbance model to second order ARMA and fixing B at zero. Instead of assuming that there isa constant MD offset from the target as in the moisture model, this state is absorbed in the ARMAmodel. The ARMA model is(i— aq1 a2q)uk = (i +b1q’ + b2q2)wk (2.6)In state-space for the basis weight model becomesXk+1 = Axk + WkUk = Cxk + Wk (2.7)yjj = ;: + Uk + Vk12Chapter 2: Process Modthwhere nowa 1 b1A= , Wk= Wk,a 0(2.8)C=[1 0]As for the moisture model, R and q, the variances of (vi and (w) are assumed to be known, as are theparameters of the disturbance al, a2, b1, and b2. In Wang’s algorithm these parameters are estimated.Otherwise these parameters maybe estimated independently and used in the EIBMC algorithm.13Chapter 3: The EIBMC AlgorithmChapter 3The EIBMC AlgorithmIn this chapter the EIBMC algorithm, developed to de-couple CD and MD components frommeasurements taken from scanning sensors, is presented. The version of the algorithm describedhere is a modification of the EIMC algorithm in that the optimal estimate of the MD component isextracted. The EIMC algorithm does not use the MD estimate as the algorithm has until now beenused only for CD estimation and as such it is only required to store a prediction of the next MD valuerather than the optimal estimate based on the most recent measurement. The next section describesthe how the indexing of the process model relates to timing of the scanning sensor. A descriptionof the EIBMC algorithm including a basis weight interpretation follows.3.1 Measurement Indexing and the ScannerA data indexing scheme was introduced in Chapter 2 for describing the process model uponwhich the EIBMC algorithm is based. The scanning sensor traverses the sheet with uniform velocityand takes measurements uniformly in time. Hence, the measurements will be spaced evenly across thesheet and evenly in time. Each measurement taken across the sheet in a single traverse correspondsto a unique databox n E {1, 2,3, ..., N}. Each new measurement also corresponds with a MD ortime index k e {..., —1, 0, 1, .. . }. Thus, with each new measurement the MD index increases by oneand the CD index either increases, for forward scans, or decreases, for reverse scans. Figure 3.6illustrates how the indexing scheme is interpreted for scanned data.3.2 The EIBMC Algorithm for MoistureThe EIBMC algorithm identifies the model parameter B as well as estimating the CD and MDcomponents of paper properties. Calculations are carried out in two time frames. The non-linearterm, B, is updated at the end of each scan whereas the CD and MD components are updatedon a measurement by measurement basis. A bootstrap estimation technique is employed involving aKalman filter for MD and a modified version of the recursive least squares with forgetting and resetting14Chapter 3: The E1BMC Algorithmalgorithm for CD. For a given measurement the MD and CD estimates are updated as follows. Thecurrent MD value is predicted based on the process model and the most recent MD estimate. This isthe prediction part of the Kalinan filter. Using this prediction and the current measurement the CDvalue corresponding to current databox is update using the Exponentional Forgetting and ResettingAlgorithm (EFRA) of Salgado, Goodwin and Middleton[15]. Then, given the updated CD estimate,the MD state is updated using the measurement update step of the Kalman filter.The moisture version of the EIBMC algorithm isXkIkl (3.9)kIk—1 = [1 1]XkIk_1 (3.10)= 1 + BU,IJC..1 (3.11)Ykk.-1 Pr—i +C1kIk_1 (3.12)1 r inirn’2v = —v—+ /3 — .5(Vz)2 (3.13)1+v()una,‘k Yk— YkIk—1j5=j3+1+(b)2v(3.14)=(1+E) [1 1] (3.15)K* —___________kCflECflT+R (3.16)CDFigure 3.6: Indexing for scanned data.15Chapter 3: The FJBMC AlgorithmE C7’CEEk = A (Ek_1— CEkCT:)AT + Q (3.17)XkIk = xkIk._1 + K(y— (p; + CXkIk_1)) (3.18)UkIk = [1 1]XkIk (3.19)The above calculations are performed for each measurement taken in a scan. The calculationsare done in the same order as the measurements are taken according to the indexing scheme. Inaddition to the above calculations, at the end of each scan, B is also updated and the mean deviationfrom target, ilk is modified to compensate for drifts in the pseudo profile average.A dynamic simulator, which will be discussed in Chapter 4, has been used to illustrate howMD estimates are improved by using the measurement update of the MD state instead of theprediction update. Figure 3.7 shows the MD estimator’s response to a MD step disturbance. Boththe measurement update and the prediction update from the old algorithm are given. Clearly themeasurement update gives the superior estimate.I I II pL --- jçJ8.4 ‘‘:..‘8.2—trueMD—— measurement update7.:prediction update7.60 2 4 6 8 10 12 14 16ScansFigure 3.7: MD estimate response to a step input, showing measurement and prediction updates16Chapter 3: The EIBMC Algorithm3.2.1 Estimation of BCalculation of B is based on how the variance of the process varies with changes in MD value.The process model isy =p+(l+Bp”)uk+vk (3.20)Letz= p + (1 + Bp)Uk (3.21)so that= 4 + V (3.22)From equation (3.21), it can be seen that(3.23)where a(n) and u(n) are the standard deviations of z and u, if the pseudo-profile is constant andBp’ is greater than —1.From equation (3.22) we get= u(n) + B (3.24)where a(n) is the standard deviations of y for databox n. Thus oz(n) can be calculated fromequation (3.24) and then (3.23) forms a set of linear equations in Ba and from which B can becalculated by the recursive least squares algorithm with forgetting.The standard deviation o(n) is updated using an exponential filter and (3.24) as followso(n,m) = pô(n,m —1) + (1— — fl)2 — B] (3.25)where the index m indicates the scan number, even meaning forward scan and odd meaning reverse,and p is the filter’s forgetting factor. The measurement mean, is also updated using an exponentialfilter.17Chapter 3: The EIBMC Algorithm3.2.2 Modification of the Mean Deviation from IkirgetAt the end of a scan, the pseudo profile has been updated but may not still sum to zero. Bydefinition, the excess represents MD variation. To correct, the average pseudo profile deviation fromzero is subtracted from the pseudo profile and added to ük.3.2.3 EFRAEquations 3.13 and 3.14 comprise the EFRA. The EFRA is a modified version of the well-known recursive least squares with forgetting factor (see Aström, 1990[2]), designed to prevent thecovariance matrix from growing either too large or too small. If limits are not put on the values inthe covariance matrix there is a danger that its values will exceed the ability of the computer to storethem. Inaccuracies could then cause failure of the estimation algorithm. The parameters of the EFRAare a, 3, A, and 6. A is the forgetting factor and a is the gain of the algorithm. /3 and 6 determine theupper and lower limits of the eigenvalues of the covariance matrix. V” is the covariance for databox n.3.2.4 Kalman FilterEquations 3.16—18, and 3.9 make up the Kalman filter. Equation 3.9 is the state prediction update.Equations 3.16 and 3.17 update the Kalman gain and covariance matrix respectively. Equation 3.18is the measurement update. The matrix Q contains the variances of the processes {u} and {}.Although an assumption of the moisture model that the mean deviation from target, {u}, is constant,the implementation must allow for changes in the mean deviation. Hence, this parameter is assigneda variance in the variance matrixqi 0Q=E{WkWI}= (3.26)OqThe variance q is assigned a small positive value, If qi is set to zero the Kalman gain for thisstate would go to zero and the algorithm would no longer respond to changes in the mean deviationfrom target.18Chapter 3: The FJBMC Algorithm3.3 The EIBMC Algorithm for Basis WeightThe basis weight version of the algorithm, modified to account for differences with the moisturemodel, isXkIk_1 = Ak_1Ik_1 (3.27)UkIk_1 = CkIk...1 (3.28)Ykk—1 = Pr—i + UkIk_1 (3.29)vn = :Vfl—+ /3 6(V)2 (3.30)a— (3.31)K*k—CECT +1? (3.32)= A (Ek_1-AT + Q (3.33)XkIk xkIk_1 + K*(y — (j3 + CXkIk_1)) (3.34)UkIk = (3.35)Q for this model isQ=E{WkWfl=b12 bcr, (3.36)where o, is the variance of Wk. Note that B has been set to zero.3.4 Parameter SelectionThe process models can be identified using any of a number of identification tools; see forexample Wang[17], or Panuska[14]. The parameters of the EFRA may be found by experimentation.Jonsson[7] has had success with the valuesa = 0.5/3 = 0.05(3.37)= 0.00005A = 0.9519Chapter 3: The EIBMC AlgorithmFor the variance of the mean deviation from target, q = 0.00001 has been used. Increasing thisvalue increases the rate at which the estimate responds to changes in the mean deviation.20Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmChapter 4Machine Direction Estimation Using the EIBMC AlgorithmAs previously noted the EIBMC algorithm was developed primarily for estimation of CDproperties but MD estimates are also available for each CD estimate. Modem scanners do notreport samples at the instant they are taken but instead deliver them in sets at the end of a sub-periodof the full-sheet scan time. The rate at which scanners report measurements limits the rate at whichMD estimates, produced by the algorithm, can be used for MD control. Additionally, the user islimited to the boxcar window filter implicit in the sensing method.Higher frequency MD estimates can be achieved by bypassing the scanner electronics and tappingdirectly into the sensor signal. Doing so provides not only high frequency estimates but also allowsfor greater flexibility in filtering. Filters superior in performance to windows can be used but at thecost of using a filter with poles. The poles will cause the signal to be delayed. The magnitude ofthe delay will vary with the frequency components in the signal and the variable delays will causea smearing of CD features present in the signal. The spatial effect of smearing will be doubledsince the scanner takes measurements in both directions across the sheet. However, by using a filterthat minimizes the variability in the delay, namely the Bessel filter, and by shifting the data vectorspatially a suitable amount in the EIBMC algorithm, the distortion of CD features is minimized.The Bessel filter is designed to have a time delay equal to the time it takes for the scanner to passover an integer number of databoxes. The raw samples are then shifted back this integer number ofdataboxes before being presented to the EIBMC algorithm. This process of shifting scanned data thenumber of databoxes corresponding to the delay of the filter, prior to being presented to the EIBMCalgorithm, will be referred to as delay compensation.In this chapter we show, with the help of a scanner/sensor simulator, the effect of poor filteringand no delay compensation on CD estimates and how poor CD estimates cause errors in the MDestimates. We show how to filter directly sensed measurements optimally and how to modify theEIBMC algorithm to compensate for delays. We show that CD estimation errors such as the onescaused by smearing result in MD errors proportional to the CD estimation error and we show how21Chapter 4: Machine Direction Estimation Using the EIBMC Algorithmto further modify the ETBMC algorithm to provide MD estimates at any period that is a multipleof the sensor sample period.The chapter is organized as follows. In section 4.1 we present the report timing of typicalcommercial scanners along with the filtering effect of the sensing method. In section 4.2 we discussdirectly sampling the sensor output. Section 4.3 describes the design of the scanner simulator. Section4.4shows how raw samples are shifted to compensate for filter time delays. Section 4.5 discussesthe effect of pre-filtering on the signal. Section 4.6 shows how to compensate for spatial smearingthrough pre-filter design and modification of the EIBMC algorithm and section 4.6 shows how toalter the rate of MD estimates.4.1 Commercial ScannersThis section describes how a typical commercial scanner would sense and report measurementto the control system and what effect the sensing method has on pre-filtering the signal.4.1.1 Partial Scan ReportingFigure 4.8 illustrates the movement of the sensor relative to the paper. The sensor moves backand forth on the 0-frame tracing out a zigzag trajectory. Samples are taken uniformly in time asthe sensor passes over the sheet and report a vector of the most recent samples at regular intervalsduring the scan. The period of time between reports will be referred to as the report period andthe set of samples reported at the end of the report period will be referred to as the mini-scan. Atypical report period for present day scanners is five seconds with a scan time as low as ten secondsand as high as sixty seconds (see Havelock[6] and Lindeborg, 1986 [11]) depending on the width ofthe sheet and the speed of the scanner. Today’s fast scanners will move between twelve and sixteeninches per second. The number of databoxes available on current scanners can vary from tens tothousands depending on the application.4.1.2 SensingLater in this chapter we show that the accuracy of the CD estimation affects the accuracy ofthe MD estimation. Before the signal is sampled, high frequencies must be eliminated to prevent22Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmFigure 4.8: Sensor-sheet orientation and measurement timingaliasing. An important aspect in the use of the EIBMC algorithm, affecting the quality of the CDestimates, is the manner in which data is filtered prior to being presented to the algorithm. Basisweight and moisture measurements are sampled by integrating their signals over the sample period.Basis weight for example is measured with a beta gauge in which the mass of fibres and water isdetermined by the degree of absorption of a stream of electrons passing through the sheet. Basisweight is a function of the current that reaches a detector on the other side of the sheet. A currentproportional to the intensity of the stream of beta particles is amplified and converted to a frequencyusing a voltage controlled oscillator (VCO). The pulses are then counted over the sample period todetermine basis weight. The effect of counting the pulses is to average the basis weight over thesample period. For example, let c(kT8)be the count over the sample period ending at time kT8. Thedifferential change in the count due to a differential change in time isdc(t) = Kb(t)dt (4.38)where b(t) is the basis weight and K is the ratio of changes in the VCO frequency to changes inbasis weight. Thus, integrating (4.38) over a sample period gives, excepting discretization errorsSensorJ.-0-frameProcessing finstant for 2---------------miniscan 123Chapter 4: Machine Direction Estimation Using the EIBMC Algorithminherent in the counting of pulses,C(kT3)= K J b(t)dt (4.39)(k1)THence the count is proportional to the average basis weight over the sample period.4.1.3 The Filtering Effect of the Sensing MethodThe filtering effect of averaging the raw signal over the sample period is equivalent to acontinuous-time boxcar window. Note that the filtering inherent in the counting of pulses to obtainsamples of the raw signal is similar to the filtering inherent in taking an average of the measurementsin a scan to obtain MD estimates. Both are boxcar filters. However, in the present case we obtainraw samples from the sensor containing both CD and MD components at the sample rate of thesensor. In the case of the scan average we are averaging the raw measurements taken during a scanto obtain an estimate of the MD component one per scan.Often only first-order filtering is used in the sensor amplifier4. The combined effect of theamplifier filtering and the averaging inherent in the counting form the anti-aliasing filtering with thisscheme. The filtering due to averaging over the sample period is poor at rejecting frequencies abovehalf the sampling rate as was shown for the boxcar filter in chapter one (see Figure 1.3). A firstorder filter in the amplifier has poor magnitude response and phase response; linear phase responseis particularly important for reasons that will be clear in the following sections.4.2 Direct MeasurementsCurrent scanners are limited in the rate at which they deliver measurement reports and do notallow for flexibility in filtering. An alternative is to tap directly into the sensor output and do filteringand sampling directly5. By using a filter whose transfer function has poles, better filtering can beachieved in the sense of rejecting frequencies above the Nyquist rate. As we will see, it is alsodesirable to chose a filter that gives equal delays to all frequencies in its passband to minimize spatialdistortions.Fmm conversations with Karl Hahn of Measurex Inc., Cupertino, California, U.S.A.Devron-Hercules produced a product for this purpose in 1990— the High-speed Universal Measurement Module (HLJMM)24Chapter 4: Machine Direction Estimation Using the EIBMC Algorithm4.3 Scanner SimulatorThe dynamic system simulation software MatrixXlSystemBuild was used to construct a computersimulation of a typical scanner for the purposes of investigating the effects of filtering on CD andMD estimates. The scanner simulator is part of a paper machine simulator that also includes, forexample, consistency control, coarse and fine stock flow control, steam pressure control, and machinechest level control. The software allows functional blocks to be interconnected graphically. C codemay be used to define custom designed blocks which are then linked to the SystemBuild image ifieand interconnected as regular blocks.The simulator has been set up as the scanner for a typical newsprint machine. The scan time istwenty seconds with 18.8 seconds on-sheet. The sensor report period is 4.7 seconds, and the off-sheettime is 2.4 seconds. There are 60 databoxes and four reports per scan, each containing 15 databoxes.Figure 4.9 shows the structure of the scanner simulator components. The simulator is designedas follows. To simulate the paper sheet an MD signal is added simultaneously to each component in avector of CD values. The sum is known as the sheet vector. The number of components in the sheetvector is variable and each component may be interpreted as a databox or else several componentsmay correspond to a databox depending on the desired resolution. High resolutions may be useful inthe future for studying CD control or wave effects on the wire. There is a block that simulates therail potentiometer which gives a signal indicating the position of the scanner on the sheet. Anotherblock, known as the sheet width, indicates the locations of the edge of the sheet. The combinationFigure 4.9:25Chapter 4: Machine Direction Estimation Using the EIBMC Algorithmof the rail pot. and sheet width give flexibility for easy changes in sheet width, simulation of sheetwander, and various changes in sensor timing. To simulate the scanning action, a multiplexer choosesone of the sheet vector components. The choice of component is made based on the position of thescanner as indicated by the rail pot. The output of the multiplexer represents raw sensor data.The multiplexer output is then fed to a block simulating the anti-aliasing filter. This is a genericSystemBuild block that simulates a rational transfer function. One need only specify the coefficientsappropriately. The output of the anti-aliasing block is fed to a de-multiplexer. The de-multiplexerreads the value of the rail pot. and the ifitered signal and places a sample of the signal value in avector in a position corresponding with the position of the scanner on the sheet as determined by therail potentiometer. The de-multiplexer block represents the part of the sensor data processing thatcreates vectors of the incoming data for transmission, for example, in the form of mini-scans, to thecontrol system. A second function of the de-multiplexing block is to signal, to the EIBMC block,the ends of mini-scans and to indicate the mini-scan number.No multiplexing or de-multiplexing functions are available in the SystemBuild block library.Initially the multiplexer was implemented by simulating digital counters and analog transmissiongates. However this method proved too time consuming at nm-time. The multiplexing and demultiplexing functions were implemented in C-code blocks and linked to SystemBuild. The blocks aredesigned to determine the size of the input and output vectors automatically and process accordingly.The EIBMC algorithm is also implemented in C-code. The code was ported from an existingstand-alone program for testing and developing the algorithm. It was modified to process data inmini-scans. The block wakes up when it receives an interrupt from the de-multiplexer block. Itdetermines the current mini-scan position and processes the portion of the de-multiplexer outputvector corresponding to the current mini-scan, producing the CD and MD estimates for that mini-scan. The raw values as well as the estimated CD and MD values are logged to the hard drive bythe block for latter analysis. The estimates are then fed to a graphic driver that scales and translatesthe values for display with a graph block that has been created for this simulator.Figure 4.10 shows a captured display from the scanner simulator. In this simulator there aretwo sensors on one scanner: one for basis weight and one for moisture. Four graphs are shown.26Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmClockwise from the top right are the CD estimate, the MD estimate, the raw data, and the true CD.The graphs are updated in real-time. The simulation runs at about the same rate as a real scanneron a Sun SPARC Station 2.4.4 Data Shifting in the EIBMC AlgorithmFigure 4.11 shows, for the example of the simulator, the scheme for shifting data for a forwardtraverse as the mini-scans are made available from the scanner. After the first mini-scan, the currentsamples are shifted back two databoxes. The first two samples are not used. Only 13 databoxesare used for estimation since samples corresponding in space to the last two databoxes have beendelayed in time and will not be available until the next mini-scan is reported. All 15 samples of eachof the remaining mini-scans, 2, 3, and 4, are used for estimation but are shifted back two databoxes.There are no samples corresponding to the last two databoxes, but there will be samples availablefor estimation for these databoxes on the reverse traverse. The reverse scan is treated similarly.For traditional CD control, data shifting scheme is most important for systems for which theratio of the number of databoxes to the number of actuator zones is low. If the ratio is high then thepre-filter will of necessity have a delay that is small relative to spacing of the actuators and so therelative spatial phase distortion will be small. Of course, for high frequency MD estimates, accurateCD estimates will be required at the spatial frequency corresponding to the temporal MD estimatefrequency. The CD estimates can then be mapped on to actuator zones as for traditional CD control.4.5 The Effect of Pre-Filtering on Spatial SmearingThe purpose of the pre-filter is to prevent aliasing. However, all filters will also introduce a timedelay that may vary for the different frequencies that make up the signal. Figure 4.12 shows timedelays, as expressed spatially in databoxes, for sixth order Bessel and Butterworth filters, along withmagnitude plots. The transfer function of the Bessel filter shown is, from Williams[18]26.76 5 4 3 2 1 (4.40)() + 7.79(e) + 28.8(.-) + 64.O(_) + 88.8(-) + 72.3(_) + 26.727Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmIJiIi. Lr!i.dr!r!r!rwIIII.CJ.OrIIOIJOII333 333‘r•‘11231 31I” 12 1121!II21112 33 2221111111222.222131311SI_Lr[iII !ri1___ ___________BryiI LiI_ _ __ _ _ ___‘____1____•—.r:riIllhiIr_____ _____ru__Figure 4.10: SystemBuild scanning sensor simulator.28Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmForward Scantimemini scan mini scan mini scan mini scan1 2 3 4I...I/_____1// ///// ///[.measurement_____________________________________ _____________after shift___ _ ______ _ ___I IProcessing after 1st after 2nd after 3rd after 4th nottime trigger trigger trigger trigger usedFigure 4.11: Data shifting scheme for forward sensor traverseswhere, for this example,-‘B = 4.32 and for a Butterworth filter16 5 4 3 2 1 (4.41)(.) +3.88(t) +7.5o(.) +9.19(t) +7.5o() +3.88(-) +1where ‘B = 6.21. These filters were designed to give group delays corresponding to two databoxes.The details of the filter design are given in section 4.6.2. The time delay is expressed spatially interms of databoxes rather than in time units for reasons that will be clear later. Note that in the case ofthe Butterworth filter, the time delay due to the phase shift can vary significantly over frequencies inthe passband. The effect on profile estimation is that features in the profile are, to a degree dependenton the frequencies in the features, smeared spatially. As the scanner passes back and forth over thesheet, smearing occurs in both directions preventing the estimation routine from capturing an accuraterepresentation of CD properties. One can avoid spatial smearing by choosing a filter that delays allfrequencies the same amount and then shifting the data the number of databoxes corresponding to thedelay before the data is presented the to the EIBMC algorithm. The Bessel filter gives the best phaseresponse since this filter gives the maximally linear phase response. Note that although the Besselfilter gives superior phase response, the magnitude response is not as good as for the Butterworthfilter. The Bessel filter magnitude rolls off more slowly at the cut-off frequency.4.5.1 CD DistortionTo illustrate the effect of not compensating for the effect of the delays and variability of the delays29Chapter 4: Machine Direction Estimation Using the FJBMC Algorithm10’ 100 101 102Frequency (rad!sec).E -1‘U‘U0-1.510_i 10° 101 102Frequency (rad/sec)Figure 4.12: Sixth order Butterworth and Bessel filter Bode Plots with time delay expressed in databoxescaused by the anti-aliasing filter, the following experiment has been conducted. In the simulator a CDprofile has been set up that includes a spatial step and pulse, intended to approximate the response toa basis weight or moisture actuator change, as shown in the bottom of figure 4.13. The MD signalis a sine wave with a period of two hundred seconds or ten scans. In this experiment a Bessel filteris used for anti-aliasing. Two cases are considered: with delay compensation and without.The top graph in figure 4.13 shows the output of the anti-aliasing filter for a forward and reversetraverse referenced to the position of the scanner. In this instance the reverse traverse preceded theforward traverse so that the time begins on the solid line at the right hand side of the page andproceeds to the left, turning around at the end and continuing on the dashed line. The lower graphshows the true CD profile which clearly influences the raw measurements. The influence of the MDcomponent, which in two scans passes through a fifth of a cycle, can also be seen in the top graph.In the two scans represented here the slope of the MD component is negative.Figure 4.14 shows a detail of the spatial impulse in Figure 4.13, illustrating the extent of the timedelay of the filter. For easy comparison the mean has been removed from the measurements. It is seenthat the reverse traverse shifts features to the left two databoxes and the forward traverse shifts datatwo to the right. The shifts correspond to the two databox time delay designed into the Bessel filter.30Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmFIeroutputdataboxTnieCDdataboxFigure 4.13: Spatial shifting of raw measurements due to filter delayDetaf of feature shiftingC —— ——,———tn,eCD--revemeH f ward.1r______L_————48 49 50 51 52 53 54 55 56 57 55databoxFigure 4.14: Detail of the Irue CD and raw measurement with mean removedThe CD estimates are shown in Figure 4.15. For the case without delay compensation, the effectcan be seen in two lobes on either side of the step, and the pulse has been smeared in both directions.With delay compensation the lobes and smearing disappear.As has been shown, the Bessel filter has a more linear phase response than the Butterworth filter.Because of the variability in the phase response of the Butterworth filter, CD features are smearedaccording to their frequency content. Such a distortion will make control of the higher frequencies31Chapter 4: Machine Direction Estimation Using the EIBMC Algorithm000020 30 60CD PositionFigure 4.16: CD response using sixth order Bessel and Butterworth filters with data shiffingimpossible because that information is lost in the smearing. In Figure 4.16 CD estimates, for whichthe Bessel and Butterworth filters were used, are compared. Delay compensation is used in bothcases. It can be seen that the estimates for the Bessel filter are superior to those of the Butterworthfilter, which produce overshoot at the step. There is also a slight smearing of the pulse.30CD PositionFigure 4.15: CD response using sixth order Bessel filter with and without data shifting32Chapter 4: Machine Direction Estimation Using the EIBMC Algorithm4.5.2 The Effects of CD Estimation Error on MD EstimationThe effects of poor CD estimates on MD estimates can be seen in figure 4.17. Large spikesare present at times that correspond with the step in CD and smaller spikes correspond with the CDpulse. Spikes occur when shifts are not compensated for since CD estimates tend to average the CDportion of the measurements but, because the measurements are shifting back and forth with eachscan, there is a discrepancy wherever the CD profile is changing. The EIBMC algorithm assignsmost of this discrepancy to the MD estimate and the effect is that the MD estimates contain a portionof the CD gradient. The figure also shows that the MD estimation is greatly improved by correctdelay compensation. Note that there is a phase shift and an amplitude reduction in the estimate dueto the inertia of the MD portion of the estimation algorithm.Figure 4.18 shows MD estimates for the first fifty scans when the Bessel filter and data shiftingare employed. This figure illustrates that the MD estimates tend to settle within 10 scans after alarge profile change or from initial conditions.Tmie, secFigure 4.17: MD response using sixth order Bessel filter with and without data shifting33Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmFigure 4.18: First 50 scans of MD response using sixth order Bessel filter with data shifting4.6 Compensation for Pre-Filter DistortionsThe pre-filters compared are both sixth order and the delay of each corresponds to exactly twodataboxes. A backward of two databoxes for both the forward and reverse traverses will allow forthe alignment of CD features before the samples are presented to the ETBMC algorithm.4.6.1 Maximally-Linear-Delay FilterThe smearing problem can be alleviated by using a pre-filter with a linear phase response so thatthe filter will have a constant delay for all frequencies. In this way the true shape of features withfrequency content below the cut-off will pass through the filter unchanged except for a time delay. Itis also required that samples be properly oriented spatially, taking into account the time delay, beforebeing presented to the estimation routine, so that samples from a given position on the sheet, taken onforward and reverse traverses of the sensor, properly align. The filter that has the most linear phaseresponse, and hence the one to best preserve spatial features, is the Bessel filter (see Williams[18]).Figure 4.12 shows how the time delay varies with frequency for Bessel and Butterworth ifiters. Inthis example a second order Butterworth has been chosen for comparison as it has the same order asScans34Chapter 4: Machine Direction Estimation Using the EIBMC Algorithmthe simplest Bessel filter. A practice in some commercial scanners is to use a first order filter at thesensor output which has poorer phase and magnitude response6.In order to compensate easily for the filter’s delay, it is chosen to be an integer multiple of thesample period. Note that as the scanner traverses the sheet the sample period corresponds spatiallywith one databox width. By choosing the delay to be an integer multiple of the sample period, we cancompensate for the delay by shifting the measurements before processing with the EIBMC algorithm.The delay of the filter is a function of the filter order and the cut-off frequency. A common criterionfor choosing the cut-off frequency is to specify the attenuation at half the sample rate. The attenuationfactor, however, will depend on the magnitude of process disturbances above the half the sample raterelative to those below. When this factor is selected, the ifiter order and cut-off frequency can beadjusted to arrive at a delay that represents an integer multiple of the sample period.The time delay or group delay can be calculated as (w = 0) where ç& is the phase of the filter’stransfer function and w is the frequency in radians per second.4.6.2 Bessel Filter Design ExampleWe illustrate the relationship of the ifiter order and cut-off to time delay with an example showingthe design of a sixth-order Bessel filter to give the appropriate integer delay. We then discuss howto determine the attenuation properties of the filter.The generalized form of a sixth-order Bessel filter is26.76 5 4 3 2 1 (4.42)() + 7.79(e) + 28.8(-) + 64.0(-) + 88.8(_) + 72.3(-) + 26.7where cs.’B is the cut-off frequency. For low frequencies the filter can be approximated asH(jw) 26.7 (4.43)72.3(f) + 26.7The phase of H(j)at low frequencies is/72.3warg(H(jw)) = arctan\26.7wB/2.71w= arctan I\ WB6 Fmm conversations with Karl HaJin of Measurec Inc., Cupertino, California, U.S.A.35Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmThe time delay as a function of frequency is thusTd(w) = _{arg(H(jw))}d I f2.71= —arctan(dw I \ W (4.45)2.71i+()2Evaluated, the group delay at DCTd = —-i (4.46)WBFor the scaimer simulator, sixty databoxes are sampled during a 18.8 second scan. Thus thesample rate is T8 = 0.313 s. If we choose to have our anti-aliasing filter cause a two databoxdelay thenTd=2Ts=2 71WB (447)= = 4.32 rad/sChoice of the filter order is a trade off between attenuation, and group delay and filter complexity.We can form a table of attenuation for various filter orders and delays corresponding to integernumbers of databoxes. Let ,j = be the number of integer delays. Then2.71 2.71 (4.48)Td iiTHalf the sample rate, i.e., the Nyquist frequency, is = -. In equation (4.42) may be replacedby= zç- =i.E =jl.1677 (4.49)For all orders of Bessel filters the normalized frequency at the Nyquist frequency may be expressed as= jat where a will be referred to as the orderfactor. The cut-off frequency is then expressed asWB = —s- (4.50)36Chapter 4: Machine Direction Estimation Using the EIBMC AlgorithmlOrder 12 14 16 IOrder factor 2.42 1.5 1.16Table 4.1 Bessel filter order factorsInteger delay Filter orderjj 2 4 61 0.2363(1.3) 0.4243(2.1) 0.6218(2.7)2 0.0664(2.6) 0.0554(4.2) 0.1057(5.4)3 0.0302(3.9) 0.0120(6.3) 0.0129(8.1)4 0.0171(5.2) 0.0039(8.4) 0.0025(10.8)Table 4.2 Bessel filter attenuation at half sample rate with the ratioof Nyquist frequency to cut-off frequency, w/w, in brackets.Some order factors for low order Bessel filters are given in table 4.1. To find the attenuation of thefilter at half the sample rate we substitute (4.49) into (4.42).Table 4.2 summarizes the attenuations for various filter orders and integer delays along with theratio of the Nyquist frequency to the cut-off frequency. The attenuation and cut-off frequency givethe most important characteristics of the filter as they indicate how much information will be lostby filtering and how much noise will be added by aliasing. Design of the appropriate filter theninvolves choosing from the attenuation table the delay and ifiter order according to the attenuation andbandwidth appropriate to the application, determining the order factor, a, for the chosen filter orderfrom table 4.1, calculating WB using (4.50), and then parameterizing the filter using the normalizedform for the given filter order. For example, (4.42) is the normalized form for the sixth-order Besselfilter.4.7 MD EstimationThe EIBMC algorithm produces one MD estimate for each databox during a scan. For thepurposes of control, however, this sample rate may be faster than is desired when one considerssuch factors as simplicity in the control design for the process dynamics and dead-time, and wearon the actuators. This section shows how to produce MD estimates with a period that is any integermultiple of the sensor sampling period. In the case of using the reports that commercial scanners37Chapter 4: Machine Direction Estimation Using the EIBMC Algorithmprovide several times per scan, the minimum MD estimate rate usable for MD control would be onceper report. In decimation, i.e.,changing the effective sampling rate of a discrete-time signal, it is notenough to simply choose every n’th sample because the original signal may contain frequencies thatwould be aliased in the new sequence. Thus, a digital anti-aliasing filter is required to attenuate thosefrequencies above half the new effective sample rate. Once again the choice of filter will depend onthe relative amount of energy in original sequence above half the effective sample rate of the newsequence. An example based on the simulator parameters follow.Suppose that we are limited to using the data reported by the scanner each 5 seconds. For eachmini-scan we have 15 values corresponding to 15 databoxes. But, because data is only availableevery 5 seconds, it is only possible to take a new control action every 5 seconds. Hence it will benecessary to compress the 15 MD estimates corresponding to the 15 new values down to one MDestimate representing the current 5 second interval. The new MD estimate is arrived at by using adigital filter to reduce aliasing and then sampling the filtered stream. The digital filter, obtained byapplying a bilinear transform to a second order Butterworth filter, designed with cutoff at one eighthof the new sample rate, is:0.000661 + 0.001322q1+ 0.000661q2U,k= 1 — l.927q’ + 0.9293qUk (4.51)lii our case we take the last value in the filtered stream as the MD estimate for that 5 second intervalas it contains the most recent information. This filter is easily implemented in the EIBMC algorithm.The last two MD estimates and filtered MD estimates are stored. When a new MD estimate isavailable the new ifitered estimate is calculated according to the filter parameters. The current valuesare stored and the cycle repeats. A routine to implement the filter can be attached to the EIBMCalgorithm with little modification to the EIBMC code. For example, if the programming language isC, a call to a filter function could be placed directly after the point at which the MD estimate for thecurrent databox is available. The filter routine would store past values in static variables within thefunction so that the main EIBMC would not even have to have new variables declared.Figure 4.19 shows the dynamics of the filtered and re-sampled MD estimates in terms of theestimator’s response to a MD step change. The time constant is about two thirds of a scan. There is38Chapter 4: Machine Direction Estimation Using the EIBMC Algorithmovershoot as is to be expected for a Butterworth filter. The spikes in the high frequency MD estimatesare from CD estimation errors caused by large CD gradients as shown in previous examples.8.c8.48.28 —trueMD—— tittered and resampledhigh frequency estimation7.87.60 2 4 6 8 10 12 14 16Figure 4.19: MD filtered and re-sampled every five secondsAlthough a Butterworth filter is used in this example, there are no restrictions in the choice offilter type or filter order. One must observe, of course, that aliasing is eliminated by choosing afilter that attenuates sufficiently beyond the Nyquist rate. A simple and effective design method isto start with a continuous-time filter design and convert the filter into discrete-time using a bilineartransform. Astrom[2] may be consulted for bilinear transform methods.4.8 SummaryWe have shown the shortcomings of present day scanners in the way that measurements areprocessed and delivered to the control system. We have shown that poor filtering causes CD estimationerrors that in turn cause MD estimation errors proportional to the gradient of the CD errors. We haveshown that in order to avoid CD estimation errors a linear phase filter should be used. The filterwith optimally linear phase, the Bessel filter, was presented. A scanner simulator was described andsimulations illustrating the effect of poor and good filtering shown. It was shown how to compensatein the EIBMC algorithm for the time-delay due to the Bessel filter and how to design the Bessel39Chapter 4: Machine Direction Estimation Using the EJBMC Algorithmfilter to have a delay corresponding to an integer number of databoxes so that compensation withthe EIBMC filter was simplified. Finally, we discussed how to filter the MD estimates so that theymay re-sampled at lower rate.40Chapter 5: Continuous-Thne Variance Ana1ysL of Sampled Time-Delayed SystemsChapter 5Continuous-Time Variance Analysis of Sampled Time-Delayed SystemsIn devising a scheme for increasing the rate at which MD estimates are provided, we must ask ifany benefit can be expected. In this chapter we describe the development of algorithms for analyzingthe performance of a controlled, sampled, time-delayed, continuous-time stochastic process as thesample rate is varied. For non-delayed processes with finite energy disturbances, the output variancecan be reduced arbitrarily by choosing a sufficiently high sample rate. However, with a time delaythere will be a limit to output variance reduction; there will be a sample rate beyond which noadditional variance reduction can be achieved.In addition to determining the optimal sampling rate, we also address the effect of process designchanges on the output variance. For example, in order to reduce the time-delay in the basis weightprocess, it has been proposed to sense basis weight either after the dry line or at the slice lip, theopening of the headbox from which the stock is deposited on the wire. Measuring at the lip couldreduce the time delay as much by as 45 seconds and, since the time delay is one of the main factorslimiting variance reduction, could allow for superior performance of the basis weight ioop.The chapter proceeds as follows. We first determine a criterion for comparing the performanceof the system at various sample rates and then show how the analysis is formulated. The continuous-time average process variance is chosen as the criterion and the dynamics of the anti-aliasing filterare included in the process model, since, for part of the range for which sample rates are considered,the anti-aliasing filter dynamics are significant compared to the plant dynamics. The analysis is thenpresented. Because not all numerical problems have been solved for the variance calculations, nonumerical results for time delayed plants are presented. An example of variance analysis is given,however, for a system with no delay, and the numerical difficulties are discussed. It is seen that thevariance at the output of the anti-aliasing filter can differ significantly from the process variance andthat for systems with no delay the process variance may be made arbitrarily small for sufficientlyhigh sampling rate.41Chapter 5: Continuous-fine Variance Analysis of Sampled lime-Delayed Systems5.1 Criteria SelectionIn studies to determine the effect of sample rate on control system performance, McGregor[12]investigates the performance of discrete time minium variance controllers over varying sample rates.But his analysis does not consider the activity of the process between control actions since he usesas a criterion for comparison the variance of the measurement at the sample points. There are twodeficiencies with this approach. Firstly, the process variance can vary significantly between thesample instants and can exceed by large amounts the variance at the sampling times. The criteriondoes not account for this behavior and will not accurately reflect the quality of the process response.Secondly, this approach does not consider that an anti-aliasing filter must be used before samplingto attenuate frequencies above half the sample rate, since a ifiter will hide disturbances from thecontroller and, hence, from the variance calculations. The frequencies of the disturbances hiddenthis way will vary as the sample rate is varied and thus the criteria biases the variance according tothe sample rate and the frequency distribution of the disturbances. In addition, the dynamics of theanti-aliasing filter may be significant compared with the dynamics of the plant and may dominate.This is the case for low sample rates such as scan averages for both basis weight and moisture.To illustrate these effects, consider that 30 seconds is a typical time constant for the basis weightresponse to thick stock valve changes. Suppose the scan time is 20 seconds and we required 20dBattenuation at half the sample rate to avoid aliasing. Using a sixth-order filter would require a cut-offfrequency of 1/60 Hertz and the anti-aliasing filter dynamics would dominate those of the plant.Hence, we must consider the dynamics of the anti-aliasing filter when we compare the performanceof the control system over a range of sample rates that approach the scan average.A suitable criterion for comparing control performance over varying sample rates should be basedon continuous-time and on the process output, not the measurement. Lennartson[81 has recommendedthe use of the continuous-time average variance, which we will adopt for the analysis. The criterioncan be expressed as the expected value of(k+1)T9_ f {y(t)}2d (5.52)kT8where T3 is the sample period and y(t) is the process output.42Chapter 5: Continuous-Thne Variwzce AnaiysL of Sampled Thne-Delayed SystemsZk5.2 Overview of the Variance CalculationsThe goal of our analysis is to determine the minimum achievable continuous-time variancefor a range of sample rates. In this way we have a tool for identifying optimal sample rates forexisting process configurations and for predicting the performance of new process configurations.The approach we take is to design a controller to minimize the continuous-time average varianceand then calculate the continuous-time average variance for the closed-loop system. A predictivecontroller will be used as recommended by Lennartson[8]. The state of the process is estimated by aKalman filter and then the predictable state, i.e., the state of the process after the time delay (fractionaldelays are allowed), is predicted according to the process model and the most recent control actions.A control signal is then calculated as for a system with no delay. In this way controller design isindependent of the time delay. An advantage to this choice of controller is that the poles of the closedloop system are those of a system with the same dynamics but no delay plus those of the Kalmanfilter plus poles at the origin corresponding to the number of lags in the delay. Another advantage isthat the analysis involves matrices of the dimension of the system without the delay. Thus, numericalwork is simplified. Traditional approaches, as described in Aström [2], to simplify the varianceanalysis, reduce the problem to a system without time delay by expanding the state to include pastcontrol actions. That approach involves matrices on the order of the system size plus the number oflags in the delay. Figure 5.20 shows the structure of the controller. The predictable state is defined asx(t) = x(t + r) (5.53)where r is the time delay.Figure 5.20: Controller structure43Chapter 5: Continuous-Thne Variance AnalysLc of Sampled Time-Delayed SystemsThe analysis proceeds as follows. First, a continuous time state-space model that includes theplant, anti-aliasing filter, and disturbance is expressed such that the process output and anti-aliasingfilter output are both easily extracted. With this formulation the statistics of interest, namely thevariance of the process and measurement are easily calculated. Next, the continuous-time modelof the plant, the anti-aliasing filter, and the continuous time control criterion are sampled to createa discrete-time model of the plant, anti-aliasing filter and disturbance. The predictor and Kalmanfilter are presented and a controller is designed to minimize continuous-time average variance. Thechapter concludes with the variance calculations.5.3 Continuous-Time ModelThe plant and controller block diagram is shown in figure 5.21. From the diagram, u(t) is thecontrol signal, y(t) is the process output, ym(t) is the output of the anti-aliasing filter, Ym,k is theoutput of the anti-aliasing filter at the sampling instants, T3 is the sample rate, vk is measurementnoise, zk is the sampled measurement, n(t) is the process disturbance, and e(t) is continuous-timewhite noise and as such its time integral,w(t)=(5.54)has independent increments which for this system will have variance Rdt where 2irRis the spectraldensity of e(t). The measurement noise, vk has covariance R2The formulation of the state space equations is conducted as follows. Two state space representations are created. One includes both the plant and disturbance and the other includes the anti-aliasingfilter. To extract the output of the anti-aliasing filter and the process output, y(t), for the purposesof simplifying the variance calculations, the process state, W,1, is expressed in observable canonicalform and the filter state, Waa, is expressed in controller canonical form. The two components arethen expressed in a single state space equation.44Chapter 5: Continuous-lime Variance Analysis of Sampled lime-Delayed SystemsYm(t)Figure 5.21: Mathematical Model of Plant and Controller5.3.1 Plant and DisturbanceThe continuous-time model of the plant, from the signals e(t) and u(t) to y(t), expressed intransfer function form isAt,(s)Y(s) = Btf(s)U(s) + Ctf(s)E(s),Atf = S +a1s’ + ... + a,_s + an,(5.55)Btf =b18’— +b2s+ + b,_1s+ b,Ctf = C1S” +C2S’ + + C,_1S + C,zwhere the polynomials B, and C have no common factors. Y, U, and E are the Laplacetransforms of y(t), u(t), and e(t) respectively. The state space model for the two-input system isobtained by first forming observable canonical forms for the systems that map U onto Y and E ontoY. The choice of observable canonical form makes it easy to add the two systems using the principleof superposition to obtain the state space model for the plant and disturbance. With—a1 10••. 0 b1—a2 0 1 0 c2(5.56)—an_i 0 0 ... 1 Cn_i—a 00.. 0 b CnC=[i 0 ... 01,the system can be represented asdW1(t) = AW1(t)dt + .u(t — T)dt + .kdw(t)(5.57)y(t) = CW1(t)45Chapter 5: Continuous-Thne Variance AnalysL of Sampled fine-Delayed Systems5.3.2 Anti-Aliasing FilterIn a similar manner, a state space equation representing the inputloutput relationship of theanti-aliasing filter is expressed in controllable canonical form.Waa = GWaa + Hy(t)Ym = SWaa,—gi —g2 —gm—i —gm 11 0•• 0 0 0 (5.58)G= 0 1 0 0, H—0,00 0 1 0 0S = [s 82 8m—1 Sm]5.3.3 Combined State Space Description for Process and Anti-aliasing FilterNow, if we let,kA0 0 0A = , B = • , K = •LG0 010•..0(5.59)00.•.00000C=[1 0 •.. 0],Cm = [0 “. 0 s 82 83 Sm]the complete system can be represented bydX(t) = AX(t)dt + Bu(t — r)dt + Kdw(t)y(t) = CX(t) (5.60)ym(t) = CmX(t)This representation allows for easy calculation of the process variable y as well as the measured Ym.46Chapter 5: Continuous-Time Variance Analysis of Sampled lime-Delayed Systems5.4 Discrete-Time Model of Process and Anti-Aliasing FilterThe system is now sampled. Let the time delay ber=T3d+A (5.61)where 7’ is the sample period, d is the integer part of and A is the fractional delay.Using the multiplier e4t on (5.60) givese_AtdX(t)— e_AtAX(t)dt = e_AtBu(t— ‘r)dt + etKdw(t)(5.62)d(e_AtX(t)) e_AtBu(t— r)dt + etKdwIntegrating over the sample period givest+Te_A(t+T)X(t + 7’s) — e_A(t)X(t)= / e_A(8)Bu(s — r)ds + / e_A(3)Kdw(s)t-f-T8 t+TX(t + T3) = eATSX(t) +e4(t+Ts) / e_1’()Bu(s — ‘r)ds + et+Ts) / e_A(8)Kdw(s)(5.63)Noting that the control signal is held constant over the sample interval,t+AX(t + T3) = eATBX(t) + t+Ts) / e_A8dsBu(t — (d + 1)T3)t+T t+T (5.64)+e4(t+Te) / e_A8dsBu(t — dT3) +e4(t+Ta) / e_MKdw(s)Define the following variablese’’2’,= eA(t+T8) /e_48dsB = eA(Ts_A) edsB,t+Ta T8AFo = eA(2’) / e_AsdsB= / eA3dsB, (5.65)t+T‘I’(t + T3) =e4(t+Ts) / eA8 Kdw(s)47Chapter 5: Continuous-Thne Variance Analysis of Sampled Thne-Delayed SystemsThen equation (5.64) becomesX(t + T8) = X(t) + I’1u(t — (d+ 1)T3) + Tou(t — dT3) + W(t + T8). (5.66)Note that P is a stochastic process. Its statistics can be detennined as follows (See Aström [2]).Since dw has zero mean, I’ also has zero mean. The covariance ist+TR1 E[T] = E eA(t+T8_3)Kdw(s) / KTe4T(t+T8_r)dw(r)] =t+T.= / eA(t+T$_8)KRKTeAT(t+T8_8)ds (5.67)= f eAsKRKTeATSd8Using an integer index set to indicate the sample times, in which t = kT3, the discretized systembecomesXk+1 = Xk + Fluk_d_1 + rouk_d + Pk+1.Yk = CXk (5.68)Ym,k = CmXkThis discrete equation will allow us to calculate the state and control signal variance at thesample points. From these we can integrate the stochastic differential equation, (5.60), to obtain thevariance of the process and measurement between the sample points.The matrix exponential integrals in (5.65) can be calculated using methods given in [16].5.5 Estimation and PredictionAccording to the separation theorem, the tasks of control and estimator design can be handledindependently. Lennartson[8] has shown that restating the system so that the delay is shifted to theoutput rather than the input can lead to great simplifications in variance calculations. In essence, thisscheme moves the burden of dealing with dead-time from the control design to the estimator designwhere it is more easily handled. Lennartson introduces the ‘predictable’ stateX(t) = X(t + r) (5.69)48Chapter 5: Continuous-7Tme Variance Analysis of Sampled Thne-Delayed SystemsThen, 5.60 becomesdX(t) = AX(t)dt + Bu(t)dt + Kdw(t + r)y(t) CX(t — r) (5.70)Ym(t) = CmXp(t — r)Xp,k+1 = Xp,k + FtLk + k+d+1Uk = CX,k(5.71)Ym,k = CmXp,kZk = Ym,k + VkwhereF = r(T8) = feA3dsB (5.72)afld Ym,k is the sampled measurements, and vk is the measurement noise, assumed to be whiteGaussian with variance R2.The estimation of the current state is done with a Kalinan filterXk+lIk = .KkIk + rlUk_d_1 + rouk_dXk+lIk+1 Xk+lIk + Kk÷1 (Ymm,k+1— Cm±k+lIk)(5.73)Kk+1 = Pk+lCm[CmPk+lC + R2J’Pk+1 = [Pic — PkC [CmPkC + .R2] CmPk] T + R1where Pk is the covariance matrix of the estimation error Xk — Xkk. Prediction of X(kT3) =X(kT3+ r) is done by integrating (5.60) from kT8 to kT8 + r and taking the expectation givenX(kT3). Thus,t(kT3)=t(kT8+rI= + FIuk_d_1 + FU_1(5.74)where(r) = eAT (5.75)49Chapter 5: Continuous-Time Variance Aiwlysis of Sampled Time-Del.2yed Systems5.6 Controller DesignSince our task is to compare the best performance possible for controllers on a particular processover a wide range of sample rates, only a measure of performance based in continuous-time makessense. As the primary goal of most process controllers is to produce a product with the least variance,the time-average variance will be chosen for minimization. Since a stable controller will result incyclo-stationary statistics, i.e., the statistics of the controlled process are periodic, the cost may becalculated over a single sample period. Thus, the cost function is(k+1)Ta +TV (5.76)kT8+rwhere t kT3 is the current time. The square of the output is integrated from the beginning ofthe time of influence of the control action Uk until the beginning of the point of influence of thenext control action.5.6.1 Sampling the Continuous-Time Cost FunctionIt is convenient to sample the cost function so that once a sample rate has been chosen theanalysis may proceed entirely in the discrete—time domain. Equivalently, as in [2], the controllerminimizes the discrete-time performance indexV’ = XT(kTS + r)Q1X(kT5+ r) + 2XT(kT8+ r)Q12u(kT3+Q2u(kT3)= X’(kT8)Q1(3+ 2Xf(kT3)Q12u(kT8+Q2u(kT8) (5.77)= XkQ1Xp,k + 2X7’kQ12uk + Q2uwhere(k-I-1)T8Qi f (s)CCm(s)ds,kT9(k+1)T8Q12 J (s)CCmF(s)ds, (5.78)kT.(k+1)T8Q2= J [r(s)ccm (s)]dskT850Chapter 5: Continuous-Thne Variance Analysis of Sampled Thne-Delayed Systems5.6.2 Control LawIn this section we develop the control law corresponding to the cost function of the last section.The control law for the discrete-time system (5.71) and cost function (5.77) is now easily shown to beUk = —LX,k (5.79)whereL = (5.80)and iCp,kis the of the current predictable state as provided by the estimator.5.7 Variance CalculationsThis section outlines Lennartson’s [8] method for calculating variances at the output of theprocess and at the measurement.The steady state a priori error covariance matrixP = E{ (Xk XkIk_1) (Xk— xkIk_1)} (5.81)from 5.73 can be found by solving the Riccati equationp =— PC [CmPC +R2F’CmPI T (5.82)Let Xk = Xk— XkIk be the a posteriori estimation error for the measurement update and =X,k—Xp,kIk be the measurement update for the predictable state. And, define= E{XkXfl=Px = E{XP,kXlk}P = E{ (Xk XkIk) (Xk — XkIk)} (5.83)= E{ (xP,k—Xp,klk) (x,k — Xp,klk)}Px = E{Xu}P, = E{u2}51Chapter 5: Continuous-Thne Variance Analysis of Sampled Thne-Delayed SystemsThen, it can be shown for systems with optimal estimators thatP1 = P PC [CmPC + R2] ‘CmP (5.84)andPj = (r)P,jT(r) + R1(r) (5.85)Now the covariance matrix for the predictable state can be found by solving the Lyapunov equationP,ç, = (c— FL)Px(4 — FL)T +(5.86)— TL)Fç,( — TL)T + R1Because we are using a Kalman filter, X, and Xp are statistically orthogonal. Hence,=- P (5.87)This leads simply toPx,u = -PLPU=LPLT(5.88)= ((T3 — A) — r(T3—A)L)P ((T5—A) — r(T3 — A)L)T— A)P(T8— A)T + Ri(T8 — A)5.8 Time-Average VarianceHaving calculated the statistics of the system at the sample points, we may now proceed tocalculate the value of the cost function, namely the continuous-time average variance. This sectionpresents an algorithm for calculating the time average variance of the process output and the outputof the anti-aliasing filter.Lennartson has shown that the cost function, Jcan be calculated as follows.Jc = -trace[QiPx + 2Q + QPV] + Q2P (5.89)52Chapter 5: Continuous-lime Variance Analysis of Sampled lime-Delayed SystemswhereT9Qi= f4(s)Q(s)dsQ12 = JT(s)Qr(s)dsO2s (5.90)Q2 = JrT(s)Qr(s)czs= f3 — s)(s)RT(s)dsThe matrix determines the signal for which the cost is being calculated. For Q = CTC, the cost willbe that of the process output. This is the true process continuous-time variance. For Q = CCm,the cost will be that at the output of the anti-aliasing filter. This is the variance of the output ofthe anti-aliasing filter. Lennartson[81 has given accurate and efficient algorithms for calculating Qi,Q12 Q2, and P,5.9 The Calculation of Matrix Exponential IntegralsLennartson et al.[9] and Van Loan[16] have given algorithms for calculating the many matrixexponential integrals found in the formulation of the variance analysis. All of the integrals in theanalysis can be found in one of the following forms.F(s)=/8)Bds (5.91)R1()=(5.92)Q12()= 1TQcrds (5.93)= IrTsQcTsds (5.94)53Chapter 5: ContinuousThne Variance Analysis of Sampled lime-Delayed Systems= f ( — s)(s)QT(s)ds (5.95)Recall that (s) = eA8.Van Loan gives an algorithm for calculating matrices (5.91) through (5.94). In his formulationthe following matrix is formed.AT I 0 00 _AT Q 0C = (5.96)0 0 AB0 0 00The exponential of C is then taken and expressed in the following partition.F1() G1() H1Q) K1()0 F2(i) G2(L) H2()e = (597)0 0 F3(iS) G3(z)0 0 0 F4(i.)It is then shown thatF(s) =R1(I) = F1(z)G2(I),(5.98)Q12(z) = F1(I)H2(IQ() = [BTF1(I)KiQ)] + [BTFi(I)Ki(I)jIt remains to calculate (5.97). Van Loan gives an algorithm for this matrix exponential based onthe well known Padé approximation. This is the algorithm used by MathworksTht to make the expmfunction in Matlab, the function for calculating matrix exponentials. The Matlab expm function,however, was found to be insufficiently accurate for the variance analysis. It would therefore bepreferable to code the algorithm directly in order to have control over the error bounds. In this waywith experimentation accurate solutions could be achieved.Lennartson gives an algorithm for calculating (5.95). We have coded this algorithm and foundit to be accurate. A Matlab version of the algorithm, coded for this project follows54Chapter 5: Continuous-fine Variance Analysis of Sampled lime-Delayed Systemsfunction Qw = lennqw(A,Q,h)%function Qw = lennqw(A,Q,h)% LENNQW% This module calculates the integral%% h/I (h_s)exp(As)*Q*exp(As)ds/% 0% as shown in the paper lnter—sample Behaviour as% Measured by Continuous% Time Quadratic Criteria by B. Lennartson,% T. Soderstrom, and S. Zeng-Qi% mt. J. control, 1989, Vol. 49, No.6, 2077—2083.%%% INPUTS: h: a scalar, usually a sample period.% A: a square matrix.% Q: a square matrix of order(A).%% OUTPUT: Qw: a square matrix of order(A).%% AUTHOR: Scott Morgan% LAST REVISION: Dec 1, 1992%% FUNCTIONS CALLED:%%tol = le—lO;% Determine doubling factor.= 0;del =normA = norm(A*de1,fro);while normA > 5,j =del = hI(2j);normA = norm(A*del,fro);end% Summation.M = del*Q;N = 0;S =T = N;k = 1;while (norm(M,’fro’) > tol) + (norin(N,’fro’) > tol),N = del*k*M/(k+1);M = del*[A*M+M*A]/(k+1);S = S + M;T = T + N;k = k + 1;55Chapter 5: Continuous-lime Variance Analysis of Sampled lime-Delayed Systemsend% Doubling.for i = [j:—l:l],del = h/(2i);T = T + exp(A*del)*(T + del*S)*exp(A*del);S = S + exp(A*del)*S*exp(A*del);endQw = h*S -returnThe tolerance variable tol may be adjusted to determine the accuracy of the solution.5.10 Process and Measurement Variance for the Case of No Time DelayWe now present an example of variance calculations for a system with no time delay. Variancecalculations for the case of no time delay are more simple than for the time-delayed case and tend to bemore stable numerically. To simplify the analysis further we have set the measurement noise to zero.The plant is first order with transfer function(5.99)s+5The following first order noise shaping filer is used(5.100)s + 12A first order noise model is often used for moisture disturbances as indicated in Chapter 3. The poleis chosen to give disturbance energy at frequencies beyond the pole of the plant. This is consistentwith typical paper machine processes (see Bialkowski[3].)The process and measurement variance for the above system are plotted against sampling ratein figure 5.22. The variance at the measurement is always less than that of the process output sincethe anti-aliasing filter has removed high frequency variations. At low sampling rates the anti-aliasingfilter removes most of the disturbances. Thus the measurement variance is low. Because the controlsystem does not see those disturbances, it cannot reject them and so the process variance is high. Asthe sampling rate is increased the anti-aliasing filter passes higher frequencies. Thus the measurementvariance increases as more disturbances are being seen. The process variance decreases since the56Chapter 5: Continuous-Thne Variance Analysis of Sampled fine-Delayed Systemscontrol system can now reject the new disturbances it sees. As the sampling frequency is furtherincreased the process variance continues to fall as it sees more and more of the disturbance. But,because the controller is rejecting more of the disturbances, the measurement variance begins to fall.The disturbance energy falls off at higher frequencies; at higher sampling rates the measurement isthus reporting all of the disturbances to controller and so the process and measurement variancesconverge and continue to drop to zero as the sampling frequency goes to infinity.25201510510’Figure 5.22:10’ 10’ 10’ 10samp’e ra (radslsec)Process and measurement variance for a non-delayed plantIn an industrial control system one has only the measurements from which to gauge the qualityof the product. The plot of measurement variance for the above example is shaped as a hump. Themeasured variance is low for both high frequencies and low. Clearly the measurement variance isnot always a good measure of product quality as it can be seen in the figure that it is possible for theprocess to have high variance while the measurement system reports low variance.Taking high frequency MD estimates by the methods outlined in this thesis is one way todetennine how the variance reported by the control system compares with the variance as measuredusing higher frequency estimates. If the higher frequency estimates give higher measurement variancethan that reported by the control system then it is an indication that the measurement system is not— poc•ss—— measurement/ \/57Chapter 5: Continuous-Thne Variance Analysic of Sampled Time-Delayed Systemsgiving an accurate picture of product quality. One would have to sample at a frequency above the bandencompassing most of the disturbance energy in order to get an accurate estimate of process variance.There is an awareness among CD control designers that reducing the spatial sampling frequencyreduces the apparent output variance, a matter often alluded to in jaundiced comparisons of CD-controlvendor’s specifications. The awareness is likely due to the fact that with CD control hardware it is easyto experiment with sampling frequencies and filtering. It is not clear that there is a correspondingawareness for standard control problems.For systems with no time delay we expect that we may reduce the output variance arbitrarilyby choosing a sufficiently high sampling rate. This will not be the case for systems with time delaysince high frequency variations may come and go before the influence of control actions take effect.Thus for all time delayed systems there will be a variance, say VL, below which the process variancecan not be reduced no matter how high the sampling rate. In choosing a sample rate it would notbe worthwhile to sample at frequencies beyond the point at which the process variance approachesVL It is expected that for longer time delays the optimal sampling frequency decreases. The authorexpects that optimal sampling period, without considering the increase in control effort, would beon the order of 1/4 of the time delay.5.11 Some remarks on the difficulties encountered with the variance calculationsAlthough the variance calculations described in this chapter were implemented, they failed mostof the time. The reasons may be a combination of the following. The sampling process may beproducing unstable zeros that are being cancelled by the control law that minimizes the continuous-time average variance. Cancellation of these zeros would cause the control signal to be unboundedand result in the computations failing. For discrete-time minimum variance control laws this problemis solved by reflecting the unstable zeroes into the unit circle and placing poles at these locations(See Aström[2]). A similar approach may apply to the case of continuous-time minimum variancecontrol laws.The algorithms were implemented using the matrix computation software, Matlab. It isapparent that the accuracy of some of the prepackaged Matlab routines is not sufficient for our58Chapter 5: Conjinuous-7Tme Variance Analysis of Sampled 7Tme-Delayed Systemspurposes. In particular the exponential matrix function proved to be inaccurate. Lennartson[9] andVan Loan[16] have given accurate and efficient algorithms for the matrix exponential integrals. Someof these we have been successfully implemented in Matlab files. Coding the remaining integrals wouldease the difficulties of insufficient computational accuracy.The canonical forms used to express plant and filter dynamics are well known to give poorlyconditioned matrices. These perhaps should be transformed, if the formulation of the problem allowsit, into a better conditioned form before the calculations are carried out.Debugging the program was made difficult since cross-checking along the path of calculation wasnot possible. Thus it was difficult to verify that the steps were leading to correct or accurate answers.5.12 SummaryIn this thesis we have shown how higher frequency MD estimates can be achieved, either foruse in control or for monitoring the performance of the control system. The work described in thischapter is a step towards validating quantitatively the use of higher frequency estimates to reduceproduct variance. Distinct from previous analyses of this type we have considered the impact of theanti-aliasing ifiter on variance reduction. We have shown that this consideration is necessary becausefor low sample rates the anti-aliasing filter dominates closed-loop dynamics. An example has beengiven of variance analysis for a system with no time delay. It was seen that the variance may bereduced arbitrarily for sufficiently high sample rate and that for sample rates below the pole of thedisturbance, the measurement variance is a poor indication of control system performance.59Chapter 6: ConclusionsChapter 6Conclusions6.1 Summary of Results and ContributionsThe thesis describes how the EIBMC algorithm can be used to obtain MD estimates and how tofilter MD estimates in order to use the optimal sample rate for the given process. Variance analysiswas presented for determining the optimal sample rate given the process dynamics and time delay.A computer simulation of a scanning sensor has been developed for testing the effects of the filteringand data-shifting on the estimation of MD properties. It was shown• that sensors which convert their output to a frequency and then count over the sample periodin effect average the sheet property resulting in a boxcar moving average filter which may givepoor anti-aliasing.• that poor CD estimates result in poor MD estimates.• that CD estimation errors result in MD estimation errors in time proportional to the gradient ofthe spatial CD errors.• that anti-aliasing filters with non-linear phase response smear sheet features spatially to a degreedependent on their frequency distribution and that this smearing causes CD distortion and henceMD estimation errors.• that an anti-aliasing filter with optimally linear phase response, the Bessel filter, with data shiftingbefore estimation was the optimal way to filter the raw scanner data before estimation with theEIBMC algorithm.• how to choose the anti-aliasing filter order and cut-off frequency so that the delay becomes aninteger multiple of the sample frequency.• how to extract from the EIBMC algorithm MD estimates at any submultiple of the sensorsample rate.• how to formulate the variance analysis to allow for calculation of both the process variance andthe variance at the output of the anti-aliasing filter.60Chapter 6: Conclusions6.2 Suggestions for Further ResearchMany avenues of further research are suggested by the thesis, mostly stemming from the varianceanalysis. In particular, once the numerical difficulties of the variance calculations are solved, thetools thereby gained may be applied to the problem of determining the optimal sample rate for papermaking or other industrial processes. The tools could also be used to study the potential productquality improvement by altering the process by, for example, reducing the time delay by sensingbasis weight at the lip or at the dry line. Since control effort tends to increase with sample rate andvariance reduction, one must ensure that the control effort is tolerable when choosing a sample rate.It would be useful for the purposes of future process design to learn how control effort is related tothe minimum achievable continuous-time average variance for existing processes in working mills.For a time-delayed system, the number of delay lags in the sampled system increases with thesample rate. Since the time delay is a limiting factor affecting the minimum achievable continuous-time average variance, we can expect that we may determine a number of lags, say N, such that ifthe number of lags in a model is greater than N, then the system has been sampled at too high arate in the sense that little improvement in variance reduction is to be achieved by sampling at afrequency that gives more than N lags. It would be of interest to determine the number of lags, N,to be used as a heuristic in controller design.In this thesis we have applied the ELBMC algorithm to moisture and basis weight processes.However, the algorithm could well be applied to other processes such as caliper and the controlof paper additives. To apply the algorithm to a new process would require determining a suitabledynamic model for the process and modifying the Kalman filter part of the algorithm accordingly.This is an area requiring investigation.61References[1] K. Astrom. Computer control of a paper machine - an application of linear stochastic controltheory. IBM Journal, pages 389—405, 1967.[2] K. Astrom and B. ‘Wittenmark. Computer Controlled Systems. Prentice Hall, 2nd edition, 1990.[3] W. Bialkowski. Newsprint uniformity and process control - the untapped competitive edge. InCPPA Annual Meeting, Montreal, pages A255—269. Canadian Pulp and Paper Association,1990.[4] G. Dumont. Control techniques in the pulp and paper industry. Control and Dynamic Systems,37:65—114, 1990.[5] G. Dumont, M. Davies, K. Nataragjan, C. Lindeborg, F. Ordubadi, Y. Fu, K. Kristinsson, anda. Jonsson. An improved algorithm for estimating paper machine moisture profiles usingscanned data. In IEEE Conference on Decision and Control, 1991.[6] G. Havelock. Process Control. Pira Reviews of Pulp and Paper Technology. Pira International,1993.[7] I. Jonsson. Estimation and identification of moisure content in paper. Master’s thesis, Universityof British Columbia, 1991.[8] B. Lennartson. Sampled-data control for time-delayed plants. mt. J. Control, 49(5):1601—1614,1989.[9] T. Lennartson, B; Soderstrom. Intersample behaviour as measured by continuous-time quadraticcriteria. mt. J. Control, 49(6):2077—2083, 1989.[10] C. Lindeborg. A process model for the moisture variations at the dry end of a papar machine.Control Engineering Laboratory Report R 83-01, Chalmer University of Technology,Goteborg, Sweden, 1983.[11] C. Lindeborg. A process model of moisture variations. Pulp and Paper Canada, 87(4):T142—147, 1986.62[12] J. MacGregor. Optimal choice of the sampling interval for discrete process controL Technometrics, 18(2):151—60, 1976.[13] K. Natarajan, G. Dumont, and M. Davies. An algorithm for estimating cross and machinedirection moisture profiles for paper machines. In IFAC/IFORS Symposium, 1988.[14] V. Panuska. A new form of the extended kalman filter for parameter estimation in linearsystgems with correlated noise. IEEE Transactions on Automatic Control, AC-25:229—235,1980.[15] M. Salgado, G. Goodwin, and R. Middleton. Modified least squares algorithm incorporatingexponential reseting and forgetting. mt. J. Control, 47(2):477—491, 1988.[16] C. F. Van Loan. Computing integrals involving the matrix exponential. IEEE Transactions onAutomatic Control, AC-23(3):395—404, 1978.[17] G. Wang, G. Dumont, and M. Davies. Modeffing and identification of basis weight variationsin paper machines. IEEE Transactions on Control System Technologies, June 1993.[18] A. Williams and F. Taylor. Electronic Filter Design Handbook: LC, Active and Digital Filters.MacGraw-Hill, 2nd edition, 1988.63