Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Alignment, modeling and iterative adaptive robust control of cross-directional processes Farahmand, Fazel 2009

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

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2009_fall_farahmand_fazel.pdf [ 3.12MB ]
Metadata
JSON: 24-1.0070915.json
JSON-LD: 24-1.0070915-ld.json
RDF/XML (Pretty): 24-1.0070915-rdf.xml
RDF/JSON: 24-1.0070915-rdf.json
Turtle: 24-1.0070915-turtle.txt
N-Triples: 24-1.0070915-rdf-ntriples.txt
Original Record: 24-1.0070915-source.json
Full Text
24-1.0070915-fulltext.txt
Citation
24-1.0070915.ris

Full Text

ALIGNMENT, MODELING AND ITERATIVE ADAPTIVE ROBUST CONTROL OF CROSS-DIRECTIONAL PROCESSES by Fazel Farahmand B.Sc., Iran University of Science and Technology, 2001 M.Sc., University of Manitoba, 2005 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Electrical and Computer Engineering) The University of British Columbia (Vancouver) June, 2009 c© Fazel Farahmand 2009 Abstract A key issue in paper-machine Cross-Directional (CD) Control is alignment. Typi- cally, this mapping problem is a non-linear and slowly time-varying phenomenon for individual machines. The first part of this thesis specifically focuses on different causes of the mis- alignment and reviews recent developments in CD mapping. It summarizes the results of a comprehensive survey of different alignment schemes that have been used in CD control processes. This dissertation presents a novel, deterministic, tensor-based modeling of the closed loop controlled CD process and an alignment method. First, we link the CD data to the tensor model. Exploiting this link, we derive a deterministic blind PARAFAC decomposition as an alignment method. The proposed PARAFAC cap- italizes on the physical location of the actuators, scanning databoxes and their tem- poral diversities. Its performance is verified in several simulations then tested, evaluated and implemented as an alignment tool on real industrial paper machine. In the second section we present a new novel Adaptive Robust Control ap- proach to the multivariable CD process of continuous web manufacturing. We have applied discretized Internal Model Control(IMC)-based classical Windsurfing to each individual separated spatial frequency. This approach allows the dynamical bandwidth of the closed-loop system to be increased progressively at each spatial frequency through an iterative control relevant system identification and control de- sign procedure. The method deals with both model uncertainty and measurement noise issues. We also have applied Discrete-time H∞ Windsurfing to each individual sepa- rated spatial frequency, starting with an initial model and a robust stabilizing con- troller at each spatial frequency. This modified approach provides robust stability through iterations. The approach is evaluated through a number of simulation ex- ii Abstract periments. Finally, in a closed-loop approach to modifying the existing industrial CD con- troller, the objective is turned to the modification of the Windsurfing method. An input signal is designed for the system identification and control design. The ν- gap approach to robust control is fed into the Windsurfing model. Since the ν-gap metric establishes the upper and lower bounds, this approach reduces the number of iterations before the design is completed. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 The Industrial Paper Machine . . . . . . . . . . . . . . . . . . . 1 1.1.1 Sheet Properties . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Dynamic Limitations . . . . . . . . . . . . . . . . . . . 3 1.1.3 Spatial Limitations . . . . . . . . . . . . . . . . . . . . . 4 1.2 Objectives and Contributions of the Work . . . . . . . . . . . . . 6 1.2.1 Motivation and Objectives . . . . . . . . . . . . . . . . . 6 1.2.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1 Signal Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.1 Process Model . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.2 Cross Directional Response Model . . . . . . . . . . . . 11 2.2 Cross Directional Mapping and Misalignment . . . . . . . . . . . 12 iv Table of Contents 2.3 System Identification and Controller Design . . . . . . . . . . . . 13 2.3.1 Ill-conditioned Multivariable System . . . . . . . . . . . 13 2.3.2 Uncertainties . . . . . . . . . . . . . . . . . . . . . . . . 13 3 A Survey on Cross Directional Alignment . . . . . . . . . . . . . . . 15 3.1 Effects of Misalignment . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Shrinkage Modeling . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 Bump Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4 Image Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.5 On-Line Closed-Loop Non-Invasive Mapping Methods . . . . . . 31 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4 Tensor Modeling and Alignment . . . . . . . . . . . . . . . . . . . . 36 4.0.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.1 Tensor Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1.1 Correlation Cancellation . . . . . . . . . . . . . . . . . . 40 4.1.2 CD Tensor Model . . . . . . . . . . . . . . . . . . . . . 42 4.2 CD Alignment by Tensor Decomposition . . . . . . . . . . . . . 45 4.2.1 Alternating Least Squares Regression . . . . . . . . . . . 46 4.2.2 Identifiability and Uniqueness . . . . . . . . . . . . . . . 47 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.3.1 Simulation Results . . . . . . . . . . . . . . . . . . . . . 50 4.3.2 Experimental Results . . . . . . . . . . . . . . . . . . . 55 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5 Multivariable Windsurfer Adaptive Robust Control . . . . . . . . . 65 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.1.1 Cross-Directional Process Model . . . . . . . . . . . . . 68 5.2 Classical Windsurfer Adaptive Robust Control . . . . . . . . . . 69 5.2.1 Spatial Frequency Decomposition . . . . . . . . . . . . . 70 5.2.2 Controller Design . . . . . . . . . . . . . . . . . . . . . 73 5.2.3 System Identification . . . . . . . . . . . . . . . . . . . 75 5.2.4 Multivariable Windsurfing Algorithm . . . . . . . . . . . 80 5.2.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 v Table of Contents 5.3 MultivariableH∞ Windsurfing . . . . . . . . . . . . . . . . . . 89 5.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 92 5.3.2 Multivariable Iterative Adaptive Robust Control . . . . . 94 5.3.3 H∞ Controller Design . . . . . . . . . . . . . . . . . . . 95 5.3.4 System Identification . . . . . . . . . . . . . . . . . . . 99 5.3.5 Modified Multivariable Windsurfing Algorithm . . . . . . 101 5.3.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.3.7 IMC versusH∞-based Windsurfing . . . . . . . . . . . . 110 5.4 Cautious Bandwidth Control with ν-gap . . . . . . . . . . . . . . 112 5.4.1 ν-gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.4.2 Iterative Controller Design . . . . . . . . . . . . . . . . . 118 5.4.3 Enhanced Multivariable Windsurfing Algorithm . . . . . 119 5.4.4 Example . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.5 Input Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.5.1 Sinusoid on CD . . . . . . . . . . . . . . . . . . . . . . 122 5.5.2 Chirp Signal on MD . . . . . . . . . . . . . . . . . . . . 124 5.5.3 Example . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.7 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6 Summary and Future Work . . . . . . . . . . . . . . . . . . . . . . 129 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Appendices A Alignment Programs . . . . . . . . . . . . . . . . . . . . . . . . . . 155 B Windsurfing Programs . . . . . . . . . . . . . . . . . . . . . . . . . 156 B.1 IMC-based Windsurfing . . . . . . . . . . . . . . . . . . . . . . 156 B.2 H∞ Windsurfing . . . . . . . . . . . . . . . . . . . . . . . . . . 156 B.3 ν-gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 vi List of Tables 1.1 Physical Characteristics of Current Paper Machines. . . . . . . . . 5 4.1 Mapping Change Observed with PARAFAC Decomposition and Justified with the Bump Test. . . . . . . . . . . . . . . . . . . . . 59 vii List of Figures 1.1 Schematic view of a typical paper machine showing positions of the actuators and sensors. . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Three common CD response models. Where the thin line shows the Gaussian response, the thick solid line shows a response with negative side lobes and the shade curve corresponds to the Bimodal responses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Illustration of the scanner’s zig-zag measuring path. Black circles on each path correspond to the databoxes. . . . . . . . . . . . . . 6 2.1 The spatial interaction matrix G0 with a band diagonal structure (left) and under affection of linear shrinkage (right). . . . . . . . . 11 2.2 Feedback controller configuration. . . . . . . . . . . . . . . . . . 14 3.1 Typical symptom of mapping misalignment known as picketing. . 16 3.2 CD shrinkage profile across three different paper machines. . . . . 18 3.3 Linguistic variables of the Fuzzy Logic and their interdependence. 20 3.4 A block diagram illustrating the key working components embody- ing the discussed alignment method. . . . . . . . . . . . . . . . . 23 3.5 A block diagram illustrating key components of working embodi- ment of the alignment based on AFO. . . . . . . . . . . . . . . . 24 3.6 Experimental identification of the alignment of a paper machine based on bumping a few actuators. . . . . . . . . . . . . . . . . . 25 3.7 FFT spectrum of an image obtained from the full width sample viewed (Arkwork from Hoole et al. [92]). . . . . . . . . . . . . . 30 3.8 An excerpt of the time-frequency spectrum of the material-dependent part. (Artwork from Kaestner et al. [103]) . . . . . . . . . . . . . 31 viii List of Figures 3.9 Typical transmitted-light image by a CCD camera and used in es- timating the CD shrinkage profile of a paper machine. (Artwork from I’Anson et al. [98]) . . . . . . . . . . . . . . . . . . . . . . 32 3.10 2D FFT spectrum obtained of the image shown in Figure 3.9. (Art- work from I’Anson et al. [98]) . . . . . . . . . . . . . . . . . . . 32 4.1 Actual output of the actuator i. . . . . . . . . . . . . . . . . . . . 41 4.2 Schematic representation of the CD process concerning CD sheet shrinkage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3 Vector decomposition of the CD process corresponding to Equa- tion 4.12. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.4 Matrix decomposition of the CD tensor process at each instant of time, k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.5 Closed-loop response of actuators. . . . . . . . . . . . . . . . . . 51 4.6 The solid line shows the forced misalignment. The dash-dot line shows the original mapping neglecting the misalignment. The dot- ted line represents the databox aligned with the actuator at each moment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.7 The scanning databoxes of a system with no misalignment but with high model uncertainties. . . . . . . . . . . . . . . . . . . . . . . 53 4.8 The actuator mapping matches with the expected uniform distribu- tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.9 Identification of the CD response of an actuator using PARAFAC and ALS. The initial model has high level of model uncertainties while has no misalignment. . . . . . . . . . . . . . . . . . . . . . 54 4.10 An image of scanning databoxes with local instability resulted from both model uncertainties and misalignment. . . . . . . . . . . . . 55 4.11 Theapping obtained from PARAFAC shows obvious shift on the response of the last three actuators resulted from shrinkage on that area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.12 Identification of the CD response of an actuator using PARAFAC. The initial model has high level of model uncertainties while the system suffers from misalignment. . . . . . . . . . . . . . . . . . 56 ix List of Figures 4.13 Local instability (picketing) around databoxes 120 to 180. . . . . . 57 4.14 New alignment of the actuators based on the bump test. Black ellipsoid shows the location of the manual bump tests. . . . . . . . 58 4.15 Detected sudden alignment change of the actuator no. 24. . . . . . 58 4.16 PARAFAC results. . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.17 PARAFAC results. . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.18 PARAFAC results. . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.19 Alignment changes while performing bump test. The overshoot and undershoot a bump test responses of the actuator is not ob- served at the same databox. . . . . . . . . . . . . . . . . . . . . . 63 4.20 Identified matrix A. . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.1 Iterative adaptive control system. . . . . . . . . . . . . . . . . . . 70 5.2 Nonzero elements of a banded symmetric circulant spatial interac- tion matrix (left), the associated band-diagonal Toeplitz symmetric matrix (center), and the difference ‘ears’ (right). . . . . . . . . . . 72 5.3 Internal model control (IMC) structure. . . . . . . . . . . . . . . 74 5.4 Bandwidth increasing of the controller design. The entire shaded zone is the region Ωh. Its “boundary” is the curve where |G̃(ν, eiω)| = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.5 Identification of Rj,i and Ψj,i [118]. . . . . . . . . . . . . . . . . 78 5.6 Block diagram of the CD control applying multivariable Wind- surfer adaptive robust control. . . . . . . . . . . . . . . . . . . . 79 5.7 Iterative closed-loop output-error for model number i at the jth con- troller update. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.8 Comparison between the Singular Value Decomposition (SVD) of the toeplitz matrix and its circulant extension. . . . . . . . . . . . 84 5.9 Nominal model of the plant. . . . . . . . . . . . . . . . . . . . . 86 5.10 Iterative system identification on each spatial frequency. . . . . . 87 5.11 Response to a step reference input with CD square shape (n=20). . 88 5.12 Nominal and identified actuator response model of the basis weight. 90 5.13 Response to a step disturbance input with CD square shape. . . . . 91 5.14 Disturbance rejection of the system after 200 scans. . . . . . . . . 92 x List of Figures 5.15 Standard feedback configuration. . . . . . . . . . . . . . . . . . . 95 5.16 Standard sampled-data system. . . . . . . . . . . . . . . . . . . . 97 5.17 Discrete-time system. Kd is to be designed. . . . . . . . . . . . . 98 5.18 Identification of Rj,i and Sj,i [118]. . . . . . . . . . . . . . . . . 100 5.19 Iterative system identification on each spatial frequency. . . . . . 105 5.20 Response to a step reference input with CD square shape under noise free conditions (n=20). . . . . . . . . . . . . . . . . . . . . 107 5.21 Response to a step reference input with CD square shape corrupted with coloring noise (n=20). . . . . . . . . . . . . . . . . . . . . . 108 5.22 Rejection of the disturbance with CD picketing shape. . . . . . . . 109 5.23 Disturbance rejection of the noise-free system after 200 scans. . . 110 5.24 Identified model by Windsurfing based on IMC versusH∞ result. 111 5.25 Power Spectral Density of the identified model by Windsurfing based on IMC versusH∞. . . . . . . . . . . . . . . . . . . . . . 111 5.26 Normalized coprime factorization: C stabilizes P = M−1N . . . . 113 5.27 Projection onto the Riemann sphere of the Nyquist plots of two transfer functions and chordal distance between them. . . . . . . . 117 5.28 Identified model by enhanced Windsurfing based on ν-gap. . . . . 122 5.29 Enhanced Windsurfing response to a step reference input with CD square shape corrupted with coloring noise (n=20). . . . . . . . . 123 5.30 PSD for chirp and step inputs. While step input has clear excitation at zero frequency, the chirp excites the entire design bandwidth almost with equal power. . . . . . . . . . . . . . . . . . . . . . . 124 5.31 An example of chirp signal. . . . . . . . . . . . . . . . . . . . . . 125 5.32 CD model identification by chirp and step-input. While step input creates a model with corrupted corners, the chirp input generates a better response with higher resolution. . . . . . . . . . . . . . . . 126 5.33 Spectrum of the model identified by chirp input versus step input. 127 xi Acknowledgements This thesis is a result of a close collaboration between the industrial partner Can- for Pulp and Paper, Honeywell Industrial Solutions and the University of British Columbia. I would like to thanks my supervisor, Professor Guy Dumont for his guidance, advice, support, understanding and caring in difficult times throughout the course of this project. The efforts of my committee members, Professor Philip Loewen and Professor Michael Davies are gratefully acknowledged. I would like to thank members of control group at the Pulp and Paper Centre, Mohammed Ammar, Soroush Karimzadeh, Ali Soltanzadeh and Xinmei Shi for their valuable discussions and friendship. A number of people at the University of British Columbia have been very friendly and supportive. In particular, I would like to thank Lisa Hudson, Ken Wong and Tim Paterson for being there to help with any problem that I encountered. Further, I wish to thank Professor Rabab Ward for her support and kindness, during my early years in graduate school. Over the past years, I have benefited from the technical discussions and infor- mation supplied by a large number of people. I am grateful of Grant Emmond, Dr. Gregory Stewart and John Ball. xii Dedication I would like to dedicate this dissertation to my parents for their self-devotion, con- sistent guidance, stimulating and precious support, and most of all for their unwa- vering, patient love. I would never have reached my goals, if there had not been the steady support of my parents. Also this dissertation is dedicated to my beloved Farnaz who has been my source of motivation, inspiration and unconditional love. xiii Chapter 1 Introduction 1.1 The Industrial Paper Machine Industrial sheet and web-forming processes produce a wide array of products, ranging from plastic film, metal sheet and foil, fibreglass mats, multi-zone crys- tal growth furnaces to paper. A typical paper machine produces a sheet of paper 10 meters wide that emerges from the machine at speeds up to 120 km/h, under very strict quality specifications. Paper properties are measured by a scanning sensor at up to 2000 locations across the sheet and controlled by up to 300 independent actuators distributed across the width of the moving sheet. It typically takes about 20 seconds for the scanning sensor to traverse a 10m sheet (Figure 1.1). The overall machine can be divided into four sections [128], as follows: 1. The Sheet-forming section in the Wet end of the paper machine, where the mixture of water (99.5%) and pulp fibre is delivered onto the fast-moving forming fabric. The output of the wet end is a paper sheet with 20% fi- bre concentration. The sheet-forming section contains drainage elements of forming board, hydrofoil assemblies, table rolls, vacuum boxes, etc. 2. At the Press section, the sheet of paper is dewatered in heated steam-boxes and by counter-rotation rolls. The output of this section is a sheet that has approximately 40% fiber concentration. 3. The Dryer section of the paper machine consists of a series of large-diameter, rotating, steam-filled cans that evaporate the sheet’s moisture down to a 5-9% water content. Dryer sections have two typical configurations [148]: single- tier and two-tier. A single-tier dryer section is one in which the heated drums are arranged in a single-tier. Consequently they are long and require a large 1 Chapter 1. Introduction amount of mill area. A two-tier section generally requires less space and has an upper tier of spaced apart dryer drums and a lower tier of spaced apart dryer drums. 4. The Post-drying operations section is the final section of the paper ma- chine. In this section, the paper sheet thickness, or caliper, is controlled and paper properties are monitored by a moving scanner. As it emerges from the paper machine, the sheet of paper is wound onto a large reel. Figure 1.1: Schematic view of a typical paper machine showing positions of the actuators and sensors. There are two significantly different families of problems encountered in paper- machine control, corresponding to the two perpendicular axes most naturally as- sociated with the process: machine-directional (MD) control and cross-directional (CD) control. This work focuses on cross-directional control, where uniformity across the width of the sheet is the central objective. 1.1.1 Sheet Properties The most important paper properties that must be controlled are [50]: 2 Chapter 1. Introduction • Basis Weight: the sheet weight or mass per unit area in terms of grams per square meter (gsm). CD control of this property is accomplished by the slice lip dilution flow valves actuators at the headbox. • Moisture: the sheet moisture content, which is expressed as a percentage of sheet weight (%). The moisture of the paper can be corrected by steam boxes, infrared dryers and rewet shower. • Caliper: the thickness of the finished sheet, expressed in µm. The caliper of the paper is controlled by passing the sheet between rotating rollers in the calendar stack. • Coat: coat weight and thickness variations in coating processes. 1.1.2 Dynamic Limitations The open-loop dynamics in a CD control system are well-described by a first order system with dead-time. The transport delay between the actuators and the scanning sensor is called the dead-time Td. The closed-loop cut-off frequency is defined as the lowest frequency of the sensitivity function (i.e. 0dB point) by Haggman et al. [82]. Although this value depends on the type of applied algorithm, it is common to relate it to the closed-loop time constant λ. Note that λ is usually set around 1.5 times the natural process time constant, Tr. It is practically accepted [82, 166] that a typical CD control system attenuates disturbances slower than the cutoff period Tc: Tc ≈ k(Td + 1.5 Tr) (1.1) where k is a constant factor, chosen equal to 6. It is suggested by Bialkowski [11] that the reasonable value in which the con- troller removes 90% (−20dB) of effective disturbances at high-performance fre- quencies is: T90 ≈ 10 Tc (1.2) 3 Chapter 1. Introduction 1.1.3 Spatial Limitations The cutoff wavelength, Xc, of a practical CD control systems is determined only by the bump response shape. As shown in Figure 1.2, based on this shape the CD actuator response of different actuators can be divided into three groups. 200 400 600 800 1000 1200 1400 100 200 300 400 500 600 Figure 1.2: Three common CD response models. Where the thin line shows the Gaussian response, the thick solid line shows a response with negative side lobes and the shade curve corresponds to the Bimodal responses. GAUSSIAN RESPONSE The first class of CD response is considered for the consistency profiling for weight control [164] and rewet showers [187] and has a bell-shaped or Gaussian response with no side lobes (less than 15% of the response peak). NEGATIVE SIDE LOBES RESPONSE The spatial response of actuators with negative side lobes is shown in Figure 1.2. Slice lip actuators for paper grades of light to intermediate, steam box actuators 4 Chapter 1. Introduction that control moisture and induction heating of rolls for caliper control have negative side lobes response. Depending on the type of machine and the type of product this shape might be observed in a wide response [89] or narrow response [167]. BIMODAL RESPONSE The last but not least bump response has a bimodal or ‘two-humped’ appearance. Actuators with this shape of response are the most difficult to control [166]. The bimodal response shape is only observed in slice lip actuators of heavy grade paper, sack and board processes. The separation of peaks is likely due to the fact that in heavy grade paper product the machine speed is low. This speed allows propagation of waves in the pulp slurry before the weight profile is frozen at the dryline [66]. Physical properties of typical paper machines are listed in Table 1.1. Table 1.1: Physical Characteristics of Current Paper Machines. Paper machine Controlled Industrial Range of Number of Spacing Response Property Name Variation Actuators Width Weight Slice Lip 45-450gsm 50-118 7-20cm varies Dilution control 45-450gsm 60-80/150-200 6-7cm 7.5-12.5cm Moisture Steam Box 5-9% 55-171 7.5-15cm 30-60cm Rewet Showers 5-9% 50-120 7-15cm 15-18cm Caliper Calendar Stack 70-300µm 100-150 7.5cm 30-40cm Since the sheet of paper is moving much faster than the sensor, the scanner’s measurement pattern follows a zig-zag path on the sheet. As a result, scanned measurements are a combination of machine-direction and cross-direction mea- surement variables. 5 Chapter 1. Introduction Figure 1.3: Illustration of the scanner’s zig-zag measuring path. Black circles on each path correspond to the databoxes. 1.2 Objectives and Contributions of the Work 1.2.1 Motivation and Objectives The first objective of this thesis is to develop an automatic, on-line and non-intrusive cross-directional alignment method between CD actuators of a paper machine and their corresponding databoxes in order to improve paper quality. Most current methods require ‘bump tests’, in which a few actuators are excited, and the peaks in the observed scan data are assigned to the excited actuators. A major drawback of these methods is that they need to be manually initiated and thus require the CD control system to be in manual mode. Although there has been a fair amount of academic work for developing alignment methods, very few industrial applications have yet been reported. Currently there are no reported applications of adaptive robust control on CD processes. In addition there are some disadvantages with current methods includ- ing the lack of guaranteed robustness in presence of the inevitable model uncer- tainty, too conservative and/or sluggish controller designs and consequently satis- faction of low robust performance. The second part of this thesis aims at bridging the gap between existing theoretical work on adaptive control, robust control, iter- ative identification and control of CD processes. 6 Chapter 1. Introduction 1.2.2 Contributions The major contributions of this work can be highlighted as follows: • A complete survey of the industrial and academic state of the art on CD alignment of sheet making processes. • The development of a novel automatic nonintrusive online alignment method for CD processes. This method is applicable to any multivariable CD sys- tem with any number of actuators. It is incorporated into the industrial CD alignment. The results show that it perfectly maps the physical location of actuators to the scanning databoxes. • Application of the Windsurfing Iterative Adaptive Robust Control as a tool for the design of feedback controllers for paper machine cross-directional processes. The Windsurfing approach to adaptive robust control is shown to greatly improve the control performance while preserving robustness. • Extending the modified H∞-based Windsurfing iterative robust control to a multivariable process. Some other minor contributions are also summarized • The development of an enhanced Windsurfing strategy based on the ν-gap that does not require the static steps of classical Windsurfing. • Input design to the Windsurfing control of multivariable processes. 1.3 Thesis Overview An outline of the structure of this thesis is as follows: Chapter 2 describes the cross directional mapping problem. It also discusses CD controller design requirements of robust performance, robust stability as well as their constraints. Chapter 3 presents an overview of existing techniques for CD actuators’ map- ping. In addition the effects and impacts of misalignment on control of sheet pro- cesses are studied. 7 Chapter 1. Introduction Chapter 4 introduces a tensor model of sheet making processes. This tensor model is decomposed to find the CD model of actuators and CD shrinkage profile. The PARAFAC decomposition technique is used to decompose the large Multi- Input Multi-Output (MIMO) nonlinear time variant process into a fully decom- posed nonlinear shrinkage, CD mapping and response model in matrix format. Chapter 5 is devoted to the Windsurfing as a new approach to CD control de- sign. In the first section this chapter the classical Windsurfing approach to iterative adaptive robust control design is applied to the sheet making processes. Through iterations of adaptive robust control design and system identification on a simulated paper making process better controllers and plant models are obtained. The imple- mentation of the H∞-based Windsurfing for sheet making processes is the subject of the second section of this chapter. The Internal Model Control (IMC) design of the classical Windsurfing is modified and replaced with an H∞ (sub)optimal controller design. Two simulated examples based on industrial paper machine CD processes are given to demonstrate the effectiveness and efficiency of this method. Finally, in Section 5.4 a new improvement on Windsurfing is derived and evaluated through a number of simulation examples. The ν-gap is embedded into the H∞- based Windsurfing to reduce the number of iterations before a final multivariable controller is designed. An input signal is also designed to excite the actuators in closed-loop on MD and CD directions. A summary, conclusions and suggestions for future research directions are pre- sented in Chapter 6. 8 Chapter 2 Problem Statement In first section (Section 2.1) of this chapter, the industrial Cross-Directional process model is introduced. Then, the major difficulties of the mapping in these processes is described in Section 2.2. Finally, the requirements for this work are specified in Section 2.3. 2.1 Signal Model 2.1.1 Process Model The common nominal model for the weight process actuator dynamics was ad- dressed by Stewart in [164]. In the presented transfer matrix, it is assumed that the actuator responses can be represented in a square symmetric, band-diagonal Toeplitz matrix. The problem with this modeling is down-sampling of the scanned databoxes to the number of the actuators. It comes from the fact that a sheet of paper shrinks nonlinearly, while down-sampling is usually done in a linear fashion. Also, the model is not suitable for identifying misalignment, since shrinkage is not considered. Moreover, regarding the misalignment problem, the nonlinearity of shrinkage across the sheet [100], makes linear down-sampling inherently limited. The above-mentioned model was modified by Fan [52] to eliminate down- sampling. Based on this modification a new matrix of dimension I×R for the spatial interaction matrix was developed (Equation (2.4)), where I stands for the number of scanned data boxes and R represents the number of corresponding CD actuators. The associated MIMO model becomes: Y(z) = G(z)U(z) + D(z), (2.1) 9 Chapter 2. Problem Statement where G(z) = H(z) G0, (2.2) with H(z) = diag(hi(z))I×I , hi(z) = z−Td 1− aiz−1 , (2.3) G0 =  g1 gm · · · gq ... . . . . . . 0 . . . 0 0 g2 ... . . . ... 0 . . . . . . ... . . . ... ... ... g2 . . . ... gq . . . . . . 0 . . . ... ... gm g1 . . . g2 . . . . . . . . . gq . . . 0 ... ... g2 . . . g1 . . . . . . . . . ... . . . gq 0 gq ... . . . g2 . . . . . . . . . g2 . . . ... gq 0 gq . . . ... . . . . . . . . . g1 . . . g2 ... ... 0 . . . gq . . . . . . . . . g2 . . . g1 gm ... ... . . . 0 . . . . . . gq ... . . . g2 ... ... ... . . . ... . . . . . . 0 ... . . . ... g2 0 0 . . . 0 . . . . . . ... gq · · · gm g1  ∈ RI×R (2.4) and the measurement points Y (z) ∈ RI , the actuator set-points U(z) ∈ RR and the noise D(z) ∈ RI . h(z) is the Z-transform of the dynamic response of a first- order process with dead-time Td. G0 ∈ RI×R is the band-diagonal interaction matrix with the given parameters. These parameters, g1, . . . , gq, describe the spa- tial response of the process, where q corresponds to the half-width of an actuator’s response and m represents to the number of scanning spaces between adjacent ac- tuator centers. With this model, the problem of down-sampling is solved, but nonlinear shrink- age is still not considered. The difference in shape between the presented band diagonal interaction matrix and the realistic interaction matrix with shrinkage is illustrated in Figure 2.1. The interested reader is referred to [52, 164, 179] for further details on the 10 Chapter 2. Problem Statement 10 20 30 40 20 40 60 80 100 120 10 20 30 40 20 40 60 80 100 120 Figure 2.1: The spatial interaction matrix G0 with a band diagonal structure (left) and under affection of linear shrinkage (right). modeling and CD control of sheet and film processes. 2.1.2 Cross Directional Response Model The above-mentioned matrix is used in the profile model of Gorinevsky et al. [71]. In Cross Directional control, separability of the MD and CD models is commonly assumed. It is also assumed that all of the actuators’ CD response models are the same. In a typical paper machine different types of the actuators have different response curves [166]. In other words the shape of the response depends upon the type of actuator employed. Further, varying response shapes result from different grades of paper having being produced. The parameters that define the shape of the above cited model are w, the “width” of the actuator CD response, a, the size of response negative lobes or “attenuation” and δ the “divergence” parameter which stands for the presence of two maxima in the response and the distance between these two maxima. Currently these shapes can be derived/found only by forcing a “bump test” and performing 11 Chapter 2. Problem Statement system identification. High resolution examples of different models are shown in Figures 1.2. The most common spatial response shape for CD processes is the Mexican-hat with negative side lobes which is associated with slice-lip control of light to intermediate grades, steam box control of moisture and induction heating of rolls for caliper control. A bimodal response can observed on slice lip control of heavy grade paper, sack and board processes [166]. 2.2 Cross Directional Mapping and Misalignment There are quite a number of factors that significantly affect mapping between the measurements and the actuators in a CD process; some of them have been ad- dressed in [71, 122]. The factors cited by Gorinevsky et al. in [71] and Mast et al. in [122] include: 1. Initial geometric alignment of the CD actuators and the measurement device. 2. Position of the individual actuators within the CD actuator array. 3. Wandering of the paper web perpendicular to the machine direction. 4. Paper shrinkage through the drying process. 5. Drying variation across the web. 6. Machine speed changes. 7. Flow pattern of the extruded liquid paper stock in the initial stage of the paper-making process (paper-machine wire). 8. Grade changes. 9. Variation of sheet draw tension. Although the first two factors can be determined from accurate measurements, the complicated, time-varying and non-linear behaviour of the other factors creates significant problems [157, 188, 192]. 12 Chapter 2. Problem Statement The recent commercial introduction of full-width scans are the latest develop- ment in CD control [13]. These are expected to resolve many of the cross direc- tional control probleams [13, 105]. However, even with this type of sensors the CD mapping is still an unsolved issue. 2.3 System Identification and Controller Design 2.3.1 Ill-conditioned Multivariable System The interaction matrix of actuators in a cross-directional process is inherently ill- conditioned [60, 91, 111, 150]. It mostly happens if the spectrum of the actuator response shape is narrow and goes to zero within the Nyquist frequency range [78]. In another words the physical response of an actuator is wide. On the other hand if the physical response of the actuators is small, an accurate CD mapping and alignment becomes more difficult. Consequently a system would be inherently ill-conditioned and more sensitive to mapping. For example a small misalignment of just one actuator mis-mapping can destabilize a controlled but ill-conditioned sheet making system [48, 172]. Therefore it can be seen that no matter whether if the actuator response is wide or if it is narrow, they make an inherently ill-conditioned interaction matrix. 2.3.2 Uncertainties Different sources of model uncertainties can be listed as follows: • Alignment. • MD-CD cross coupling. • CD response shape. • Actuator failure. • Process dynamics. • Actuator nonlinearities. 13 Chapter 2. Problem Statement Currently no modeling strategy can capture all the above-mentioned features. Consequently all of the practical CD controller designs tolerate process model un- certainties [164]. The commonly accepted configuration for a closed-loop CD con- trol system is illustrated in Figure 2.2 Figure 2.2: Feedback controller configuration. The model uncertainty for such a system in the simplest form can be rep- resented in an additive unstructured perturbation on the nominal transfer matrix model of Equation 2.5. Πg := { G(z) + ∆G(z) : σ(∆G(ejω)) < l(ω) } ∀ω ∈ [−pi, pi] (2.5) where G(z) is the process model, ∆G is the additive model uncertainty, l(ω) is a positive scalar function which limits the set of all perturbed plants Gd(z) ∈ Πg to a neighborhood of the nominal model G(z) at dynamical frequencies ω and direc- tions m. Even for a simplified additive type of model uncertainty in a steady-state many singular values are vanishingly small. This implies that there is a singular value at a spatial frequency of m where σm(G(ejω)) < l(ω). In other words an additive model uncertainty of larger than the singular value of a specific spatial frequency is equivalent of %100 uncertainty in phase at that singular value. 14 Chapter 3 A Survey on Cross Directional Alignment In the first chapter we discussed the structure of paper machines. In this chapter we will discuss the state-of-the-art on mapping of a paper machine’s cross directional actuators into scanned databoxes. As the most important causes of misalignment is shrinkage, we will focus a major part of this chapter on the available methods of shrinkage’s identification. In this chapter CD mapping techniques are divided into the two main categories, open-loop and closed-loop. Intrusive and nonintru- sive alignment methods are other aspects of the CD mapping that we will discuss. Alignment techniques based on image processing are another possible approaches that we will also discuss. This chapter is divided into the following sections: As a starting point, Section 3.1 covers the effects and impacts of misalignment on CD control of paper ma- chines. Section 3.2 gives an insight into the physics of paper shrinkage. Recent contributions on the identification of shrinkage are also presented in this section. In Section 3.3 intrusive actuator-mapping techniques based on the bump test of ac- tuator(s) across the sheet are given. Image-based techniques of CD alignment are the subject of Section 3.4. Recent developments on closed-loop and non-invasive alignment methods are the subject of Section 3.5. This section also provides exper- imental results from a real-life example of misalignment and remapping. Some of the basic contributions of this thesis regarding nonintrusive alignment methods will be presented in this section. Section 3.6 is the conclusion with some final remarks and open issues. 15 Chapter 3. A Survey on Cross Directional Alignment 3.1 Effects of Misalignment The effects and impacts of misalignment on the control of sheet making processes have been studied independently by Fu [61] and Nuyan [136]. Fu et al. [61] have shown that the gradual degradation of a sheet’s profile is a function of a gradual increase in misalignment. We find that this phenomenon first appears at high fre- quencies and, at even higher levels of misalignment, instability results. It is closely related to the positive feedback resulting from uncertainties larger than 180-degree changes in phase of the actuator’s response [135]. This localized sinusoidal varia- tion, [128] or sawtooth pattern, [166] which appears in both the scanning databoxes and the actuator profile with roughly the same frequency and spatial location, is the most common result of CD misalignment. This pattern is often referred to as a “picket fence”, “actuator picketing” or “picketing” [17, 122, 126]. Figure 3.1 shows an industrial example of a typical picket fence around scanned databoxes be- tween 120 and 170 and actuator numbers between 19 and 29 observed on a modern Canadian paper machine. 0 50 100 150 200 250 300 350 400 450 500 550 −4 −2 0 2 4 CD Location Measured profile with mapping misalignment 0 10 20 30 40 50 60 70 0 20 40 60 80 100 Actuator Location Actuator setpoint with mapping misalignment Figure 3.1: Typical symptom of mapping misalignment known as picketing. 16 Chapter 3. A Survey on Cross Directional Alignment Although a well-tuned controller can compensate for most slight mapping mis- alignments [122], it is widely accepted that a mismapping larger than one-third of the spacing between CD actuators can cause instability of the CD control systems [48]. Fu et al. [61, 62] have shown that a mapping error of larger than half the spac- ing between actuators can cause instability of CD control systems. This threshold on the magnitude of mis-mapping errors in the actuator mapping has been exam- ined by Taylor et al. [172]. They have defined a shift mis-map function, adopted from Gendron et al. [64], by using the shift property of the Fourier transform. In their definition the phase of the Fourier transform specifies the CD location of each actuator. They have shown that the widely accepted mis-map threshold of one third of the spacing between CD actuators does not apply to all real-life control systems. By defining the full actuator spacing error as a new mis-map threshold, some im- provement in robustness has been achieved. However, a mapping error higher than this magnitude still causes a control system to destabilize and its associated CD actuators to eventually become saturated [172]. Currently the best CD control designs are based on a uniform distribution of actuators’ responses across the web without any misalignment or mapping change [52, 53, 146, 169]. Recent developments in sheet and film process control have been reviewed in a book by Featherstone et al. [60] and also in a paper by VanAntwerp et al. in [179]. Based on the literature review, it is apparent that online, nonintrusive alignment is still an active research topic in CD control. It can also be concluded that although there has been great progress in CD model identification and con- trol design, very few studies have been performed on closed-loop mapping and alignment. 3.2 Shrinkage Modeling Since a sheet of paper is under tension during most of the drying procedure, it is restrained from shrinking in the machine direction [100]. However, this ten- sion typically causes the paper to become narrower across its width and to have an anisotropic shrinkage [191]. Nonlinear cross-directional dimension changes affect nearly all paper machines [93]. As Figure 3.2 shows, this effect is significant, with shrinkage ranging from an extreme of 10% at the edges to close to 0% at the center 17 Chapter 3. A Survey on Cross Directional Alignment [100]. The resulting shrinkage profile is currently measurable only in the labora- tory or through manual tests. Shrinkage in paper machines causes misalignment of the sheet, also known in the industry as “mis-register”. I’Anson et at. [100] found that misalignment is extremely magnified in the newest, fastest machines(see Fig- ure 3.2). Figure 3.2: CD shrinkage profile across three different paper machines. It is well accepted that shrinkage affects various paper properties, e.g., dimen- sional, surface and compression properties [119, 131, 134]. Magnified variations in these properties can cause later problems in printing. A simulation of the CD shrinkage profile is based on the properties of the multiply board (e.g. the geome- tries of the free draws, where there is no CD forces acting on the edge of the sheet between each pair of dryer cylinders, in the dryer section [191]). The effect of free draws on the dimensional stability of a paper machine was demonstrated by Brecht and Wanka [16]. They obtained estimates of CD shrinkage profile by measuring hydroexpansion at a range of positions across paper machine. Later Gallay [63] used the same principle and demonstrated that the degree of paper expansion, when it is exposed to either vapor or liquid water is highly 18 Chapter 3. A Survey on Cross Directional Alignment correlated with the amount of shrinkage that occurs at the drying section. More direct estimates of the CD shrinkage profile were made by Wadhams et al. [188]. The dimensional changes of forming fabric marks during drying were measured in order to estimate the shrinkage profile. A few years later Viitaharju et al. [181] used a two-dimensional Fast Fourier Transform (FFT) and made measurements in the frequency domain to improve Wadhams results. Hoole et al. [93] adopted the method originally developed by Praast [143] and showed that free draws affect the CD shrinkage profile. Measure- ments of CD shrinkage on a commercial newsprint machine proved that the profile is not caused by the press section. Heaven et al. have developed an integrated tool for identification, alignment, shrinkage and dynamics of CD processes named IntelliMap [89, 90]. Based on the CD model suggested by Gorinevsky et al. in [75], height, width, attenuation and divergence are used to describe the actuator model. IntelliMap considers shrinkage as a linear phenomenon and uses actuator offsets to do the mapping as following: Shrinkage = Width at Actuator −Width at Scanner Width at Actuator (3.1) The alignment model is updated during the CD control system operation by us- ing real time sheet edge data from the scanner and from the wet end of the machine in order to compensate for sheet wandering [89, 90, 125]. Gorinevsky et al. in [73, 75, 76] have used ‘fuzzy logic’ to model shrinkage. They have assumed that nonlinear shrinkage has a known parametric form with a few unknown parameters. The input variable of the model is the CD coordinate x, with which three linguistic labels are associated: low edge, middle and high edge. Then the model is defined by the following simple inference rules: If [x is low edge] then s = slow If [x is middle] then s = s0 If [x is high edge] then s = shigh This method gives a parametric model to find nonlinear shrinkage [71]. In total six parameters should be found: one alignment parameter, three ‘shrinkages’ 19 Chapter 3. A Survey on Cross Directional Alignment parameters, and two ‘widths’ parameters on the edge nonlinear shrinkage zones. These parameters can be identified by means of an iterative minimization of the model-fit error. This approach, however, neglects the commonly-observed local instability caused by any local misalignment. This method also uses a bump test to find shrinkage parameters. At the same time Adamy [1] independently developed a fuzzy logic model for shrinkage. Compared to Gorinevsky’s model, Adamy’s proposed model considered more variables for shrinkage than just position alone (See Figure 3.3). Because of these large numbers of variables, fuzzy logic will be difficult to implement. There- fore decision logic is applied to a number of smaller, manipulable fuzzy subsys- tems. Consequently ‘fuzzification’ and ‘defuzzification’ are required at each level of Figure 3.3. Figure 3.3: Linguistic variables of the Fuzzy Logic and their interdependence. Wahlström et al. [189] use a finite element method to simulate in-plane stresses and strains in the paper during the free draw between the dryer cylinders, where no forces act on the edges of the web. Kniivila [109] also has demonstrated the effect of free draws on shrinkage and the MD tension profile. 20 Chapter 3. A Survey on Cross Directional Alignment Experimental curve-fitting is another approach to the identification of shrink- age in newsprint machines proposed by Philips et al. in [141]. They found that the best fit model is a reciprocal function, having three fitted parameters. This was done by dividing the CD profile into 2 halves and fitting each half to the following function: y = ax+ b x+ c (3.2) with a > 0 and b > c as the three fitted parameters. In 2003, Wahlström et al. [190] showed a steep CD-shrinkage gradient at the edges of the sheet of a two-tier paper machine. This behaviour was studied in more detail in a laboratory scale apparatus where the length, the width, and the anisotropy of the paper were varied. I’Anson et al. [101] have proposed another technique. In their method, the displacements of the headbox slice lip or dilution levels at positions across the ma- chine, along with the grammage profile, are used to estimate real-time CD shrink- age after calibration with a laboratory measurement of the profile. The core of their method is based on Equation 3.3, ∆W = 1 (1− S)(1 + E) , (3.3) where E stands for a fractional machine direction strain, the fractional CD shrink- age is represented by S and the change in grammage is shown by ∆W (gm−2). 3.3 Bump Test To perform alignment and detect nonlinear shrinkage, it used to be common in the industry to dye the paper at the CD actuators and detect its appearance on the scan- ner side [73, 76]. The dye test is a lengthy, cumbersome and imprecise procedure that results in the production of large amount of scrapped paper production and requires considerable manual labor [73, 76]. In the early 1990’s it was common to obtain the spatial response of a single actuator either offline, using a bump test, or online using a perturbation signal [49, 108]. The relationships among CD actuator settings, the grammage profile on the 21 Chapter 3. A Survey on Cross Directional Alignment wet end and the dry-line on a Fourdrinier table are discussed by Kjaer [107, 108]. Currently in addition to CD or 2D model identification, bump tests are used for mapping and alignment through the following three steps [88]: • Bump at least one of the CD actuators. The bump excitation signal might have different types, e.g. step-response test or, by using a pulse-sequence as the test input or, by utilizing a reception method [127]. • Measure a profile of the sheet at a substantial distance from the bumped actuators. In order to do that the signal is usually filtered and denoised. Also MD and CD components must be separated. • Determine the location of the effect of each of the bumped actuators on the sheet. Duncan has examined the 2D sheet variation measurement and based on Karhunen- Loeve (KL) expansion has shown that a 2D sheet variation can be decomposed into CD and MD variations [49]. In [49] he also has used this technique to estimate alignment. Metsälä et al. [126] have developed a mapping test. The test is performed by utilizing a reception method which employs correlated variance between the bumped actuators and the measurements.They also developed an alignment method through performing their mapping test [127]. They have assumed that the misalign- ment can be separated into a linear shift and nonlinear shrinkage. The mapping test result is used in the method to form a table of nonlinear shrinkage. Another aspect is that a linear mapping error is generated by measurement of the change in posi- tion of the paper web edges. Then, the shrinkage profile and linear error are used to develop a mapping model. A block diagram of their patented method is illustrated in Figure 3.4. Corscadden et al. [34] have used an Artificial Neural Network to estimate the effective locations of actuator responses in web forming processes as an open-loop. A Levenberg-Marquardt second-order nonlinear optimization is the method used to implement the artificial neural network estimator. This algorithm was subsequently improved by Corscadden et al. in [36]. By employing radial basis functions a suitable reduced-order estimator was defined 22 Chapter 3. A Survey on Cross Directional Alignment 200 400 600 800 1000 1200 1400 1600 100 200 300 400 500 600 700 800 900 1000 Figure 3.4: A block diagram illustrating the key working components embodying the discussed alignment method. in [36]. The basis functions in the hidden layer were described using discrete Chebyshev polynomials (DCPs). Then a recursive prediction error method was used to solve this nonlinear optimization problem. The above alignment algorithm has been implemented on a real machine [35]. Duncan’s method of CD modeling [49] was improved by Chen et al. [25]. A recursive algorithm was generated to update the 2D model. When a specified tol- erance for the convergence of estimates was achieved, iterations were terminated. It was also recommended to use wavelet decompositions [24] to extract the actual CD response, depending on the level of random variations in 2D sheet variations. The method is also tested on a paper machine. The idea of applying wavelet analysis to identify actuator responses in a CD process has been brought by Ghofraniha et al. [67]. To identify the location and shape of actuator response a matched filtering scheme or a continuous wavelet transform of the bump response with respect to the mother wavelet was found. Chen et al. have developed an automated optimization technique they called 23 Chapter 3. A Survey on Cross Directional Alignment the Adaptive Fuzzy Optimizer (AFO) method [28]. The method determines the locations of various mapping misalignments and establishes an effective perfor- mance indicator to measure the impact of this mapping misalignment. As shown in Figure 3.5, the method applies a searching technique, embodied in fuzzy logic, to search for and identify an improved CD mapping and then applies the improved CD mapping to the machine. Figure 3.5: A block diagram illustrating key components of working embodiment of the alignment based on AFO. Instead of bumping actuators Chen et al. [26] have used a Pseudo-Random Binary Sequence (PRBS) as a probing signal [22] for identification. Then the op- timal CD response model is found through a recursive algorithm. At each iteration the Nelder-Mead technique is used to perform the optimization. After defining some tolerances the iteration is terminated. In this method the common idea of Karhunen-Loeve (KL) expansion is used to decompose 2D sheet variations into temporal and spatial variations. Kalman filtering is used by Chen in [23] to channel and record the contribu- tions of CD and MD variability of sheet properties into the appropriate scanning 24 Chapter 3. A Survey on Cross Directional Alignment databoxes. Chen’s method provides improved estimates of the locations of CD variability. Meanwhile alignment in industrial CD control systems is based on an ‘identifi- cation experiment’ that involves bumping several CD actuators [66]. An open-loop bump test is performed while the control system is turned off in order to identify the spatial response of those actuators in this system [77]. Gorinevsky et al. have assumed that the shrinkage has a uniform model as follows: cj = α1 + α2xj (3.4) where cj is the center of the downstream response to the jth actuator; xj represents the position of the jth actuator and α1 and α2 are the parameters of the linear model. Gorinevsky’s method [77] is used to obtain the actuator response centers by filtering and/or averaging the profile data in the machine direction. An alignment function is then constructed by aligning the peaks of the responses with the corre- sponding actuators [72, 73] (see Figure 3.6). Figure 3.6: Experimental identification of the alignment of a paper machine based on bumping a few actuators. 25 Chapter 3. A Survey on Cross Directional Alignment An alternative approach to averaging the profile data in the MD is given by Gorinevsky et al. in [74, 76]. The algorithm is called Advanced Shrinkage Profile IDentification System (ASPIDS). Quadratic model fit error for the entire array of data is minimized to find the best model fit. Alignment parameters are identified by finding the best fit (maximal correlation) of the modeled system response versus the actual measured response. This technique provides the potential to perform the on-line monitoring and identification of paper alignment in closed-loop. A major drawback of existing techniques based on bump tests is that they need to be manually initiated and therefore require the CD control system to be in man- ual mode. Moreover a clear and reliable observation of bump centers requires large magnitude and the long term excitations of actuators. Finally bump tests are not able to adapt to frequent and continuous mapping variations [122]. In a recent study of shrinkage, I’Anson et al. [101] present a direct relation- ship between fractional CD shrinkage and paper grammage. The shrinkage cor- responding to each actuator is calculated by assuming a constant paper weight across the paper’s cross-directional profile. I’Anson et al. [101] assume that the shrinkage profile is broadly symmetrical about the central point of the dry paper. However, misalignment is commonly observed that does not follow the assumption of symmetrical shrinkage. Consequently, this assumption might not be appropri- ate. Moreover, I’Anson presumes a constant weight across the entire sheet, which is not the case when there is a misalignment. Constantino et al. have proposed an important model of shrinkage in [32]. They suggest that CD shrinkage at a certain point in the paper sheet depends on the dis- tance of that point from both edges of the sheet, and also depends on the length of the unsupported draws in the dryer section. The model presented provides a reason- able mapping that fits the experimental data. It also considers a constant factor for MD variations. Although this model is useful for the modeling of shrinkage, its de- pendence on so many unknown factors, that are also a function of grade/conditions of papermaking, makes it inappropriate for alignment identification. This method also suffers from the same assumption of symmetrical shrinkage. Taylor et al. [173] have used Bayesian methods for CD identification by an- alyzing the results of a bump test. The Bayesian approach uses a posterior dis- tribution and considers non-linear parameters of the model, for example actuator 26 Chapter 3. A Survey on Cross Directional Alignment mapping. The difficulty with Bayesian analysis is the requirement to solve certain intractable integrals that define the likelihood ratios shown in Equation 3.5. This issue has been resolved by applying the Gibbs sampling and Metropolis Hastings (MH) sampling forms of the Markov Chain Monte Carlo (MCMC) methods [145]. f(θ|y) = f(θ)f(y|θ)∫ f(θ)f(y|θ)dθ , (3.5) where f(θ|y) is the likelihood function of the conditional posterior probability distribution, θ stands for the parameter(s) of the model, y is the given data and the prior distribution, f(θ). Recently, Stewart [165] has patented a closed-loop method of alignment called the “Reverse Bump Test”. In this method, the closed-loop system is excited with a probing/perturbation signal to the scanner measurement rather than the actuator setpoints. Alignment is studied in steady state conditions and the phase change of the actuator profile is used to drive the required alignment update. As an assump- tion it is an essential that at low temporal and spatial frequencies, Equation 3.6 will be satisfied. [1 +KG]−1K ≈ G−1. (3.6) They have concluded that the feedback assists identification rather than biasing it. They also concluded that at low spatial frequencies the phase of the actuator profile response contains exactly the alignment change expected at that location of the paper sheet. Finally, the capability of the conventional bump tests is heavily limited to the scanning sensor. The scanning sensors measure the extremely sparse data of the sheet properties profile every 20-60 seconds. In addition the scanned data are highly corrupted with noise. Consequently the conventional bump tests have to be extended for multiple scans before one can observe a reliable actuator response. A bump test not only suspends the normal control operation, but also causes more severe product deviations for an extended period of time. 27 Chapter 3. A Survey on Cross Directional Alignment 3.4 Image Analysis Optical measurement of paper’s characteristics dates back to 1976 and has been described in previous reviews and textbooks [12, 15, 37, 68, 96, 97, 138, 142, 153– 155]. During manufacture, paper is in contact with different elements that leave a faint impression on it. Among these elements, the most important elements are the forming fabric, the dandy roll (if applicable), felts, presses and drying cylinders [80]. The imprints made by the forming fabric are called wire marks. Periodic measurement of wire marks in the sheet can be used to determine the shrinkage profile. In an image processing method pioneered by I’Anson et al. [99] those imprints are used to estimate the shrinkage profile. In their method frames of a microscopic image of paper are repeatedly taken. It is important to have these images being taken with illumination at a grazing or low angle of incidence. Then the images are superimposed on one another, with a slight offset. When the offset distance just matches that of the repeating pattern in the forming fabric, the “wire mark” becomes sharply apparent. At the same time, any non-repeating features such as the effects of fiber flocs, tend to cancel each other out, when averaging the repeated images. Di Mauro et al. [123] have offered four different methods that use wire marks to find the shrinkage profile. Their first three algorithms model the paper machine as a stationary process, while their fourth algorithm, instead, assumes the paper signal to be non-stationary. Again the input signals are collected by line-scan CCD cameras. Specifically the algorithms are as following [123]: 1. Fourier Analysis: This identifies the peaks in the frequency domain due to the repetitive patterns. The position of the main peak in the 2D Fourier spectrum corresponds to the frequency of the main knuckle separation of the mark in the line-scan image. 2. Autoregressive (AR) model: This structures the process model without a priori information by a recursive least-square estimation. 3. Eigen-analysis methods: This extracts information about the frequency of the paper marks, and at the same time, estimates the variance of noise in the 28 Chapter 3. A Survey on Cross Directional Alignment image. 4. Phase-Locked Loop: The pattern differential shrinkage is seen as a frequency modulation process. The unshrunken pattern is the carrier and the shrinkage profile is the modulating signal. Niemi et al. [132] have proposed another approach by using a high-speed video camera to monitor the dryline on the machine. This technique provides a general indication of basis weight but is also influenced by drainage characteristics on the wire and limited to Fourdrinier type machines [27]. The wire marks are often recorded by measuring the intensity of light transmit- ted through the paper, using a CCD camera [93]. The wire distance can then be de- termined by analyzing the calculated image of the wire pattern. A two-dimensional Fast Fourier Transform (FFT) on images at a series of positions was conducted by Guesalaga in [80, 81] and later by Hoole et al. in [93]. Hoole et al. [92] have named the peaks, which are very close to the vertical axis of the FFT spectrum, the ‘distortion peaks’. They are justified to be produced by the mark of the CD yarns in the wire. They also have assigned the horizontal axis peaks to the “shrinkage peaks”, which have been used to measure the relative CD separation of an aspect of the wire mark close to the MD (See Figure 3.7). Guesalaga et al. [79] have also presented another image processing approach to identification of shrinkage profile. They consider the effect of the forming fab- ric as ‘water marks’ or ‘signature’ on the image, where the frequency spikes are originated by these marks. These spikes, also shown in Figure 3.10 and numbered 1-4, tend to move in the spectrum plane causing a shift to higher frequencies as the paper shrinks during the drying stages. Kaestner et al. [103] also have used a FFT image processing technique with a one dimensional red fluorescence sensor [133] recording (shown in Figure 3.8). They have named the method Fluorescence-based Shrinkage Measurement (FSM). Three different types of frequency estimators based-on time-frequency spectra have been evaluated for estimating the shrinkage profile of newsprint. The method de- tects the variations in lignin content in the forming-fabric mark. These variations do not appear to be and should not be counted as forming-fabric distortion. The profile estimation is divided into the following subtasks: Preprocessing to extract 29 Chapter 3. A Survey on Cross Directional Alignment 50 100 150 200 250 300 50 100 150 200 250 300 Figure 3.7: FFT spectrum of an image obtained from the full width sample viewed (Arkwork from Hoole et al. [92]). the region of interest from the full spectrum, Frequency estimation to estimate the location of the web-related harmonics and Profile estimation to fit a polynomial to the harmonics. Three different estimators are evaluated, as follows: • Locating the harmonic with maximum intensity (MAX). • Thresholding based on a Neyman-Pearson criterion (NP). • Correlation-Energy Maximization (CEM). Villforth et al. [182] used fuzzy clustering methods to accurately find the loca- tion of the peaks in the FFT spectrum without operator intervention. The above-mentioned method in [101] recently has been refined by I’Anson et al. [98] for laboratory estimation of CD shrinkage profile. As in the established method this method uses FFT to detect the dimensions and orientation of the fabric marks. It also uses a geometrical method to produce the profile. Figure 3.10 shows the 2D FFT spectrum of the paper image shown in Figure 3.9 with the highlighted key peaks. I’Anson et al. [98] have used the bottom 30 Chapter 3. A Survey on Cross Directional Alignment Figure 3.8: An excerpt of the time-frequency spectrum of the material-dependent part. (Artwork from Kaestner et al. [103]) forming-fabric peaks (i.e. 1 and 2) to find the model parameters, while peaks 3 and 4 are the corresponding maxima for the top forming-fabric marks [98]. 3.5 On-Line Closed-Loop Non-Invasive Mapping Methods State-of-the-art in the identification of sheet misalignment is based on the detection of picketing by the human eye or by defining some threshold for the maximum variation of the sheet’s properties. Unfortunately it is not practical for a human to observe web variations on an uninterrupted basis [176]. Moreover, detection of misalignment with current methods is possible only after misalignment problems have persisted for a significant period of time with a significant value [176]. For a number of years, researchers have been exploring methods to find CD mapping based either on open-loop procedure or through intrusive methods. The next couple of sections will discuss a number of these new, non-intrusive alignment 31 Chapter 3. A Survey on Cross Directional Alignment Figure 3.9: Typical transmitted-light image by a CCD camera and used in estimat- ing the CD shrinkage profile of a paper machine. (Artwork from I’Anson et al. [98]) Figure 3.10: 2D FFT spectrum obtained of the image shown in Figure 3.9. (Art- work from I’Anson et al. [98]) 32 Chapter 3. A Survey on Cross Directional Alignment methods. Saffar et al. [151] have proposed an approach to develop indirect closed-loop identification in an unconstrained model-predictive-control (MPC) framework for the CD processes. The MPC technique developed originally by Phan et al [140] uses a state-space system and controller, in a combined two-step process. An ob- server is used within the identification equations. Then the closed-loop observer Markov parameters are identified and converted into a step-response model. The results of the simulation confirm the validity of this method in identifying a sim- ulated, nonlinear shrinkage of up to 10% at the edges of the paper sheet. This method also assists in finding the correct CD alignment. Karimzadeh has proposed an online picketing detection scheme based on the monitoring of the frequency components derived from actuator profiles [104]. The idea is based on the fact that the occurrence of the phenomena of picketing results in the growth of the high frequency components of the actuator profile. The change in the power of this component is then detected by a Cumulative Sum (CUSUM) algorithm. A newly proposed, online technique that is non-intrusive has been suggested by Gopaluni et al. in [69, 70]. The proposed technique theoretically is another approach to the nonparametric system identification technique given in [121]. It is based on estimating the cross-correlation between actuator movement and the measurement points to find a mapping. To find this correlation, the reference sig- nal, including all of the inputs and outputs should be passed through a whitening filter. The peak of the correlation is then chosen as the mapping, independent of the shape of correlation function. We see in Figure 3.1 a phenomenon called ‘picketing’. Caused by misalign- ment, it commonly appears as the saturation of adjacent actuators driven to their minimum and maximum values. During the misalignment, feedback from a mis- mapped actuator is chosen from a misaligned scanning databox. Hence, the mis- mapped actuator will also have high correlation with wrong databox. In such a case cross-correlation will not give a unique peak for this actuator located at the area of the misalignment and in fact, quite a number of peaks might appear in a cross-correlation function. Stationarity of the signals is the main requirement to apply this correlation method [121]; however, if there is a misalignment we then 33 Chapter 3. A Survey on Cross Directional Alignment see an increase in both the mean and the variance of the databox measurements and the actuator profile. A new closed-loop identification method is presented in [122], based on the newly defined “Local Variability Index” denoted as “ρ(k,w)”, where k stands for the CD location and w is the local window size. Instead of using the global in- dex, 2σ, that measures the variability across the entire sheet; the maximum local variability is used for quantification of the mapping misalignment at that location. 2σ = 2 ( n∑ i=1 (p(i)− p)2 n ) 1 2 (3.7) where p(i) 1 6 i 6 n is a profile, p is profile average and n is the size of the profile. ρ(k,w) = 1 σ ( k+w∑ i=k−w (p(i)− pk)2 2w + 1 ) 1 2 (3.8) where pk is the local profile average over the specified window. It is claimed that this maximum always occurs because of local instability caused by misalignment. This method is currently used as the core of an indus- trial control system [129]. However, in addition to misalignment there are several other factors that significantly effect the variability of the paper sheet [71]. 3.6 Summary A conceptual introduction to the effect of misalignment on CD control followed by an overview of existing alignment or mapping analysis techniques has been presented in this chapter. As the main reason of alignment change, some physical properties of the sheet shrinkage have been explained. Depending on the type of excitation and control, the alignment methods were classified into two different groups: open-loop or closed-loop. Then, depending on the presence or lack of external excitation for mapping, the methods were classified into intrusive or non- intrusive groups. In those cases of alignment where a bump signal is employed as a reference 34 Chapter 3. A Survey on Cross Directional Alignment input to the closed-loop system, or as a direct excitation of the actuators, the corre- sponding method is categorized as a bump test method. They are performed either in open-loop or under closed-loop CD control. Vision-based techniques is another approach to the alignment and shrinkage identification that has been studied in this chapter. In vision based methods, image processing techniques are used to detect wire mark or a machine’s physical effects on the paper, compared to the other methods that consider just physical properties of paper with respect to the actuators. The last group of alignment methods is the non-interrupting or non-intrusive techniques. In this class of mapping techniques, the alignment method won’t in- terrupt the closed-loop system and also in order to do the alignment no external excitation signal is needed. 35 Chapter 4 Tensor Modeling and Alignment A key issue in paper-machine Cross-Directional (CD) Control is alignment, i.e., accurate spatial mapping of the actuators to the scanning points. Typically, this mapping problem is non-linear phenomenon that slowly varies over time. Most current methods require bump tests, in which a few actuators are excited, and the peaks in the observed scan data are assigned to the excited actuators. A major drawback of these methods is that they need to be manually initiated and thus require the CD control system to be in manual mode. The problem statement of Chapters 2 and 3 requires that a nonintrusive and efficient CD mapping and alignment method be designed for a large-scale CD pro- cess to deal with these nonlinear time-varying mapping changes. The main goal of this chapter is to present a novel, deterministic, tensor-based model of the CD process and secondly to offer an alignment method that continues to function while the closed-loop CD control system is running. The motivation for using wireless communication techniques for mapping between scanning databoxes to their cor- responding actuators comes from a well-known decomposition technique. First, we link the CD data to the parallel factor model (PARAFAC). Exploit- ing this link, we then derive a deterministic blind PARAFAC decomposition as an alignment method with performance close to non-blind minimum mean-square er- ror (MMSE). We show that a blind alignment follows from simultaneous matrix decomposition. The proposed PARAFAC capitalizes on the physical location of the actuators, the scanning databoxes and their temporal diversities. Its perfor- mance is then verified in several simulations for different actuator models. The discussed algorithm is then tested on industrial paper machine data and evaluated as an identification tool. In previous chapters we discussed the structure of paper machines and the state- of-the-art in CD identification and alignment of the actuators. This chapter will 36 Chapter 4. Tensor Modeling and Alignment introduce a novel deterministic alignment method. The work is divided into the following sections: In section 4.1 we develop a new multi-linear algebraic model. The PARAFAC decomposition (Parallel Factor) or CANDECOMP (Canonical De- composition) developed independently by Carroll et al. [20] and Harshman [85] is used as an alignment tool in Section 4.2. It is done by extracting the sequences generated by the CD actuators. Section 4.3 illustrates the performance of the pro- posed technique by means of simulations. This section also provides experimental results in a real-life example of misalignment. Section 4.4 is the conclusion of this chapter. 4.0.1 Notation Let us first define some of the notations that we will use. A calligraphic letter Y denotes a third-order tensor. A bold-face capital denotes a matrix (Y). Vectors are written in italic capitals (Y) and Yk indicates the kth column of matrix Y. Italic capitals are also used to denote index upper bounds (i=1,2,...,I). Scalars are iden- tified by lower-case letters (y). The scalar yi indicates element i of vector Y, the scalar yi,j denotes the element of the ith row and the jth column of matrix Y and the scalar Yi,j,k denotes the ith element of the first dimension, the jth element of the second dimension and the kth element of the third dimension of tensor Y . The norm ‖ · ‖ is the Frobenius norm. Furthermore, the operator vec(·) builds a vector from a matrix by stacking the columns of this matrix one above the other; more specifically, element yi,j of the (I×J) matrix Y becomes the element at posi- tion i+(j-1)I of the vector vec(Y). The symbol⊗ represents the Kronecker product, A⊗H def=  a11H a12H · · · a21H a22H · · · ... ...  . The symbol  represents the Khatri-Rao or column-wise Kronecker product: AH def= ( A1 ⊗ H1 A2 ⊗ H2 · · · ) . 37 Chapter 4. Tensor Modeling and Alignment The logo ∗ symbolizes the convolution function, (f ∗ g)(m) = ∑ n f(n)g(m− n), and ◦ denotes the outer product, defined by (U ◦ V)i,j def= uivj . The notation can be summarized as follows: j preshrinkage CD location, 1 6 j 6 J r actuator number, 1 6 r 6 R k time instant / scanning number, 1 6 k 6 K i databox number, 1 6 i 6 I ur,k reference input of actuator r at time k cr,j CD response gain of actuator r at location j sr,k temporal response of actuator r at time k hk system’s dynamic (system’s impulse response at time k) ar,i fading factor of actuator r to scanned databox i xr,j,k paper property before shrinkage yi,j,k paper property measurement in databox i, time k, due to preshrinkage location j vi,j,k measurement noise 4.1 Tensor Modeling Current models of paper machines consider only the number of scanning databoxes as the only CD dimension of the sheet. All of them ignore the initial CD dimension of the paper’s sheet at the actuators location−i.e. before shrinkage happens. This simplification has made online mapping an unsolvable problem. Having higher order of dimensions (i.e. time, initial width of the sheet, CD databoxes at the scanner and actuators) suggests developing a tensor model for the papermaking processes. 38 Chapter 4. Tensor Modeling and Alignment Moreover, the interaction matrix of the actuators in paper processes are inher- ently ill-conditioned and their input signals are highly dependent. Therefore the CD response of the actuators produces a sparse, non-full-rank matrix [164]. Con- sequently classical Least Square (LS) methods of system identifications cannot be applied to decompose a multivariable system. The existing approaches to array signal processing and system identification mostly rely on matrix (2D array) models. Usually only space and time dimensions are considered. Array signal processing is generally used in wireless communi- cation systems to mitigate multiuser (co-channel) interference, inter-symbol inter- ference as well as to benefit the spatial diversity available in the wireless channel. However, the inherent lack of uniqueness in working with a matrix model in sig- nal recovery is the main limitation [40]. Current blind algorithms generally use special and specific properties of the data such as orthogonality, constant-modulus, finite-alphabet and/or cyclostationarity to overcome the non-uniqueness of matrix decompositions and successfully perform multiuser signal separation equalization [139, 171, 177]. One of the most studied decompositions of 3D, and higher dimensional ten- sors, is called PARAFAC. It is a subset of a multi-way analysis that can be viewed as multi-linear algebra. PARAFAC relies on a fundamental result of Kruskal [112] regarding the uniqueness of low-rank three-way array decomposition to simulta- neously recover all information-bearing signals and associated system parameters only from output data [160]. As a result PARAFAC does not require even the spa- tial model of actuators. PARAFAC can be seen as generalizing joint-approximate diagonalization ideas [178]. Compared to the decomposition of a matrix in rank-1 terms that is not unique, the power of PARAFAC stems from its unique properties. In matrix decomposition, orthogonality of the components, which is unique up to trivial indeterminacies, is usually imposed. The singular value decomposition (SVD) and the symmetric eigenvalue de- composition (EVD) can be seen as ways to decompose a matrix into mutually- orthogonal rank-1 terms. However, orthogonality conditions may not be satisfied by the actual underlying components that arise in paper making models. Therefore those decompositions do not provide useful insights in our case. 39 Chapter 4. Tensor Modeling and Alignment Independent Component Analysis (ICA) and Principal Component Analysis (PCA) are very popular but offer non-interpretable results when a higher-dimensional structure is hidden in the data. PARAFAC is most often used as a powerful tool in wireless communication to decode transmitted signals. It is a deterministic technique and does not require the estimation of statistics. PARAFAC is also blind and needs no training sequences. In this thesis, a new tensor-based model is developed for the actuator response of the paper making processes. Then considering the unique power of the PARAFAC decomposition, this technique is used to decompose our developed three-dimensional tensor− aiming to recover an online mapping (i.e. alignment of CD actuators). We think that this method is appropriate when compared to wireless applications since both the input and the scanned output signal are known. Further, this method can also be used in an alternative fashion to identify the actuator model. 4.1.1 Correlation Cancellation A mapping method based on correlation cancellation between input noise and the actuator set-point has been proposed by Farahmand et al. [54]. The measurement noise satisfies the Gauss-Markov assumption [170]. Then the Best Linear Unbiased Estimator (BLUE) will be found. The parameters of this estimation provide the actuator model and its mapping. This method is a batch algorithm, which means a rather large number of scans are needed to apply this technique. The Second Order Statistics (SOS) was then used to do the mapping, based upon the technique of Blind Source Separation (BSS). This technique is the basis of spatial, temporal and spatio-temporal decorrelation methods employed in sig- nal processing. Spatial decorrelation, or prewhitening on FIR systems, (also called correlation canceling), is often considered as a necessary (but not sufficient) condi- tion for a stochastic independence criteria. This spatio-temporal and time-delayed decorrelation can be used to identify the mixing matrix and perform BSS of colored sources [31]. Rewriting Equation (2.1) for the closed-loop model shown in Figure 4.1 yields 40 Chapter 4. Tensor Modeling and Alignment the following: Y (k) = GC 1 +GC R(k) + Dt 1 +GC W (k) (4.1)   100 200 300 400 500 600 700 50 100 150 200 250 300 350 Figure 4.1: Actual output of the actuator i. Equation (4.1) can be simplified as follows: Y (k) = HR(k) +D(k) (4.2) where Y(k)∈ Rm and R(k)∈ Rn are related by the above-mentioned linear transfor- mation, H is an unknown mixing-matrix and e(k) is a vector of zero-mean global error, interference or noise. Generally, vectors R(k) and Y(k) are correlated, i.e. RY R = E{Y RT } 6= 0 but the error (or colored noise D) is uncorrelated with R. Hence, our objective is to find the matrix H such that a new pair of vectors D(k) = Y (k)−HR(k) and R are no longer correlated with each other, i.e., RDR = E{DRT } = E{(Y −HR)RT } = 0. (4.3) 41 Chapter 4. Tensor Modeling and Alignment The cross-correlation matrix can be written as RDR = E{Y RT −HRRT } = RY R −HRRR. (4.4) Hence, the optimal mixing matrix for correlation cancelation can be written as Hopt = RY RR−1RR (4.5) or Hopt = E{Y RT }(E{RRT })−1. (4.6) Since we have direct and full access to the reference input, R, as well as databoxes, Y , then calculating their autocorrelations and cross correlation is possi- ble. Calculation of Equation (4.5) leads us to the mixing matrix, H. In CD control, this mixing matrix is the interaction matrix. The problem that makes current BSS techniques impractical in CD processes is the statistical independency or nongaussianity of the components. However, with regard to Cross Directional processes we are dealing with a distributed closed-loop system, where the inputs are highly correlated. 4.1.2 CD Tensor Model The paper’s characteristics before shrinkage at CD location j and time k produced by actuator r, j ∈ [1, J ], k ∈ [1,K], r ∈ [1, R], is given by xr,j,k = sr,kcr,j , (4.7) where cr,j is the jth CD response gain of actuator r before shrinkage and sr,k is the normalized temporal (Machine Direction) response of the rth actuator at time k defined by sr,k = ur,k ∗ hk, (4.8) in which ur,k is the reference input of the rth actuator at the time k, convolved with the system’s dynamic, hk. 42 Chapter 4. Tensor Modeling and Alignment The standard separability assumption of the MD and CD responses is used here. There is no assumption either regarding the CD response of the actuators or regarding the orthogonality of the columns of the interaction matrix. However, we assume the width of paper machine, J, is known or has been defined. Then the measured value in scanned databox i, i ∈ [1, I], for CD location before shrinkage j and time k, can be written as the sum over all the actuators of the product of the fading factor between actuator r and scanned databox i (ar,i) and the CD signal value excited by actuator r (xr,j,k): yi,j,k = R∑ r=1 ar,ixr,j,k, (4.9) yi,j,k = R∑ r=1 ar,icr,jsr,k. (4.10) A schematic tensor representation of the CD process is given in Figure 4.2. Equation (4.10) is known as an R-component PARAFAC decomposition of the I × J × K third-order tensor Y [39, 40]. The element yi,j,k can be seen as the element in position (i, j, k) of tensorY ∈ CI×J×K . For convenience, no noise term is shown in Equation (4.10). But of course noise is a major concern in applications, and should appear in all the equations above [113]. Now we move directly to the case where the noisy scanned data is represented as follows: yi,j,k = R∑ r=1 ar,icr,jsr,k + vi,j,k, (4.11) where vi,j,k is the measurement noise. The standard PARAFAC model for a three-way (3D) array expresses the orig- inal tensor as a sum of three-way factors, each one of which is an outer prod- uct of three vectors. Tensors consisting of the outer product of some vectors, are called rank-one. By analogy with the definition of matrix rank, the rank of a third- order tensor is defined as the minimum number of rank-one, three-way components needed to decompose it into a sum [112]. 43 Chapter 4. Tensor Modeling and Alignment 100 200 300 400 500 600 700 800 900 100 200 300 400 500 600 700 Figure 4.2: Schematic representation of the CD process concerning CD sheet shrinkage. 44 Chapter 4. Tensor Modeling and Alignment 4.2 CD Alignment by Tensor Decomposition Let us call Ar, Cr and Sr the three vectors of size I, J and K respectively repre- senting the fading coefficients ar,i, the CD response coefficients cr,j and the input coefficients to actuator number r. Equation (4.10) can be written in tensor-notation as follows: Y = R∑ r=1 Ar ◦ Cr ◦ Sr, (4.12) or, with noise, Ỹ = R∑ r=1 Ar ◦ Cr ◦ Sr + V, (4.13) Figure 4.3 illustrates the PARAFAC vector decomposition of tensor Y . Figure 4.3: Vector decomposition of the CD process corresponding to Equation 4.12. Finally the PARAFAC decomposition can be represented in matrix notation [39]. Let us define three matrices A, C and S with the respective dimensions I×R, J × R and K × R and elements ar,i = [A]r,i, cr,j = [C]r,j and sr,k = [S]r,k. We can slice the scanned signal tensor Y ∈ CI×J×K along all three dimensions. Now let us define a set of J ×K matrices Yi.., a set of K × I matrices Y.j. and a set of I×J matrices Y..k. Based on these definitions, the tensor model given in Equation (4.12) can be decomposed along any one of its three dimensions as follows: 45 Chapter 4. Tensor Modeling and Alignment Yi.. = R∑ r=1 ar,iCrS T r = CDi(A)S T , i = 1, · · · , I, Y.j. = R∑ r=1 cr,jSrA T r = SDj(C)A T , j = 1, · · · , J, Y..k = R∑ r=1 sr,kArC T r = ADk(S)C T , k = 1, · · · ,K, (4.14) where the operator Di(A) forms a diagonal matrix with the ith row of A. An illus- tration of this decomposition is shown in Figure 4.4. Figure 4.4: Matrix decomposition of the CD tensor process at each instant of time, k. 4.2.1 Alternating Least Squares Regression Although there are many published algorithms for decomposing a PARAFAC model, it is difficult to gain an overview of their relative performances [51]. Since the 46 Chapter 4. Tensor Modeling and Alignment scanner outputs are highly corrupted by noise, the principle of Alternating Least Square (ALS) is recommended for such data [85]. Whilst in the noiseless case and under considerably stronger conditions, A, C and S can indeed be found using eigen analysis [160, 163]. Moreover, ALS allows modeling of higher-order data, as well as the incorpo- ration of constraints on the parameters and further ALS is able to efficiently handle missing data and ignore the diagonal values of the data [51]. Although ALS is sig- nificantly slower compared to the newer algorithms, but it is found that the models generated by ALS are generally of better quality than any of the alternatives. The least-squares cost function associated with Equation (4.11) is minimized by means of alternating updates of one of its matrix arguments, keeping the other two matrices fixed. Since PARAFAC is a multi-linear decomposition, each update requires solving a classical linear least-squares problem, while minimizing the as- sociated cost function. f(A,C,S) = ‖Y − R∑ r=1 Ar ◦ Cr ◦ Sr‖2 (4.15) def= ∑ ijk ∣∣∣ yi,j,k −∑r ar,icr,jsr,k ∣∣∣2 (4.16) where Ar, Cr and Sr the three vectors of size I, J and K respectively representing the fading coefficients ar,i, the CD response coefficients cr,j and the input coeffi- cients to actuator number r. 4.2.2 Identifiability and Uniqueness The PARAFAC decomposition (Equation 4.14) gives us the A matrix for each instant of time k. As can be seen from its dimension this matrix represents the mapping between actuators and scanned databoxes. Linking the CD process to PARAFAC leads to a powerful identifiability result. A detailed proof of PARAFAC’s uniqueness properties is given in [40]. Uniqueness of PARAFAC decomposition was first studied Harshman et al. in [85, 86]. The first formal proof was provided by Kruskal [112]. Kruskal derived sufficient conditions for uniqueness of third- 47 Chapter 4. Tensor Modeling and Alignment order PARAFAC decompositions of real-valued tensors. Recall that for a given matrix A ∈ CI×R, the rank r(A), is equal to r iff A contains at least a collection of r linearly independent columns. The Kruskal-rank (k-rank) of A, shown with k(A), is the maximum number k such that every set of k columns of A is linearly independent [112]. The Kruskal-rank concept is the foundation of PARAFAC. Theorem 4.2.1 (Kruskal [112]) The matricesA,B andC satisfying a set of matrix- slices, which were defined in Equation (4.14), Yi.. = R∑ r=1 ar,iCrS T r = CDi(A)S T , i = 1, · · · , I are unique up to common permutation and (complex) scaling of their columns, if k(A) + k(C) + k(S) > 2(R+ 1). (4.17) Kruskal’s proof has been recently revisited and clarified by Jiang et al. in [102]. Condition 4.17 is a sufficient condition. In [102, 174], Kruskal’s condition has been shown to be not only sufficient but also necessary, if R 6 3. They provided both necessary and sufficient uniqueness conditions for PARAFAC, when one factor matrix has a full column-rank. An extension of PARAFAC’s uniqueness condition to the complex case is developed in [158, 159]. This condition is generalized to tensors of an arbitrary order in [120]. If the elements of A, C and S are drawn from a continuous distribution, the k-rank and the conventional rank coincide with a probability of 1. The sufficient condition for the PARAFAC uniqueness can be stated as [102]: min(I,R) + min(J,R) + min(K,R) > 2R+ 2. (4.18) And finally two necessary conditions for the uniqueness of general PARAFAC decomposition are: min(r(AC), r(C S), r(SA)) = R, (4.19) and min(k(A), k(C), k(S)) > 2. (4.20) 48 Chapter 4. Tensor Modeling and Alignment The PARAFAC algorithm has an advantage over matrix-based methods in that the decomposition is done under very mild conditions without requiring orthogo- nality or independence among components. Finally, PARAFAC has many statistical merits such as optimal unbiased esti- mations of the final results in ALS sense. In addition, the true and estimated load- ing matrixes A, C and S coincide when the number of factors has been correctly estimated [30]. CD Mapping and Identification In CD process control systems the actuator profile along with the normalized tem- poral response of the actuators, S, is usually known. If the shape of the CD response of the actuators is also known or if it can be effectively estimated (e.g. initial bump test), then C (or a good estimate thereof) can be assumed to be known too. In such a case we can initialize ALS with the given C and S. This significantly speeds up the convergence of finding A and improves performance because of smaller MSE. Each subsequent update reduces the value of f and improves the fit in Equation (4.15), therefore the convergence is monotonic. However, it should be mentioned that even with large model uncertainties a small number of alternating updates suf- fices to identify both the dynamic and the CD response of the actuators. In CD control of paper machines, databoxes usually outnumber actuators by factors of 5 to 10. Therefore the selection of k > 2 will obviously satisfies above- mentioned necessary and sufficient conditions for PARAFAC uniqueness given in Equations 4.18-4.20. 4.3 Results Part 4.3.1 of this section presents the result after applying our proposed method to align a simulated mis-mapped paper machine. In part 4.3.2 a field trial of our method for the alignment of the basis weight actuators of a paper machine in a Canadian paper mill is described. 49 Chapter 4. Tensor Modeling and Alignment 4.3.1 Simulation Results In the following examples at each instant of time, the designed tensor is filled with the scanning databoxes while still missing some values. In PARAFAC those missing values are easily handled using iterative estimation methods like ALS or Expectation Maximization (EM) Algorithm [21, 106, 147, 175]. The missing elements are consistently replaced with their corresponding esti- mates and the iteration is continued until no significant changes occur in the esti- mates of missing elements and the overall convergence criterion is fulfilled. It is also possible to handle missing values by a weighted regression, setting the weights of missing values to zero. This will add some bias to the decomposition. The more missing data, the higher value of the bias [175]. However, as long as the neces- sary and sufficient conditions of PARAFAC are satisfied our estimated models are unique. Example 1 In this industrial process application, the basis weight of a sheet of paper is con- trolled and the corresponding model of the actuator has been identified. The weight of the paper is controlled through an array of R = 7 identically-designed and dis- tributed slice lip actuators. The same as the model shown in Equation 2.1, the process model is given by y(z) = Gt(z)u(z) +Dt(z)d(z) (4.21) where the time step is measured in units of the sampling time (Ts = 25s), the mea- surement points of y(z) ∈ R50×1, the actuator set-points u(z) ∈ R7×1, Dt(z) ∈ R50×7 represents the transfer matrix coloring filter of the vector of white gaussian noise disturbance d(z) ∈ R50 with standard deviation of 0.1 and the dynamic of 1−0.8221z−1 1−0.999z−1 , Gt(z) ∈ R50×7 is the process transfer matrix containing both the dy- namic and the spatial responses of the system to the actuator array, with the dead time d = 3 and MD response of z −3 1−0.8221z−1 . The simulated closed-loop paper machine is controlled with a Dahlin Con- troller manipulating 7 actuators. The output of the simulated system formed 500 50 Chapter 4. Tensor Modeling and Alignment scans. Each scan contained 50 databoxes. As is true with real industrial paper machines, the effective part of scans are truncated as the unaffected edges are thrown away (arbitrarily 5 scanning points on either side of the CD process in Figure 4.5). A nominal CD process model of the actuator is generated based on the model given by Gorinevsky et al. in [71] with the following parameters: a = 10, w = 15, g = 1, δ = 2.5 100 200 300 400 500 600 700 800 100 200 300 400 500 600 Figure 4.5: Closed-loop response of actuators. The system is then forced to have a sudden sheet-wandering type of misalign- ment on actuator number 5 at scan 250. In this example it is assumed that the model of the actuator response is known. By entering the actuator’s input to the model a noise-free CD response before shrinkage is generated. A simple LS (in- stead of ALS) is applied to the above-mentioned simulated CD process to identify the location and the moment the misalignment occurred or, in another words, do the mapping. It is already proven that this mapping is unique. Figure 4.6 shows the dynamic mapping of actuator number 5. 51 Chapter 4. Tensor Modeling and Alignment 50 100 150 200 250 300 350 400 450 5 10 15 20 25 30 35 40 45 50 MD (time) Da ta Bo xe s (S pac e) PARAFAC Alignment Misalignment Paper Wandering Figure 4.6: The solid line shows the forced misalignment. The dash-dot line shows the original mapping neglecting the misalignment. The dotted line represents the databox aligned with the actuator at each moment. Example 2 The same actuator model and all the same CD process that were given in the pre- vious example are used in this simulation. This system has 9 actuators generating k = 250 scans. I = 60 databoxes (See Figure 4.7). Since the dimension J is virtually created, without loss of generality, it is arbi- trarily assumed that J = R = 9. In this case the CD response of the actuators has high-model types of uncertainties although they are perfectly aligned. Next, the PARAFAC method is performed on the system. The updated mapping is shown in Figure 4.8. As it can be seen from Figure 4.9 PARAFAC has resolved the significant issue of the model uncertainties. 52 Chapter 4. Tensor Modeling and Alignment MD (Time) C D  (S pa ce ) 50 100 150 200 250 10 20 30 40 50 60 Figure 4.7: The scanning databoxes of a system with no misalignment but with high model uncertainties. Alignment Ac tua tor  N um be r 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 −0.2 0 0.2 0.4 0.6 0.8 Figure 4.8: The actuator mapping matches with the expected uniform distribution. 53 Chapter 4. Tensor Modeling and Alignment 2 4 6 8 10 12 14 16 18 20 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 CD (Databox Number) True Model PARAFAC Initial Model Figure 4.9: Identification of the CD response of an actuator using PARAFAC and ALS. The initial model has high level of model uncertainties while has no mis- alignment. Example 3 In this example again the same CD process given in the previous examples is used for a system with R = 9 actuators for K = 250 scans. Here, the system was forced to have misalignment (unknown 3 databox shrinkage) for the actuators 7, 8 and 9. At the same time the actuator response is modeled with a model based upon one that has high uncertainties. As can be seen from Figure 4.10 this amount of shrinkage and these model uncertainties make the system unstable around databox numbers 30 - 40. Each scan contains I = 57 databoxes and let us arbitrarily choose J = 9. The PARAFAC method using ALS is performed on the system where there is problem- atic CD mapping combined with high model uncertainties. The following pictures show the identification and mapping results: The identified CD response of an actuator is shown in Figure 4.12. Since the system was designed to have the forced misalignment of one actuator spacing for the last 3 actuators, not surprisingly we find that the results are not as well resolved 54 Chapter 4. Tensor Modeling and Alignment as a system that was perfectly aligned. This perturbation we see around databoxes 30− 50 is negligible and, practically speaking is neglected for distributed systems. Its appearance is caused by large magnitude in misalignment and the resulting instability to the system. However, as it was discussed and shown in Example 2 after realignment a much better model will be obtained. C D  (S pa ce ) MD (Time) 50 100 150 200 250 5 10 15 20 25 30 35 40 45 50 55 Figure 4.10: An image of scanning databoxes with local instability resulted from both model uncertainties and misalignment. 4.3.2 Experimental Results We recently tested our new alignment method on a modern paper machine. The system consists of 78 dilution actuators for weight control, 40 actuators for mois- ture control and 600 databoxes. Our proposed method was used to align the slice lip actuators. In this industrial experiment, an acceptable model of the actuator response was available. Therefore, instead of applying ALS to find mapping and modeling, just LS was used. Seventy-eight actuators were spread across a 10 me- ter wide papermaking machine; each actuator has response width approximately 20cm. This paper machine runs at speeds up to approximately 600 m/min and 55 Chapter 4. Tensor Modeling and Alignment Alignment Ac tu ato r N um be r 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 −0.2 0 0.2 0.4 0.6 0.8 1 Figure 4.11: Theapping obtained from PARAFAC shows obvious shift on the re- sponse of the last three actuators resulted from shrinkage on that area. 5 10 15 20 25 30 35 40 45 50 55 −0.5 0 0.5 1 Space (CD) Ex ci ta tio n Va lu e   Ideal Model PARAFAC Initial Model Figure 4.12: Identification of the CD response of an actuator using PARAFAC. The initial model has high level of model uncertainties while the system suffers from misalignment. 56 Chapter 4. Tensor Modeling and Alignment each scan takes 25s, the delay Td ≈ 3 scans. The scanner records paper properties at I=560 points (i.e. the datapoints are about 1.5cm apart). As can be seen from Figure 4.13, there is a significant local instability around databoxes 120 to 180, which may be caused by an actuator misalignment. This problem, which is common in industrial CD control, is called “picket fence” or “picketing” [91]. Our proposed alignment method has been tested on this paper mill to find the mapping function. Subsequently, the presence of a misalignment in that area, as well as our method’s performance were verified by a manual bump test (see Figure 4.14). The exact time and cause of the misalignment were checked and clarified with the mill’s staff. The staff reported that the detected misalignment happened around the 10th scan after the operators increased the system’s speed (See Figure 4.15). MD (time) CD  (D ata bo x n um be r) 50 100 150 200 250 300 50 100 150 200 250 300 350 400 450 500 550 Figure 4.13: Local instability (picketing) around databoxes 120 to 180. Figures 4.16, 4.17 and 4.18 illustrates the realtime mapping of 9 actuators (ac- tuators: 7, 16, 24, 32, 40, 48, 56, 63 and 70) at various positions across the sheet of paper. The black oval drawn on the image 4.14 shows the bump test of actua- tors. The mapping obtained from new bump test matches with the databox number, 57 Chapter 4. Tensor Modeling and Alignment CD  (D ata bo x n um be r) MD (time) 50 100 150 200 250 300 50 100 150 200 250 300 350 400 450 500 550 Figure 4.14: New alignment of the actuators based on the bump test. Black ellip- soid shows the location of the manual bump tests. 20 40 60 80 100 120 140 160 180 140 145 150 155 160 165 170 175 180 185 190 MD (Time) CD  (D ata box es) Figure 4.15: Detected sudden alignment change of the actuator no. 24. 58 Chapter 4. Tensor Modeling and Alignment which was driven from our proposed technique. These results were authenticated with the results of two different bump tests. For these bump tests, the above-mentioned nine actuators were simultaneously bumped, each with an excitation sequence of 25% up(down) followed by 25% down(up) of their average excitation. The result obtained from this new bump test and from the initial test, plus the results from our proposed method, labeled Tensor, are collectively compared in the Table 4.1. Table 4.1: Mapping Change Observed with PARAFAC Decomposition and Justi- fied with the Bump Test.`````````````̀Databox # Actuator # 7 16 24 32 40 48 56 63 70 Initial Bump 36 96 159 231 295 359 423 479 525 New Bump 37 107 169 230 293 357 419 475 524 Tensor 36 106 167 230 290 353 416 476 526 In fact the first and last rows of this table represent some elements of the rows and columns of matrix A ∈ R78×560, discussed in Equation 4.14. There are eight scanning points corresponding to each actuator. This means that only significant alignment changes, larger than eight scanning points, will appear as a misalignment instability and any small mapping changes that are less than eight can be neglected. This fact plus the closeness of mappings obtained from Bump test and our new proposed method (see third and fourth rows of Table 4.1) provide strong indication of the validity and accuracy of our method. Figure 4.14 shows the improved outputs after applying the new alignment based on the bump test. Figures 4.16 and 4.17 also show that actuator mapping is continuously subject to change. Unfortunately, there is no way to tell if those changes are caused by actual misalignment or if they are caused by noise. Figure 4.19 shows the manual open-loop bump test of actuator number 32. In that figure the more yellowish and bluish parts corresponds to the center point of the actuator. Those two peaks were expected to be observed at just one databox. However, as Figure 4.19 shows, there was a misalignment or mapping 59 Chapter 4. Tensor Modeling and Alignment 20 40 60 80 100 120 140 160 180 5 10 15 20 25 30 35 40 45 50 55 PARAFAC Time Da ta Bo x N um be r (a) Aligned databox to the actuator no. 07 vs time. 20 40 60 80 100 120 140 160 180 85 90 95 100 105 110 115 120 125 PARAFAC Time Da ta Bo x N um be r (b) Aligned databox to the actuator no. 16 vs time. 20 40 60 80 100 120 140 160 180 140 145 150 155 160 165 170 175 180 185 190 PARAFAC Time Da ta Bo x N um be r (c) Aligned databox to the actuator no. 24 vs time. Figure 4.16: PARAFAC results. 60 Chapter 4. Tensor Modeling and Alignment 20 40 60 80 100 120 140 160 180 205 210 215 220 225 230 235 240 245 250 PARAFAC Time Da ta Bo x N um be r (a) Aligned databox to the actuator no. 32 vs time. 20 40 60 80 100 120 140 160 180 275 280 285 290 295 300 305 310 PARAFAC Time Da ta Bo x N um be r (b) Aligned databox to the actuator no. 40 vs time. 20 40 60 80 100 120 140 160 180 335 340 345 350 355 360 365 370 375 PARAFAC Time Da ta Bo x N um be r (c) Aligned databox to the actuator no. 48 vs time. Figure 4.17: PARAFAC results. 61 Chapter 4. Tensor Modeling and Alignment 20 40 60 80 100 120 140 160 180 400 405 410 415 420 425 430 435 440 PARAFAC Time Da ta Bo x N um be r (a) Aligned databox to the actuator no. 56 vs time. 20 40 60 80 100 120 140 160 180 455 460 465 470 475 480 485 490 495 PARAFAC Time Da ta Bo x N um be r (b) Aligned databox to the actuator no. 63 vs time. 20 40 60 80 100 120 140 160 180 510 515 520 525 530 535 540 545 550 PARAFAC Time Da ta Bo x N um be r (c) Aligned databox to the actuator no. 70 vs time. Figure 4.18: PARAFAC results. 62 Chapter 4. Tensor Modeling and Alignment change even while we were performing a bump test. It is interesting that our pro- posed method can even detect this small, sudden alignment change on actuator 32, around scan number 75 (see Figure 4.17(a)). The mapping of the actuators, A, obtained from Figures 4.16 and 4.17 at an instant of time is shown in Figure 4.20. This figure also shows that shrinkage is not symmetric across the sheet. This fact conflicts with the assumption of many previously mentioned articles. Figure 4.19: Alignment changes while performing bump test. The overshoot and undershoot a bump test responses of the actuator is not observed at the same databox. 4.4 Conclusions A new tensor-based model of a paper-making machine and an online non-intrusive method of aligning CD actuators and their corresponding response centers have been proposed in this chapter. The method is based on new tensor modeling of CD processes. The tensor is constructed according to the dimensions of the CD actuators, the output of actuators before shrinkage, the observed scanned databoxes 63 Chapter 4. Tensor Modeling and Alignment 0 10 20 30 40 50 60 70 0 100 200 300 400 500 R (Actuator number) I ( Da ta bo x nu m be r)   PARAFAC Uniform Distribution Initial Bump Updated Bump Figure 4.20: Identified matrix A. and time. Subsequently, this model is decomposed by the PARAFAC decomposition method giving a unique alignment identification at each instant of time. Unlike current industrial alignment methods, which are based on bump-test data, our proposed method can be applied without any intrusion into the normal operation of the system. Further, we obtained an explicit expression to account for the shrinkage [57, 58]. 64 Chapter 5 Multivariable Windsurfer Adaptive Robust Control Cross-directional processes are inherently ill-conditioned, as documented by Kristinn- son [111] and Featherstone et al. [60]. Even in steady-state conditions, many sin- gular values are vanishingly small [59, 78, 169]. This problem, compounded with model uncertainties, makes both identification and controller design a challenge for theoreticians and practitioners alike. Stewart et al. [168] describe an industrial tun- ing tool for cross-directional controllers, based on a two dimensional loop-shaping algorithm, using models identified from a bump-test experiment. This Chapter is divided into the three following main sections: • Application of the classical Windsurfing iterative adaptive robust control as a new approach for the control of CD processes. • MultivariableH∞-based Windsurfing. • New Enhancements on Windsurfing. In the first part of this Chapter we present a novel Adaptive Robust Control approach to the multivariable CD process of continuous web manufacturing. We assume that large-scale CD processes can be approximated as Square Circulant Matrices (SCMs) and then diagonalized by real Fourier matrices, allowing analysis in terms of a family of Single-Input Single-Output (SISO) transfer functions across the spatial frequencies. Then, we apply discretized Windsurfing adaptive robust control to each individual separated spatial frequency, starting with a stable initial model and a robust stabilizing controller at each spatial frequency. 65 Chapter 5. Multivariable Windsurfer Adaptive Robust Control This approach allows the dynamic bandwidth of the closed-loop system to be increased progressively at each spatial frequency through an iterative control- relevant system identification and control-design procedure. The method deals with model uncertainties, disturbances and measurement noise issues. Finally we provide a simulated result in order to illustrate the performance of the applied method. Next, the new modifiedH∞-based Windsurfing is applied to the CD processes. The combination of SFD and theH∞ allows the dynamic bandwidth of the closed- loop system to be increased cautiously and progressively at each spatial frequency through an iterative control-relevant system-identification and control-design pro- cedure. Windsurfing also provides robust stability through iterations. The method deals with both Machine-Direction (MD) and CD model uncertainties as well as measurement noise issues. 5.1 Introduction CD processes are large scale systems with hundreds of actuators, requiring thou- sands of measurements and hundreds of hard input constraints. The computational complexity rises quickly with the number of actuators [60]. Since the dimensions of CD processes is, for the most part, larger than multivariable control systems, the synthesis of CD systems can be intractable with available synthesis techniques − even using modern computing power [60] (e.g. H2 and H∞ [162, 195]). In addi- tion, the practical implementation of a controller must be considered in the design evaluation. In a centralized controller each actuator is connected to every sensor. Although this way satisfy performance requirements, connectivity will limit its im- plementation. Featherstone et al. [60] propose a technique for designing an input that will minimize a confidence-ellipsoid associated with the identified response. The re- sulting input does result in all actuators moving and is rich in high frequency con- tent. However, the method is very computationally intensive and requires several open-loop experiments. Once a model with uncertainty bounds (estimated or ar- bitrarily assigned) is obtained, robust control theory can be used to design a CD controller. 66 Chapter 5. Multivariable Windsurfer Adaptive Robust Control In this domain, an industrial practice is represented by the two-dimensional loop-shaping procedure of Stewart et al. [164, 168, 169]. Where applying the spatial frequency decomposition allows the computational complexity of the con- troller synthesis to be reduced [8, 38, 95, 193]. The loop-shaping control synthesis [124, 194] is extended to spatially-invariant, two-dimensional frequency domain systems [47, 114, 162, 193]. More recently, emphasis has been on the development of multi-actuator ar- rays and multi-property control systems based on Model Predictive Control (MPC), while dealing with high model uncertainty [7, 156]. Multivariable CD control is used by Berggren et al. [9, 10] to achieve optimal profile performance of coupled processes. An identification technique is used by Saffer et al. [152] to identify dual-headbox of a paper machine, where the identified model is used within a lin- ear programing-based MPC. Haznedar et al. [87] have also introduced a solution for the multiple property control problem. There are quite number of fundamental problems with adaptive and robust con- trol pertinent to our approach to adaptive robust control. Mostly, in traditional adaptive control the unknown plant is represented by a model in which everything is known except for the values of a finite number of parameters. Then a principle of certainty-equivalence is invoked to update the controller, which is not cautious and therefore is in danger of causing instability [180]. The fundamental problems with adaptive control can be summarized into three classes [2]. • The problem of changing experimental conditions, given a close but inexact model: The model might have good response with an initial controller, but it is not necessarily guaranteed to be good with a modified controller. • Impractical control objectives or other dificulties, e.g., very large transient signal values or difficulty in learning the controller or doubt in algorithm’s convergence. • Transient instability problems. In contrast, a robust controller is designed on the basis of a nominal model of the plant with associated parametric and unstructured model uncertainties explic- 67 Chapter 5. Multivariable Windsurfer Adaptive Robust Control itly taken into account. Only the a priori information on the model is considered and those characteristics of the plant are neglected that could be learned while it is being controlled. Consequently, in robust control, performance is sacrificed for a conservative design [2]. The main point of this chapter is to apply a multivariable adaptive robust con- trol approach with the aim to decrease the uncertainties of a given CD model and achieve a better controller design. We achieve this goal by taking advantage of the properties of the Windsurfer’s approach to adaptive control [115]. The key element in the Windsurfing approach to adaptive robust control is to increase the reference system’s bandwidth and decrease the model’s uncertainties. Section 5.1.1 presents an overview of existing techniques for two-dimensional frequency domain decomposition in order to improve the design, both computa- tionally and conceptually. Section 5.2 describes the constructive design technique that is central to this work. The modeling, design and implementation are the sub- ject of this section. Finally in Section 5.2.5, the designed multivariable iterative adaptive robust control technique is applied to a simulated industrial papermaking cross-directional control problem. 5.1.1 Cross-Directional Process Model With the measured paper profile y, the actuator settings u, and the disturbance input d as time-varying vectors in Cn, where n stands for the number of actuators, the process is modelled as: y(z) = G(z)u(z) + d(z). (5.1) We assume that the MD and CD contributions to the process dynamics are separable (this is common—see [14, 23]), with a matrix transfer function G of the form G(z) = h(z)G0, h(z) = z−Td 1− az−1 . (5.2) Here the scalar factor h captures the simple, first-order machine-direction dy- namics, with delay Td, and G0 is an n× n matrix describing the cross-directional 68 Chapter 5. Multivariable Windsurfer Adaptive Robust Control response. Each component of the actuator vector u corresponds to a column ofG0. The actuators are designed to be identical, so G0 is a band-diagonal symmetric Toeplitz matrix: G0 =  g1 g2 · · · gq 0 · · · · · · · · · · · · 0 g2 g1 g2 · · · gq 0 . . . . . . . . . ... ... g2 g1 g2 · · · gq . . . . . . . . . ... gq . . . g2 g1 g2 . . . . . . . . . . . . ... 0 gq . . . g2 . . . . . . . . . gq 0 ... ... 0 gq . . . . . . . . . g2 . . . gq 0 ... . . . . . . . . . . . . g2 g1 g2 . . . gq ... . . . . . . . . . gq . . . g2 g1 g2 ... ... . . . . . . . . . 0 gq · · · g2 g1 g2 0 · · · · · · · · · · · · 0 gq · · · g2 g1  ∈ Rn×n (5.3) We adopt Matlab’s notation, writing G0 = toeplitz{g1, . . . , gq, 0, . . . , 0} (5.4) The influence of each actuator is confined to a small spatial region, so the diagonal band in G0 is relatively narrow: q  n. 5.2 Classical Windsurfer Adaptive Robust Control We adopt the idea of an iterative adaptive robust control to resolve many issues with current adaptive control and robust control. Adaptive robust control recognizes that a plant characteristics can differ greatly from the estimated model at any time, particularly during the initial learning stage. It uses a posteriori knowledge about the plant to reduce conservatism in a robust control design [2]. Iterative adaptive robust control typically involves identification followed by control parameter adjustment. Several approaches are possible for iterative tuning. A generic iterative adaptive control system is illustrated in Figure 5.1. 69 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.1: Iterative adaptive control system. Here we use the intuitively appealing “Windsurfer approach” to iterative adap- tive robust control [5, 110, 115], which is akin to an indirect iterative adaptive control scheme in which a closed-loop bandwidth is gradually increased at every iteration [180]. Lee et al. [115] propose an iterative algorithm, which gradually and safely broadens the bandwidth of a closed-loop system. Iterative repetition of the controller design and system identification decreases model uncertainties at the excited frequencies. 5.2.1 Spatial Frequency Decomposition Because of the non-analytic nature of the singular values of a dynamical system, multivariable control problems cannot be controlled in a one-loop-at-a-time design [47]. However, Hovd et al. [94] have shown that there is an exemption for sepa- rable systems [161] and circulant systems [114, 193]. There are quite number of techniques to solve the problem [52, 114, 161, 164, 193] but the details of their discussion is beyond this context. Spatial Frequency Decomposition (SFD) meth- ods are based on the same concept, so the SFD of square matrices is chosen to be 70 Chapter 5. Multivariable Windsurfer Adaptive Robust Control explained here. A small change to G0 makes a big difference in the complexity of our problem [14, 23, 91, 114]. Consider the matrix Ĝ0 =  g1 g2 · · · gq 0 · · · 0 gq · · · g2 g2 g1 g2 · · · gq 0 . . . . . . . . . ... ... g2 g1 g2 · · · gq . . . . . . . . . gq gq . . . g2 g1 g2 . . . . . . . . . . . . 0 0 gq . . . g2 . . . . . . . . . gq 0 ... ... 0 gq . . . . . . . . . g2 . . . gq 0 ... . . . . . . . . . . . . g2 g1 g2 . . . gq gq . . . . . . . . . gq . . . g2 g1 g2 ... ... . . . . . . . . . 0 gq · · · g2 g1 g2 g2 · · · gq 0 · · · 0 gq · · · g2 g1  ∈ Rn×n (5.5) This is the “circulant extension of G0” [38], abbreviated as Ĝ0 = toeplitz{g1, . . . , gq, 0, . . . , 0, gq, . . . , g2}. (5.6) The matrices G0 and Ĝ0 differ only in the ‘ears’ of Ĝ0, as illustrated in Figure 5.2. We write δG0 = Ĝ0 −G0 = toeplitz{0, . . . , 0, gq, . . . , g2}. (5.7) Since each actuator’s response (i.e. non-zero bands in G0) on the transfer ma- trices of CD controlled paper machines is relatively narrow (except for some appli- cations to heavy paper grades), the resulting perturbation is generally small. (For details, see [164]). The first design requirement of any feedback control system is the satisfaction of its internal closed-loop stability. The small gain theorem sup- plies a sufficient condition for the internal stability of the true system [164, 195]. It is well known [149], but still remarkable, that the basis for Rn associated with the Discrete Fourier Transform simultaneously diagonalizes all symmetric circulant matrices in Rn×n. That is, there is a single n × n Fourier Matrix “F ” 71 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.2: Nonzero elements of a banded symmetric circulant spatial interaction matrix (left), the associated band-diagonal Toeplitz symmetric matrix (center), and the difference ‘ears’ (right). with real entries such that F−1 = F T and FCF T is diagonal for every symmetric circulant C. The columns of F are the “Fourier basis” [8, 95]. Thus, the MIMO system obtained from Equation 5.1 by replacing G0 with Ĝ0 can be decoupled as: Y (ν, z) = G̃(ν, z)U(ν, z) +D(ν, z). (5.8) where ν stands for the decomposed spatial frequency. Therefore, the change of variables outlined above reduces our task in solving the n independent SISO prob- lems: Ym(z) = h(z)g̃(νm)Um(z) +Dm(z), m = 1, . . . , n. (5.9) Theorem 5.2.1 (Hovd et al. [95]) If the generalized plant Ĝ(z) is composed of symmetric circulant transfer matrices, and if there exists a controller such that the closed-loop is any of (i) H2 optimal, or (ii) H∞ optimal, or (iii) H∞ admissible, then there exists a symmetric circulant controller K̂(z) that satisfies the respective criterion. 72 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Theorem 5.2.2 (Bamieh et al. [8]) [Controller Synthesis] Feedback controllers satisfying any one of the characteristics in Theorem 5.2.1 may be synthesized as a family of SISO controllers K̃(νm, z), independently for each spatial frequency νm ∈ {ν1, · · · , νn}, then constructing the multivariable controller K̂(z) = F Tdiag{K̃(ν1, z), · · · , K̃(νn, z)}F. (5.10) Theorem 5.2.3 (Stewart et al. [168]) [Controller implementation] If a family of causal finite-order rational SISO transfer functions K̃(νm, z) for m ∈ {1, · · · , n} has been synthesized (e.g. satisfying Theorem 5.2.1), then the full MIMO con- trollers K̂(z) in Equation (5.10) can be implemented in the discrete-time domain. Once n independent controllers K̃(νm, z), m = 1, . . . , n, have been se- lected, we recover a multivariable controller suitable for use in the original co- ordinates by recalling K̂(z) = F T K̃(z)F. (5.11) From this point system identification and controller design are performed on the decomposed circulant extension of the ‘true’ system. Details of the Spatial Frequency Decompositions (SFD) just outlined can be found in [14, 23, 91, 95, 114]. 5.2.2 Controller Design The Iterative Internal Model Control (IMC) method is applied in the control design step of the windsurfer approach where the reference input is a step function. See Chapters 4 and 7 of Reference [130] and References [115–118] for more defini- tions and theorems. The robust controller design will be based on a detuned H2 optimal controller given an estimated model and desired closed-loop bandwidth λb. This setup allows us to combine Windsurfing Adaptive Robust Control with Spatial Frequency Decomposition (See Theorem 5.2.3). For CD control design, control would first take place for low spatial and dy- namical frequencies where the model uncertainty is very low. This control could easily be achieved, e.g., by feeding an initial controller and increasing resolution at subsequent iterations. 73 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.3: Internal model control (IMC) structure. The iterative CD controller design process aims at increasing the closed-loop bandwidth by updating the existing controller. Some validation tests should be carried out to find whether the model uncertainties at the excited frequencies will significantly degrade the performance. The goal is to find the best possible CD model and closed loop response to make the paper sheet as uniform as possible. This tuning can be summarized as increasing the bandwidth of excitation at each and every spatial frequency (see Figure 5.4). There is a tradeoff between robust performance and robust stability. Conse- quently a region Ωh for the satisfaction of the initial robust stability at high spatial and dynamical frequencies has been designed as follows [164]: Ωh = {(ν, ω) : |G̃(ν, eiω)| < 0.1} (5.12) where ν stands for the spatial frequency and ω stands for the dynamical frequency. This region can be used for the desired bandwidth at each spatial frequency. Our goal is to increase the closed-loop bandwidth at each spatial frequency and find a better CD model and controller. This is illustrated by the horizontal arrow in Figure 74 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.4: Bandwidth increasing of the controller design. The entire shaded zone is the region Ωh. Its “boundary” is the curve where |G̃(ν, eiω)| = 0.1. 5.4. To ensure stability and good transient behavior during the tuning period of adaptive controllers, the controller parameters are adjusted at fixed time intervals instead of continuously. The time axis is segmented into intervals so that during each interval a fixed-gain controller is used, based on the model obtained at the end of the previous interval. During the first interval a known stabilizing controller is used. One is available in most practical applications. 5.2.3 System Identification The a priori assumptions for Windsurfing are as follows [2]: Assumption 1 The plant is linear and time invariant and has a stable transfer function. Assumption 2 The plant is stabilized by a known initial controller. (This con- troller almost certainly has low bandwidth and low authority.) 75 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Assumption 3 A reasonable (in particular, stable) transfer function model of the plant is available. These assumptions are intentionally modest. There is no assumption on the accuracy of the initial model, the order of the model, the relative degree of the plant and whether it is minimum phase or not [2]. Theorem 5.2.4 Let Kj,i = (1 − GiQj,i)−1Qj,i stabilize G and Gi, where Qj,i is a proper stable transfer function. Then G can be parametrized by a strictly proper stable transfer function Rj,i via G = Gi + Rj,i 1−Rj,iQj,i , (5.13) and H can also be written as H = Sj,i Di −Rj,iXj,i , (5.14) Proof See Lee et al. [115]. Theorem 5.2.5 (Lee et al. [115]) Let the controller Kj,i stabilize the plant G and the model Gi = NiDi , where Ni and Di are stable proper transfer functions, and let Kj,i = Xj,i Yj,i , where Xj,i and Yj,i are stable proper transfer functions satisfying the Bezout identity NiXj,i + DiYj,i = 1. Let Gi+1 be another model of G, also stabilized by Kj,i. Then it has a description Gi+1 = Ni + R̂j,iYj,i Di − R̂j,iXj,i , (5.15) where R̂j,i is a stable proper transfer function. Proof See Lee et al. [115]. Corollary 1 Let controllerKj,i = (1−GiQj,i)−1Qj,i satisfy Theorem 5.2.4. Then the updated model Gi+1 can also be stabilized with the same controller and be parametrized by a strictly proper stable transfer function R̂j,i. 76 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Proof It is shown in Theorem 5.2.5 that the updated model Gi+1 can be given by Gi+1 = Ni + R̂j,iYj,i Di − R̂j,iXj,i , then by replacing Xj,i with Qj,i and Yj,i with 1−Qj,iGi we obtain Gi+1 = Ni + R̂j,i(1−Qj,iGi) Di − R̂j,iQj,i . Without loss of generality let us assume Gi = Ni and Di = 1. So we have Gi+1 = Gi + R̂j,i(1−Qj,iGi) 1− R̂j,iQj,i = Gi − R̂j,iQj,iGi + R̂j,i 1− R̂j,iQj,i = Gi(1− R̂j,iQj,i) + R̂j,i 1− R̂j,iQj,i finally we find Gi+1 = Gi + R̂j,i 1− R̂j,iQj,i . (5.16) Estimation of the closed loop dynamics will be based on an associated open- loop estimation problem using dual-Youla parametrization. By using Hansen’s system identification procedure embedded in the windsurfer approach, we can transform the problem of identifying G in closed-loop into an open-loop system identification problem [116]. For this purpose we apply an algorithm such as Box-Jenkins to obtain an esti- mate of Rj,i and Sj,i, which satisfies Hansen’s open-loop representation β = Rj,iα+ Sj,ie. (5.17) 77 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Since the “input” α is independent of the measurement noise, which affects the “output” β, identification of the stable transfer function Rj,i and Sj,i (or equiva- lently Ψj,i in Figure 5.5) is an open-loop identification problem. Signals α and β are filtered by L to shape the bias-distribution of the estimate resulted from under- modeling to make the model error small in the appropriate frequency range. The appropriate signal model for this system identification is shown in Figure 5.5. Figure 5.5: Identification of Rj,i and Ψj,i [118]. A probing signal could be designed using dual adaptive control principles, e.g., by adding a filtered model uncertainty in the performance index to be minimized by the controller. This would contribute to reduce the uncertainty at frequencies higher than those currently in the bandwidth of the controller. Once the uncertainty is reduced to a satisfactory level, then a new controller with higher bandwidth can be designed and implemented for the next iteration. Figure 5.6 shows how this scheme works. The jth time period is denoted by Tj ; the transfer function estimate of the plant at the end of interval Tj is denoted by Ĝi and the controller used during time-period Tj by Ĉi,j . 78 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.6: Block diagram of the CD control applying multivariable Windsurfer adaptive robust control. 79 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 5.2.4 Multivariable Windsurfing Algorithm The multivariable system identification and control design procedure can be sum- marized in the following algorithm [55, 56]. Step 1 Use the circulant approximation to the interaction matrix to do the spatial frequency decomposition. Step 2 Set i = 1, j = 1 and ν = 1. Step 3 Factorize Gi(ν, z) as Gi(ν, z) = [Gi]m[Gi]a (5.18) where Gi(ν, z) is the νth decomposed spatial frequency model of the plant, [Gi]m is the minimum phase factor of Gi and [Gi]a is the associated allpass factor of Gi. [Gi]a = z−N h∏ j=1 (1− (αj)−1)(z − (αj)) (1− (αj))(z − (αj)−1) , αj > 1 (5.19) N is chosen so that [Gi]m is semi-proper. Step 4 Design Kj,i = Qj,i 1−Qj,iGi , (5.20) with Qj,i = [Gi]−1m Fj,i, where the augmented low pass filter Fj,i is chosen to ensure the robustness of the controller [130]. The simplest form for Fj,i(z) in the IMC-based control design is a first-order, one-parameter low-pass filter: Fj,i(z) = (1− λj,i)z z − λj,i . (5.21) Here λj,i is chosen so that Kj,i robustly stabilizes Gi. If the stabilizing controller meets the performance requirements of the closed-loop system 80 Chapter 5. Multivariable Windsurfer Adaptive Robust Control then go to the Step 5, else choose the next possible spatial frequency and jump to Step 3. Since IMC always has integral action, it is insensitive to model uncertainties at low frequencies. Moreover, the controller designed using IMC induces a natural factorization in the parametrization of the unknown transfer function of the plant. This enables the system identification problem to be solved efficiently [115]. Step 5 Let j=j+1 and set λj,i = λj−1,i −  for small  > 0 and redesign the controller Kj,i using the equations given in Step 4. Repeat this step if Kj,i satisfies the robust performance requirement. If the design produces a ro- bust stabilizing controller with the closed-loop system satisfying the speci- fied bandwidth, then choose the next possible spatial frequency; else proceed to the next step. As illustrated in Figure 5.7, Closed-loop output error ξj,i can be defined as ξj,i = yj,i − T̄j,ir. Step 5 is repeated and the bandwidth is increased cautiously while the ex- isting closed-loop system has robust performance. A stage can be reached before the occurrence of instability; where the modeling errors associated with Gi make a significant contribution to Ji: Ji = ‖νi‖22 > σ, (5.22) where νi = T̃ir is the tracking error. However, the tracking error νi cannot be measured directly. It can only be estimated from the closed-loop output error ξj,i = νj,i + wj,i, where wj,i is the effect of the noise disturbance e on the actual closed-loop output (See Figure 5.7). A closed-loop system might lose its robust performance if the bandwidth of the controller becomes too large. However, Morari et al. [130] have proved that if the existing closed-loop system has robust stability and the designed closed-loop bandwidth is increased cautiously, the robust stability is guaranteed for the new bandwidth. In other words, the constant  used in Step 5 should be sufficiently small. 81 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.7: Iterative closed-loop output-error for model number i at the jth con- troller update. Increasing the bandwidth depends on the satisfaction of the robust perfor- mance at the end of each iteration. Hansen’s indirect parametric closed-loop system identification [83, 84] has been adopted. Step 6 Perform control-relevant indirect parametric closed-loop system identifica- tion to update the model. Step 7 If Gi+1 is stable, then find the reduced-order model Gi+1 = arg min η ∥∥∥ GKj,i1+GKj,i − ηKj,i1+ηKj,i ∥∥∥∞ (5.23) otherwise choose the next possible spatial frequency (ν = ν + 1) and jump back to Step 3. Iterative algorithms can result in a dimension’s explosion; e.g., using a model Gi of degree n and R of degree 2 (an underestimate) results in a new model of degree 2n+ 2. Among the many ways to carry out model reduction, a frequency-weighted balanced-truncation scheme is used to reduce the order Gi+1 [6]. 82 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Step 8 Set Gi = Gi+1 and then return to Step 3, if Gj,i satisfies model valida- tion criteria. To reject an old model and/or accept a new model some model validation method should be applied. There are two different methods to val- idate a model: frequency domain (Power Spectrum) [117] and time domain (Correlation Based method) [121]. Correlation function estimates are too sensitive, in the sense that they tend to invalidate a model before it would be necessary and make it possible to identify a better model. Moreover, power spectrum estimates not only suggest that when a model becomes inadequate, they also indicate the frequency range where the SNR is high enough for identification [118]. Therefore, we used Power Spectrum method to validate the model. Remark 1 At each iteration all the decomposed identified models and the con- trollers should be reconstructed back to one multivariable controller and model. According to current literature on iterative control, the first few iterations are often enough to produce acceptable performance [65]. For constraint handling, the method described in Rojas et al. [146] could also prove useful, as it considers the constraints in order of increasing frequency. However, as noted by several authors, for a single-array situation, constraints are unlikely to become active for a scheme that is only attempting to control in directions that can be robustly controlled. Remark 2 Without loss of generality the same model and technique can be used for multivariable systems with lower or higher number of actuators. Intuitively speaking, the higher the number of the actuators, the narrower the additive band. The narrower the band, the smaller the size of perturbation introduced by the cir- culant extension. 5.2.5 Results In this section, we will demonstrate the effectiveness of our proposed method through two examples involving the control of a simulated CD multivariable sys- tem. The design steps outlined above are presented in terms of a simulated numer- ical example from the papermaking industry. 83 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10−3 10−2 10−1 100 Frequency (cycle/m) SV D   SVD of Toeplitz Matrix SVD of Extended Circulant Figure 5.8: Comparison between the Singular Value Decomposition (SVD) of the toeplitz matrix and its circulant extension. Example 1 In this simulated process, the basis weight of a sheet of paper is controlled and the corresponding model of the actuator has been identified. The weight of the paper is controlled through an array of n=100 identically designed and distributed slice lip actuators. The process model is given by y(z) = G(z)u(z) +D(z)d(z) (5.24) where the time step is measured in units of the sampling time (Ts = 25s), D(z) ∈ C100×100 represents the transfer matrix coloring filter of the vector of white gaus- sian noise disturbance d(z) ∈ C100 with standard deviation of 0.05 and G(z) ∈ C100×100 is the process transfer matrix containing both the dynamic and the spatial responses of the system to the actuator array. 84 Chapter 5. Multivariable Windsurfer Adaptive Robust Control The two transfer matrices in (5.24) are given by factors, G(z) = (I −Az−1)−1(B · z−d) (5.25) D(z) = (I −Hz−1)−1(E · z−d) (5.26) where the dead time d = 3. The matricesA,B,H,E ∈ C100×100 are all symmetric band-diagonal constant matrices. A = toeplitz{0.8221, 0, · · · , 0}, B = toeplitz{ḡ, 0, · · · , 0}, H = toeplitz{0.999, 0, · · · , 0}, E = toeplitz{0.8221, 0, · · · , 0}, (5.27) with ḡ = [ −0.0814,−0.0455,−0.0047, 0.0017, 0.0003 ] (5.28) The nominal CD process model of the actuator is shown in Figure 5.9. The same initial model G(z) given in Equation (5.29) has been used for all spatial frequencies. This model mismodels the DC gain and fails to capture the response: G(z) = 0.1z−3 1− 0.9z−1 (5.29) The initial controller has the format described by Stewart [164] as following: K(z) = [I + S(z)]−1C(z) (5.30) where C, S and K ∈ R100×100. The circulant extension of the above-mentioned matrices is required to design the new controller. Figure 5.10 shows the result of iterative system identification on one spatial 85 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.9: Nominal model of the plant. frequency. Iteratively, the same procedure has been done for all other spatial fre- quencies. Figure 5.11 shows the initial and updated step response of the discussed system to the reference input with spatial square shape. By comparing graphs (b) and (c) of Figure 5.11, we observe that not only has the MD response improved but also the resolution of the CD response has changed from an almost sinusoid shape to a very high resolution square response similar to the reference input shown in Figure 5.11(a). For reasonable control one needs to have a system response in steady state. On our test paper machine the slowest dynamical time constant is 0.8221 with time delay of 3 scans. In order to reach a settling time within a ±2% error, a total of 6 scans are required. Therefore, to be on the safe side 20 scans was chosen as the length for each iteration. To obtain the given results shown in Figure 5.11, ideally 6-8 iterations are needed at lowest spatial frequencies. However, after just a few iterations at high spatial frequencies (that included at least one iteration for the controller design and one iteration to test the new model), the system will be stopped. Depending on the signal to noise ratio, model uncertainties (both at high spatial 86 Chapter 5. Multivariable Windsurfer Adaptive Robust Control   100 200 300 400 500 600 100 200 300 400 500 600 Figure 5.10: Iterative system identification on each spatial frequency. 87 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 20 40 60 80 100 5 10 15 20 0 0.5 1 CD MD (a) Reference input. 20 40 60 80 100 5 10 15 20−0.2 0 0.2 0.4 0.6 0.8 CDMD (b) Initial response. 20 40 60 80 100 5 10 15 20 0 0.5 1 CDMD (c) Enhanced response. Figure 5.11: Response to a step reference input with CD square shape (n=20). 88 Chapter 5. Multivariable Windsurfer Adaptive Robust Control and dynamic frequencies) and with the size of steps in bandwidth increasing (), a different number of iterations might be needed to reach to a desired bandwidth and controller. Indeed, the number of iterations for system identification and controller design differs between spatial frequencies. The fewer the number of actuators means a greater impact of the circulant approximation. The effect of this perturbation on the system identification has been shown in Figure 5.12(a) and 5.12(c). Example 2 This example concerns the same process model as Example 1. This time a non-zero mean, square-shape disturbance is forced upon the system. This disturbance might be introduced by a variety of sources (e.g. The cross-coupling of actuation between processes might have the appearance of moisture actuation in the weight profile). It is expected that the controller would eliminate the effect of the disturbance. The presence of a disturbance illustrated in Figure 5.13 shows that even after 50 scans the initial controller cannot remove it. Figure 5.14 shows the 200th CD profile of the initial and enhanced controller. As can be seen even after 200 scans, the disturbance is still observable with the initial controller, while by applying the Windsurfing method after a few scans the disturbance is perfectly cancelled. 5.3 MultivariableH∞ Windsurfing In Section 5.2 the application of Iterative Windsurfer adaptive robust control to the continuous multivariable CD process of continuous web manufacturing was presented. The IMC method used for Windsurfing had the disadvantage that if the model has lightly damped stable/unstable poles and zeros within the closed-loop passband, then the design will be poor (e.g. the closed-loop gain will be very small near those zeros). In this Section, anH∞ design algorithm is introduced in order to make the iteration procedure more systematic and to remove the empirical aspect from the stopping criteria. The same assumptions that the large-scale CD processes can be approximated 89 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 10 20 30 40 50 60 70 80 90 100 −0.09 −0.08 −0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 CD(Space) Ex ci ta tio n Va lu e   Identified Model Actual Response (a) system with 100 actuators. (b) system with 50 actuators. (c) system with 11 actuators. Figure 5.12: Nominal and identified actuator response model of the basis weight. 90 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 20 40 60 80 100 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 1 CDMD (a) Additive disturbance. 20 40 60 80 100 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 1 CDMD (b) Initial response. 20 40 60 80 100 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 1 CDMD (c) Enhanced response. Figure 5.13: Response to a step disturbance input with CD square shape. 91 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 2 4 6 8 10 12 14 16 18 20 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 Actuator Number (CD) Ex cit at io n va lu e Enhanced response Initial response Figure 5.14: Disturbance rejection of the system after 200 scans. and diagonalized have been used to analyze the multivariable controller in terms of a family of Single-Input Single-Output (SISO) transfer functions across the spatial frequencies. Then in Section 5.3.2, we have applied discrete-timeH∞ Windsurfer adaptive robust control to each individual separated spatial frequency, starting with an initial model and a robust stabilizing controller at each spatial frequency. The modeling, design and implementation are the subject of this section. Finally, in Section 5.3.6, the designed multivariable iterative adaptive robust control technique is applied to a simulated industrial papermaking cross-directional control problem. Simulation results are given to illustrate the performance of the applied method. To conclude the results obtained from this method is also com- pared with the results obtained from classical Windsurfing. 5.3.1 Introduction Iterative adaptive robust control has emerged as a result of the extensive research on system identification and robust control design in early 1990’s and it has filled the wide gap between them [65]. One of the model-based iterative methods, the 92 Chapter 5. Multivariable Windsurfer Adaptive Robust Control windsurfing approach to adaptive robust control, was first introduced by Lee et al. in [115]. A few years later some issues of the proposed method were resolved by Lee et al. in [116–118]. The key element of the Windsurfing approach to adaptation is to increase the reference system bandwidth and to decrease the model uncertainty. The original Windsurfing method used the IMC design method in an iterative fashion. However, in cases where the identified model is unstable or has zeros on the jω-axis, the IMC design cannot be applied. Moreover, the IMC design method deals only with the output to reference input transfer function. For a sensible de- sign, the other three transfer functions (from plant input disturbance to output, from reference input to control signal and from plant disturbance to control signal) must be maintained below a certain size. In the following circumstances, the IMC method detailed in Section 5.2, cannot be properly used, or its application is limited by restrictive assumptions, or it may even fail to provide a good design [41–46]: 1. If the model has lightly damped stable poles within the closed-loop pass- band, then P (1 + PC)−1 will have large gain near the frequencies of those poles. 2. If the model has lightly damped stable/unstable zeros within the closed-loop passband, then either the closed-loop gain will be very small near those zeros or C(1 + PC)−1 will have large gain near the frequencies of those zeros. 3. If the bandwidth of F is chosen to be much larger than the bandwidth of P , then |C| will be very large at frequencies inside the bandwidth of F and outside the bandwidth of P , and again the control signal will be greatly am- plified at these frequencies. In transition from IMC design with its advantages and limitations to a robust control design, Dehghani et al. [41–46] have introduced an H∞ design algorithm as a modification to the Windsurfing algorithm. The main point of this section is to apply modified discrete H∞ Windsurfer adaptive robust control technique to a multivariable process with the aim to de- crease the uncertainties of given MD and CD model and achieve a better controller 93 Chapter 5. Multivariable Windsurfer Adaptive Robust Control design. We can achieve this goal by taking advantage of the H∞ properties of the modified Windsurfing. There are two reasons to make H∞ norm of a closed-loop transfer function small [186]: 1. To make the energy gain of the system small in order to satisfy performance objectives. Then different objectives in different frequency ranges is satisfied by using weighting functions. 2. To guarantee robustness in the face of modeling errors. 5.3.2 Multivariable Iterative Adaptive Robust Control Initially an Iterative Internal Model Control (IMC) method based on H2 optimal was applied in the control design step of the windsurfer approach; where the ref- erence input is a step function. See Chapters 4 and 7 of Reference [130] and Ref- erences [115–118] for more definitions and theorems. In our modified Windsurf- ing method [41–46] the control design is based on a H∞ admissible controller given an estimated model and desired closed-loop bandwidth λb. Satisfaction of the Theorems 5.2.1-5.2.3 allows us to apply Windsurfer Adaptive Robust Control on multivariable CD processes. In CD control design first low spatial and dynamical frequencies should be controlled, where the model uncertainty is very low. This control could easily be achieved, e.g., by feeding an initial controller and increasing resolution at every later iterations. The CD controller design and system identification in papermak- ing process can be explained as tuning of the existing controller and updating the available CD model to make the paper sheet as uniform as possible. This tuning can be summarized as increasing the bandwidth of excitation at each and every spatial frequency (see Figure 5.4). The same as the previous section for the Satisfaction of the initial robust stabil- ity at high spatial and dynamical frequencies a boundary Ωh is designed [164]: Ωh = {(ν, ω) : |G̃(ν, eiω)| < 0.1} (5.31) 94 Chapter 5. Multivariable Windsurfer Adaptive Robust Control where ν stands for the spatial frequency and ω stands for the temporal frequency. This boundary can be used for the desired bandwidth at each spatial frequency. Our goal is to increase the temporal bandwidth for each fixed spatial frequency. This is illustrated by the horizontal arrow in Figure 5.4. 5.3.3 H∞ Controller Design The desired filtering feature of the IMC is used in the H∞ design technique. The modified windsurfing satisfies two objectives: to have GiKi/(1 + GiKi) close to [Gi]aFi, and to satisfy a performance objective for three other transfer functions. 100 200 300 400 500 600 700 800 900 1000 50 100 150 200 250 300 350 400 450 500 Figure 5.15: Standard feedback configuration. [ y u ] = [ GiKi 1+GiKi Gi 1+GiKi Ki 1+GiKi 1 1+GiKi ] = H(Gi,Ki) [ r1 r2 ] (5.32) Then the admissible controller can be found by solving the followingH∞ prob- lem: γi = min Ki∈C ∥∥∥∥∥ GiKi1+GiKi − [Gi]aFi ε2(s) Gi1+GiKiε1(s) Ki1+GiKi ε1(s)ε2(s) 11+GiKi ∥∥∥∥∥ ∞ (5.33) where ‖G‖∞ = sup ω σ[G(jω)], C denotes the set of all proper stabilizing con- trollers for the plant, Gi(ν, z), [Gi]a is the all-pass part of the ith model, the fil- ter Fi(jω) can have an arbitrary roll-off rate, ε1(s) and ε2(s) are SISO, stable, minimum-phase and the proper weights. 95 Chapter 5. Multivariable Windsurfer Adaptive Robust Control The index in Equation 5.33 can be rewritten as γi = min Ki∈C ∥∥∥∥∥∥∥ Fl   −[Gi]aFi ε2(s)Gi0 ε1(s)ε2(s) Giε1(s) 1 −ε2(s)Gi −Gi  ,Ki  ∥∥∥∥∥∥∥ ∞ (5.34) where Fl(·, ·) denotes the lower Linear Fractional Transformation (LFT) of the generalized plant. Weight Design Design of the weighting functions is a part of the setting up of theH∞ performance index in 5.33 - 5.34. We will have different objectives in different frequency re- gions, depending on the application specifications and also the characteristics of the plant. Let us define our objectives, as specified in Equation 5.33: 1. Let the positive number α be the H∞ sense of desired closeness between GiCi(1 + GiCi)−1 and [Gi]aFi. That is we require ‖GiCi(1 + GiCi)−1 − [Gi]aFi‖∞ 6 α. 2. Let the positive number βc denote the required maximum tolerable gain in the appropriate frequency region for the transfer function Tur1 = Ci(1 + GiCi)−1. That is, we require the maximum singular value σ̄[Ci(1+GiCi)−1] 6 βc, ∀ω ∈ [ω1, ω2]. 3. Let the positive number βp denote the required maximum tolerable gain in the appropriate frequency region for the transfer function Tyr2 = Gi(1 + GiCi)−1. That is, we require the maximum singular value σ̄[Gi(1+GiCi)−1] 6 βp, ∀ω ∈ [ω3, ω4]. The three parameters, α, βc and βp, can be used to specify ε1(s) and ε2(s). It’s just needed to choose |ε1(s)| > α/βc near the frequencies of the lightly damped ze- ros of Gi and to choose |ε2(s)| > α/βp near the frequencies of the lightly damped poles of Gi. In other words γi 6 α implies that our three objectives are satisfied. 96 Chapter 5. Multivariable Windsurfer Adaptive Robust Control For details of the H∞ control problem, design of the weighting functions and the related theorems and discussions, the reader is referred to References [41– 46, 194, 195] DiscreteH∞ Design CD actuators have a continuous time response while the scanner scans the sheet of paper on a sampled signal (databoxes). The idealized model of the standard digital control system is shown in Figure 5.16. In the setup in Figure 5.16 the controller itself is K = HKdS. Figure 5.16: Standard sampled-data system. There are two groups (a total of three types) of Sampled Data (SD) controller design problems: direct method and indirect method [29]. Each type has its own advantages and disadvantages. We have chosen the direct approach (See Figure 5.17). The obvious advantage of this method is that it solves the problem with no approximations [29]. Discrete-time H∞-optimal controller design at each decomposed spatial fre- quency for the system shown in Figure 5.17 goes as follows: 97 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.17: Discrete-time system. Kd is to be designed. 1. Let the decomposed given state-model for Gd be ĝd(λ) = [ Ad Bd Cd Dd ] . (5.35) 2. Define the artificial continuous-time system Gc by ĝc(s) = [ Ac Bc Cc Dc ] , (5.36) where Ac = (Ad − I)(Ad + I)−1 Bc = (I −Ac)Bd Cc = Cd(Ad + I)−1 Dc = Dd − CcBd. (5.37) 3. Design an H∞(C+)-optimal controller Kc for Gc. Let us assume that the 98 Chapter 5. Multivariable Windsurfer Adaptive Robust Control state model of Kc is k̂c(s) = [ AKc BKc CKc DKc ] . (5.38) 4. Define the discrete-time system Kd by k̂d(λ) = [ AKd BKd CKd DKd ] , (5.39) where AKd = (I −AKc)−1(I +AKc) BKd = (I −AKc)−1BKc CKd = CKc(I +AKc) DKd = DKc + CKcBKd. (5.40) Invertibility of Ad + I and I − AKc is a necessary assumption for validity of the procedure. Since the bilinear transformation preservesH∞-norms, this method works (whereas it does not preserve H2-norms, for example) [29]. The detail and the proof of theorems behind the above-mentioned technique can be found in [29]. 5.3.4 System Identification Hansen’s indirect parametric closed-loop system identification [83, 84] has been adopted in Windsurfing. The a priori assumptions for modified H∞ Windsurfing are as follows [2]: Assumption 4 The plant is linear and time-invariant. Assumption 5 The plant is stabilized with a known initial controller. (This controller almost certainly has low bandwidth and low authority). Assumption 6 A reasonable transfer function model of the plant is available. 99 Chapter 5. Multivariable Windsurfer Adaptive Robust Control These assumptions are intentionally modest. There is no assumption on the accuracy of the initial model, the order of the model, the relative degree of the plant and whether the model is minimum-phase and/or stable or not [2, 41]. By using Hansen’s indirect closed-loop system-identification procedure em- bedded in the windsurfer approach, we can transform the problem of identifying G as a closed-loop system into an open-loop system-identification problem [116]. For this purpose we apply a system identification algorithm for a designed Box Jenkins model to obtain an estimate ofRj,i and Sj,i, which satisfies Hansen’s open- loop representation β = Rj,iα+ Sj,ie, (5.41) where α = Xj,ir1 + Yj,ir2, and β = Diy −Niu. The appropriate signal model for this system identification is shown in Figure 5.18. Figure 5.18: Identification of Rj,i and Sj,i [118]. Since the “input” α is independent of the measurement noise which affects the “output” β, identification of the stable transfer functions Rj,i and Sj,i is an open- 100 Chapter 5. Multivariable Windsurfer Adaptive Robust Control loop identification problem. A probing signal could be designed using dual adaptive control principles, e.g., by adding a filtered model uncertainty in the performance index to be minimized by the controller. This would reduce the uncertainty at frequencies higher than those currently in the bandwidth of the controller. Once the uncertainty is reduced to a satisfactory level, a new controller with a higher bandwidth can be designed and implemented for the next iteration. The above-mentioned iterative method has two different stopping criteria. Sat- isfaction of the first criterion provides us with the closed-loop performance objec- tives. The other stopping condition takes care of the quality of obtained controller and also evaluates the identified model. To validate a model of the plant it is re- quired that the closed-loop transfer function of the identified model and the con- troller, Gi, Ci, and the closed-loop transfer function of the plant and the controller, Gt, Ci, would be close in a frequency by frequency sense (i.e. max16k6N σ̄[H(Gt, Ci)− H(Gi, Ci)(jωk)] is small, where H(Gi, Ci) refers to the 4-block transfer function defined in Equation 5.32. Lemma 5.3.1 (Dehghani et al. [41]) For any Gt, Gi and Ci such that [Gt, Ci] and [Gi, Ci] are stable; that is, controller Ci stabilizes both Gt and Gi; then for SISO systems max 16k6N σ̄[H(Gt, Ci)−H(Gi, Ci)(jωk)] = = [ |Ci|+ 1|Ci| ] ∣∣∣∣ GtCi1 +GtCi − GiCi1 +GiCi ∣∣∣∣ (5.42) 5.3.5 Modified Multivariable Windsurfing Algorithm Our multivariable system identification and control design procedure can be sum- marized in the following algorithm. Step 1 Make the interaction matrix circulant and do the Spatial Frequency De- composition (SFD): G(z) ≈−→ Ĝ(z), 101 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Ĝ(z) SFD−−−→ G̃(ν, z). Step 2 For each decomposed spatial frequency set i = 1, j = 1 and Gi = G0, where G0 is the transfer function of an initial model of the plant. Step 3 Define the artificial continuous-time system Gc. Step 4 Find the critical frequency regions in Gi. The critical frequency regions are where the model has lightly damped poles. Step 5 Define positive numbers α, βc and βp, respectively, in the corresponding problematic frequency regions, as follows:∥∥∥ GiCi1+GiCi − [Gi]aFi ∥∥∥∞ 6 α, σ̄ [ Ci 1+GiCi ] 6 βc, ∀ω ∈ [ω1, ω2], σ̄ [ Gi 1+GiCi ] 6 βp, ∀ω ∈ [ω3, ω4]. Step 6 Design the frequency weights, ε1(s) and ε2(s): |ε1(s)| > α/βc, ∀ω ∈ [ω1, ω2], |ε2(s)| > α/βp,∀ω ∈ [ω3, ω4]. Step 7 Solve the H∞(C+) controller design problem given in Equation 5.33 and obtain γi and the admissible controller Kcj,i ; i.e. γi = min Ci∈C ∥∥∥∥∥ GiCi1+GiCi − [Gi]aFi ε2(s) Gi1+GiCiε1(s) Ci1+GiCi ε1(s)ε2(s) 11+GiCi ∥∥∥∥∥ ∞ . Step 8 Define the discrete-time controller Kj,i(ν, z), K̃(ν, z) = [ AKd BKd CKd DKd ] . 102 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Step 9 If max16k6N σ̄[H(Gt, Cj,i)−H(Gi, Cj,i)(jωk)] is still small and γi 6 α, then let j = j + 1, increase the bandwidth for Fj,i by 10% and go to Step 7. Step 10 If there was a previous controller, let i = i+ 1, then discard the last con- troller, perform indirect parametric closed-loop system identification to up- date the model and go to Step 7; else choose next possible spatial frequency and jump to Step 2. At each iteration, the balanced realization closed-loop model reduction method detailed in [137] should be employed to avoid ‘de- gree explosion’ in the identified model. Remark 3 At each iteration all the decomposed identified models and controllers should be reconstructed back to one multivariable controller and model. K̂(z) = F T K̃(ν, z)F, Ĝ(z) = F T G̃(ν, z)F. Without loss of generality the same model and technique can be used for mul- tivariable systems with lower or higher number of actuators. Intuitively speaking, the higher the number of the actuators, the narrower the additive band. The nar- rower the band, the smaller the size of perturbation introduced by the circulant extension. According to the current literature on iterative control, four to five iterations are often enough to produce acceptable performance. However, in this method steps of ten percent are chosen for the controller design. Therefore at least ten iterations are needed to reach the desired bandwidth of the controller. 5.3.6 Results This section demonstrates the effectiveness of the proposed method through two examples. Examples of the windsurfing methodology will be illustrated by ap- plying it to the control of a simulated CD multivariable system. The design steps outlined above are presented in terms of a realistic numerical example from the papermaking industry. 103 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Example 1 In this example the proposed method is used to control the basis weight of a sheet of paper which was explained in Example 5.25. Weight of the paper is typically controlled through an array of n=100 identically designed and distributed slice lip actuators. The same initial model G0(z) given in Equation (5.43) has been used for all spatial frequencies. This model mismodels the DC gain and fails to capture the response: G0(z) = 0.1 1− 0.9z−1 (5.43) The initial controller has the format described by Stewart [164] as follows: K(z) = [I + S(z)]−1C(z) (5.44) where C, S and K ∈ R100×100. The circulant extensions of the above-mentioned matrices are required to design the new controller. As discussed in Section 5.3.3 and in Steps 5 and 6 of our algorithm, we should set closed-loop performance objectives to be: (1) have closeness betweenGiCi/(1+ GiCi) and [Gi]aFi in anH∞ sense below 0.01 (α = 0.01); (2) Keep the maximum tolerable gain for the transfer function Tur1 below βc = 100(40 dB); and Tyr2 be- low βp = 10(20 dB). The frequency cost ε1(jω) is chosen to gradually reach the maximum gain of almost−80 dB (α/βc = 0.0001 = −80 dB) at 0.0042 rad/sec, where the plant model loses its bandwidth to the controller. The nature of actuator MD response has no lightly damped pole (in resonant mode), hence ε2(jω) can arbitrarily be set to be very small (Let it be α/βp = 0.001 = −60 dB). Balanced realization method is applied to do the model-order reduction. The result is truncated to retain all Hankel singular values greater than 0.01σ1. Figure 5.19 shows the result of iterative system identification on one spatial frequency. Iteratively, the same procedure has been done for all other spatial fre- quencies. Figure 5.20 shows the initial and updated step response of the discussed system to the reference input with spatial square shape. By comparing graphs (b) and (c) of 104 Chapter 5. Multivariable Windsurfer Adaptive Robust Control −30 −25 −20 −15 −10 −5 0 5 10 M a g n itu d e  (d B )   10−4 10−3 10−2 10−1 100 −720 −540 −360 −180 0 P h a se  (d eg ) Bode Diagram Frequency  (rad/sec) Plant Initial Model Updated Model 1 Updated Model 2 Figure 5.19: Iterative system identification on each spatial frequency. 105 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Figure 5.20, we observe that not only the MD response has been improved but also the resolution of the CD response has changed from an almost sinusoidal shape to a very high-resolution square response similar to the noiseless reference input shown in Figure 5.20(a). To obtain the given results shown in Figure 5.20, ideally 10-12 iterations are needed at the lowest spatial frequencies. However, at the high spatial frequencies just a few iterations (one iteration for the controller design and one iteration to test the new model) are required. Depending on the signal to noise ratio, model uncertainties both at high spatial and dynamic frequencies and size of steps in bandwidth increasing (), a different number of iterations might be needed to reach the same bandwidth and controller. Indeed, the number of iterations for system identification and controller design differs between spatial frequencies. Figure 5.21 shows the initial and updated step response of the discussed system to the reference input with spatial square shape and additive colored noise. Example 2 This example concerns the same process model as Example 1. However, this time a non-zero mean, square shape disturbance is forced into the system. We would expect the controller to remove the disturbance. The simulation results, presented in Figure 5.22 show that although 100 scans were made and a strong presence of disturbance is observed with an initial controller; the results with this updated controller have, nonetheless, been greatly improved. Figure 5.23 shows the 200th CD scanning databoxes of the initial and enhanced controller. As can be seen even after 200 scans noise components are still observ- able. While by applying the Windsurfing method after just a few scans all compo- nents of the noise are completely removed. Better CD elimination of the noise and faster response confirms our enhancement of the system identification and control design. 106 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 20 40 60 80 100 5 10 15 20 0 0.2 0.4 0.6 0.8 1 CD Reference Input MD (a) Reference input. 20 40 60 80 100 5 10 15 20 0 0.2 0.4 0.6 0.8 1 CD Initial Response MD (b) Initial response. 20 40 60 80 100 5 10 15 20 0 0.2 0.4 0.6 0.8 1 CD Enhanced Response MD (c) Enhanced response. Figure 5.20: Response to a step reference input with CD square shape under noise free conditions (n=20). 107 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 50 100 150 200 5 10 15 20 −0.5 0 0.5 1 1.5 CD Reference Input MD (a) Reference input. 0 50 100 150 200 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 CD Initial Response MD (b) Initial response. 0 50 100 150 200 5 10 15 20 −0.5 0 0.5 1 1.5 CD Enhanced Response MD (c) Enhanced response. Figure 5.21: Response to a step reference input with CD square shape corrupted with coloring noise (n=20). 108 Chapter 5. Multivariable Windsurfer Adaptive Robust Control (a) Additive disturbance. (b) Initial response. (c) Enhanced response. Figure 5.22: Rejection of the disturbance with CD picketing shape. 109 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 2 4 6 8 10 12 14 16 18 20 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 Actuator Number (CD) Ex cit at io n Va lu e   Enhanced Response Initial Response Figure 5.23: Disturbance rejection of the noise-free system after 200 scans. 5.3.7 IMC versusH∞-based Windsurfing The identified model from the Modified H∞-based multivariable windsurfing was compared to the results of IMC-based windsurfing (See Figure 5.24). As can be seen from Figure 5.25 the H∞-based Windsurfing gives a good re- sult at low frequencies, where our main interest of controller design. The H∞ controller design is more conservative then the IMC design. Therefore H∞ con- troller suppresses most of the high frequency components. Consequently during the system identification more conservative model will be found. It means that the identified model based onH∞ has less gain than the IMC at high frequencies. In a sense that at high frequencies model uncertainties are high, it is better to be more conservative at those frequencies and that is what theH∞ controller design does. IMC achieves exact reference tracking in desired bandwidth on the identified model. However, it is a the cost of creating large discrepancy defined in Equation 5.45. This is contrary to the H∞ controller design that ensures satisfaction of the defined threshold. 110 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 2 4 6 8 10 12 14 16 18 20 −0.09 −0.08 −0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 CD (Space) Ex ci ta tio n Va lu e IMC Windsurfing True Response Hinf Windsurfing Figure 5.24: Identified model by Windsurfing based on IMC versusH∞ result. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −120 −110 −100 −90 −80 −70 −60 −50 −40 −30 Normalized Frequency  (×pi rad/sample) Po we r/f re qu en cy  (d B/r ad /sa mp le) Power Spectral Density Estimate via Periodogram PSD Hinf Windsurfing Model PSD IMC Windsurfing Model True Model Figure 5.25: Power Spectral Density of the identified model by Windsurfing based on IMC versusH∞. 111 Chapter 5. Multivariable Windsurfer Adaptive Robust Control max 16k6N σ̄[H(Gt, Ci)−H(Gi, Ci)(jωk)] = = [ |Ci|+ 1|Ci| ] ∣∣∣∣ GtCi1 +GtCi − GiCi1 +GiCi ∣∣∣∣ (5.45) For the details on the advantages ofH∞-based Windsurfing method compared to the classical Windsurfing the interested reader is referred to [41–46]. 5.4 Cautious Bandwidth Control with ν-gap The application of classical and modified Windsurfer adaptive robust control to the continuous multivariable CD process of continuous web manufacturing was presented in Sections 5.2 and 5.3 respectively. In this section, we propose new enhancements on the iterative Windsurfing based on new advances in robust and optimal control. The Windsurfing approach to iterative adaptive robust control requires a series of controller designs with gradually expanding closed-loop bandwidth. Initially an IMC controller design was recommended to build a controller. As it was mentioned in the previous section the IMC was replaced with H∞ controller design. In this new approach to Windsurfing the bandwidth of the controller increases by 10% at each iteration. Any problematic frequency regions are defined based, on the application specifications and also the characteristics of the plant. As was mentioned in Sections 5.2.3 and 5.3.4, one of the main assumptions in either classical or modified-iterative Windsurfing was that the nominal model of the plant satisfied the Bezout Identity. In particular, it is desirable to find the minimal realization (the McMillan degree) of a plant in order to preclude the possibility of having hidden modes. This leads us to the concept of coprime factors. Let P = M−1N where M,N ∈ H∞ have been defined as the normalized left coprime factorization of the nominal system and P∆ be the transfer function of the perturbed plant respectively. The uncertainty can be modeled as an error of the 112 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 100 200 300 400 500 600 700 50 100 150 200 250 300 350 400 450 Figure 5.26: Normalized coprime factorization: C stabilizes P = M−1N . coprime factors (see Figure 5.26) [Description of the model uncertainty [144]] : P∆ = (M + ∆M )−1(N + ∆N ). (5.46) The set of transfer functions in Figure 5.26 from y and u to v1, v2, w1 and w2 can be written as: [ y u ] = [ P I ] (I−CP )−1 [ −C I ] [ v2 v1 ] + + [ I C ] (I − PC)−1 [ I −P ] [ w2 w1 ] (5.47) where the signals w1 and w2 represent disturbances acting on the plant and the sig- nals v1 and v2 represent noise on the compensator output and on the measurements respectively. Based on Equation 5.48 and since det(I − CP ) = det(I − PC), it is enough to only check either the first set or the second set of transfer functions. 113 Chapter 5. Multivariable Windsurfer Adaptive Robust Control [ P I ] (I−CP )−1 [ −C I ] [ v2 v1 ] + [ I C ] (I−PC)−1 [ I −P ] [ w2 w1 ] = [ I 0 0 I ] (5.48) Definition 1 (Admissible Perturbation [144]) For  > 0 a perturbation ∆ = [∆N ,∆M ] of a transfer function as described in [19] is called (−) admissible, when ‖∆‖∞ <  holds. The set of all (−) admissible perturbations is denoted by D: D := {∆ : ∆ ∈ RH∞; ‖∆‖∞ < }. (5.49) Theorem 5.4.1 (Robust Stabilization) C stabilizes the uncertain plantP∆ = (M+ ∆M )−1(N + ∆N ) shown in Figure 5.26 for all ∆ = [∆N ,∆M ] ∈ D iff 1. C stabilizes P 2. The following equation holds:∥∥∥∥∥ [ C(I − PC)−1M−1 (I − PC)−1M−1 ] ∥∥∥∥∥ ∞ 6 1/ . (5.50) This approach gives the ability to obtain the maximum stability margin in an explicit manner, compared to the H∞ loop shaping where a binary search is re- quired to arrive at the optimal value. However, this theorem is not valid for systems with high model-uncertainties, where ‖∆‖∞ > . 5.4.1 ν-gap In order to state sharper stability result in the face of plant uncertainty, Vinnicombe [184] has introduced a new metric. Definition 2 (Generalized Stability Margin) bP,C = { ‖T (P,C)‖−1∞ if (P,C) is stable 0 otherwise (5.51) 114 Chapter 5. Multivariable Windsurfer Adaptive Robust Control where T (P,C) = [ −P (I − CP )−1C P (I − CP )−1 −(I − CP )−1C (I − CP )−1 ] (5.52) = [ P I ] (I − CP )−1 [ −C I ] . (5.53) Note also that this function can also be written as bP,C = ‖T (P,C)‖−1∞ = 1/‖T (P,C)‖∞ = 1/sup ω σ̄(P,C)(jω) = inf ω σ(P,C)(jω). (5.54) A scalar measurement of the difficulty of controlling a particular plant is given by bopt(P ) = sup C bP,C , (5.55) It means that larger values of bopt(P ) mean − a better resolution of conflicting objectives of loop following, − a lower sensitivity to noise and disturbances, − an avoidance of input saturation, and − a low sensitivity to plant modeling inaccura- cies of certain types [4]. This introduction brings us to the idea of ν-gap. Vinnicombe [183, 186] introduced a ν-gap metric which defines a distance between any two linear time-invariant plants, stable or unstable, with the same input and output dimensions. Definition 3 [ν-gap metric] Define the ν-gap metric between two systems as fol- 115 Chapter 5. Multivariable Windsurfer Adaptive Robust Control lows: δν(P1, P2) := { ‖G̃∗2G1‖∞, if det(G̃∗2G1)(iω) 6= 0 and wno det(G̃∗2G1) = 0 1, otherwise (5.56) where ‖G̃∗2G1‖∞ = ‖(I+P ∗2P2)− 1 2 (P2−P1)(I+P ∗2P2)− 1 2 ‖∞ and wno(g) stands for the winding number about the origin of g(s) as “s” follows the standard Nyquist contour and it can be written as: wno(g) = ζ(g)− η(g), (5.57) it can also be written as wno det(G̃∗2G1) = wno det(I + P ∗ 2P1) + η(P1)− η̄(P2), (5.58) where η(g) denotes the number of open Right Half Plane (RHP) poles of g(s), ζ(g) denotes the number of open RHP zeros of g(s) and η̄(g) denotes the number of closed RHP poles of g. Satisfaction of the winding number condition is equivalent to the existence of a multivariable homotopy in the ν-gap metric between the two operators [18]. In contrast to the small-gain theorem, which relies on considering stable perturba- tions, homotopy arguments rely on perturbations which are bounded on the imagi- nary axis and also rely upon counting winding numbers [186]. It can be concluded that the ν-gap (δν(P1, P2)) is a general measure of the worst case difference in the performance of feedback system with P1 and feedback system with P2 respectively [185]. The ν-gap has the following unique properties that makes it a very valuable tool in control design. • It gives an “if and only if” theorem for robust stability and is easy to compute in contrast to theH2 gap. • It induces the correct graph-topology, in contrast to L2 gap. 116 Chapter 5. Multivariable Windsurfer Adaptive Robust Control • It is a closed-loop measure, in contrast to theH∞-norm. • It automatically emphasizes the cross over frequency. • The definition also holds in the discrete-time case by means of the bilinear transformation s = z−1z+1 . • κ(P1, P2) ∆= (I + P ∗2P2)− 1 2 (P2 − P1)(I + P ∗2P2)− 1 2 is called the chordal distance between P1 and P2 at frequency ω. Figure 5.27: Projection onto the Riemann sphere of the Nyquist plots of two trans- fer functions and chordal distance between them. In the SISO case [33], the ν-gap’s geometry of the chordal distance at fre- quency ω, κ(P1(jω), P2(jω)), is interpreted as the distance between the projec- tions onto the Riemann sphere of the points of the Nyquist plots of P1 and P2 corresponding to that frequency (See Figure 5.27). The stability results of ν-gap applicable to iterative system identification and control design are as follows: 117 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Proposition 5.4.2 (Stability Results) 1. Given a nominal plant P1 and a sta- bilizing compensator C, then (P2, C) is stable for all plants P2 satisfying δν(P1, P2) < β if and only if bP1,C > β 2. Given a nominal plant P and a stabilizing compensator C1, then (P,C2) is stable for all compensators C2 satisfying δν(C1, C2) < β if and only if bP,C1 > β. It can be concluded that any controller that satisfies the following constraint criterion without introducing caution would stabilize P. δν(C1, C2) < kbPi,Cj (5.59) where k is a safety constant in (0, 1) to consider the error between bP,Cj and bPi,Cj . Based on the assumption that Pi is a sufficiently accurate model, bP,Cj ≈ bPi,Cj , the stability of the actual closed loop follows from the part 2 of Proposi- tion 5.4.2. Whilst in Definition 3 it is not sufficient that δν(C1, C2) be small (i.e. C1 and C2 to be close enough that they simultaneously stabilize a suitably large set of controlled plant’s model); wno det(‖K̃∗2K1‖∞) = 0 satisfies that condition. Where ‖K̃∗2K1‖∞ = ‖(I + C∗2C2)− 1 2 (C2 − C1)(I + C∗2C2)− 1 2 ‖∞. 5.4.2 Iterative Controller Design In modified Windsurfing design scheme, during each iteration the nominal plant model of Gi is stabilized with the controller Cj . Then the bandwidth of the closed- loop response is increased by 10%. At each stage of the iteration two possible situations can arise [3]: 1. The modelGi is satisfactory, which means closed loop (G,Cj) and (Gi, Cj) are “close”. Then the controller adjustment based on Gi from Ci,j to Ci,j+1 is made. 2. ‖T (G,Ci,j) − T (Gi, Cj)‖∞ is large, which means that the model Gi is no longer satisfactory in the H∞ sense and we want to compute a new model that allows us to compute a better controller Cj+1. 118 Chapter 5. Multivariable Windsurfer Adaptive Robust Control By defining some constraints and knowing the problematic frequencies De- hghani et al. [41–46] brought some enhancements to the Windsurfing model. They recommended an increase in bandwidth at each iteration with static increments of ten percent. It is very common in occurrence (e.g., on paper machines) that the sys- tem suffers from high model uncertainties at high frequencies and low uncertainties at low frequencies. Therefore the recommended steps of the classical orH∞-based Windsurfing design method suggest an inadequate number of iterations. Our proposed approach in this section lets us increase the closed-loop band- width without going through the iterations of control design, which are needed in classical Windsurfing. Those iterations are waived while the ν-gap between the initial and updated controller satisfies robust stability margin. Consequently in- stead of increasing the bandwidth with static increments of ten percent we find the maximum admissible controller without going through the iterations. This problem can be resolved by considering part 2 of Proposition 5.4.2, in which theH∞Windsurfing algorithm explained in Subsection 5.3.5 would be mod- ified by the ν-gap metric satisfying three following criteria: • Performance objectives should be captured by a weighting function. • The ν-gap distance between updated controller and the initial controller, δν(Ci,j , Ci,j+1), should be small. • The H∞ design should be feasible for the designed weight. (maxC bPi,Ci,j should be sufficiently large. 5.4.3 Enhanced Multivariable Windsurfing Algorithm The modified Windsurfing Algorithm can be written following these steps: Step 1 Make the interaction matrix circulant and do the Spatial Frequency De- composition (SFD): G(z) ≈−→ Ĝ(z), Ĝ(z) SFD−−−→ G̃(ν, z). Step 2 For each decomposed spatial frequency set i = 1, j = 1 and Gi = G0, where G0 is the transfer function of an initial model of the plant. 119 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Step 3 Define the artificial continuous-time system Gc. Step 4 Find the critical frequency regions in Gi. The critical frequency regions are where the model has lightly damped poles. Step 5 Define positive numbers α, βc and βp, respectively, in the corresponding problematic frequency regions, as follows:∥∥∥ GiCi,j1+GiCi,j − [Gi]aFi,j ∥∥∥∞ 6 α, σ̄ [ Ci,j 1+GiCi,j ] 6 βc, ∀ω ∈ [ω1, ω2], σ̄ [ Gi 1+GiCi,j ] 6 βp, ∀ω ∈ [ω3, ω4]. Step 6 Design the frequency weights, ε1(s) and ε2(s): |ε1(s)| > α/βc,∀ω ∈ [ω1, ω2], |ε2(s)| > α/βp, ∀ω ∈ [ω3, ω4]. Step 7 Solve the H∞(C+) controller design problem and obtain γi,j and the ad- missible controller Kci,j ; i.e. γi,j = min Ci,j∈C ∥∥∥∥∥ GiCi,j 1+GiCi,j − [Gi]aFi,j ε2(s) Gi1+GiCi,j ε1(s) Ci,j 1+GiCi,j ε1(s)ε2(s) 11+GiCi,j ∥∥∥∥∥ ∞ Step 8 Let j = j + 1, increase the bandwidth for Fi,j by 10% and solve the H∞(C+) controller design problem and obtain γi and the admissible mul- tivariable controller Kci,j . Find the ν-gap between the last stabilizing con- troller and the last designed (i.e. δν(Ci,j , Ci,j+1)) and find the generalized stability margin, bPi,Ci,j . While the following inequality is satisfied repeat this step. δν(Ci,j , Ci,j+1) < kbPi,Ci,j 120 Chapter 5. Multivariable Windsurfer Adaptive Robust Control Step 9 Discard the last controller, i.e. Ci,j+1 (The one that didn’t satisfy the above-mentioned inequality.). Apply the last-accepted controller, Ci,j , to the plant. If max16k6N σ̄[H(G,Ci,j)−H(Gi, Ci,j)(jωk) is still small and γi,j 6 α, go to Step 8. Step 10 If there was a reconstructed controller, then discard it. Let i = i + 1, then perform indirect parametric closed-loop system identification to update the model and go to Step 7; else choose next possible spatial frequency and jump to Step 2. Remark 4 At each iteration in Step 10, the balanced realization closed-loop model reduction method detailed in [137] should be employed to avoid a degree explosion in the identified model. Remark 5 At each iteration all the decomposed identified models and controllers should be reconstructed back to one multivariable controller and model. K̂(z) = F T K̃(ν, z)F, Ĝ(z) = F T G̃(ν, z)F. 5.4.4 Example Our proposed method can be applied to the simulated paper machine detailed in Section 5.3.6. In this experiment the size of  was chosen to be small enough  = 0.002 to have a cautious bandwidth increase and controller design at high fre- quencies, while the ν-gap eliminates any unnecessary iterations at low frequencies. The identified CD response of the actuator is shown in Figure 5.28. It is observed that applying ν-gap on theH∞ Windsurfing reduces the number of iterations before we reach the final enhanced multivariable controller. While to design anH∞-based Windsurfing controller, based on the defined step-size and for all of the controllable spatial frequencies takes 653 iterations, we require just 43 iterations of our CD system, based on the ν-gap, to achieve our objective. The simulation results, presented in Figure 5.29 show how our enhanced Wind- surfing based on the ν-gap improves the global response of the CD system. 121 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 2 4 6 8 10 12 14 16 18 20 −0.09 −0.08 −0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 CD (Space) Ex ci ta tio n Va lu e   ν−gap Nominal Model Figure 5.28: Identified model by enhanced Windsurfing based on ν-gap. 5.5 Input Design There are different characteristics that should be considered and/or satisfied in an experiment design. First of all the “probing signal” should be rich [121]. It should excite the system and force it to demonstrate or display its possible properties of interest. An experiment is called informative enough if it generates informative data set. 5.5.1 Sinusoid on CD System identification of CD processes is usually done by bumping a few actuators across the sheet. The spectrum of the bump test has a high frequency CD profile. On the other hand the proposed Multivariable Iterative Adaptive Robust control design uses SFD and excites just one spatial frequency at each iteration. Conse- quently, instead of bumping actuator(s) it is suggested that the actuators be excited with a sinusoidal pattern at appropriated power level and spatial frequency of inter- est. This means the excitation should be carried out sequentially and in order and at the appropriate spatial frequencies. 122 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 20 40 60 80 100 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 CD Reference Input MD (a) Reference input. 20 40 60 80 100 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 1 CD Initial Response MD (b) Initial response. 20 40 60 80 100 5 10 15 20 −0.2 0 0.2 0.4 0.6 0.8 CDMD (c) Enhanced response. Figure 5.29: Enhanced Windsurfing response to a step reference input with CD square shape corrupted with coloring noise (n=20). 123 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 5.5.2 Chirp Signal on MD Hansen’s [83] framework of indirect closed-loop system identification is used in Windsurfing approach. It is already common to enter a step input at each iteration, either identification or control design. The dashed line in Figure 5.30 shows that the step signal has high power density at very low frequencies and doesn’t have good excitation at high frequencies. This means that we won’t have enough excitation at frequencies that interest us (e.g. cut-off frequency of the filter), where the model has high uncertainty and we want to identify a better model or we must increase the closed-loop bandwidth during iterations of the controller design. 0 0.5 1 1.5 2 2.5 3 10−8 10−6 10−4 10−2 100 102 104 Frequency (rad/sec) PSD Welch   Chirp Step Figure 5.30: PSD for chirp and step inputs. While step input has clear excitation at zero frequency, the chirp excites the entire design bandwidth almost with equal power. It was discussed in Section 5.2.3 that the estimation of the closed-loop dy- namics in Hansen’s framework was based on an associated open-loop estimation problem using dual Youla parametrization. Therefore an input signal should be chosen to satisfy the input design criterion for the open-loop system of Equation 5.41. Since either α or α are filtered known part of the input signal r, the input 124 Chapter 5. Multivariable Windsurfer Adaptive Robust Control design can be considered for r. 5 10 15 20 25 30 35 40 45 50 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Time(scans) Figure 5.31: An example of chirp signal. To excite all frequencies of the closed-loop bandwidth at each iteration that includes every spatial frequency, a Chirp signal is chosen as an input (See Figure 5.31). A Chirp signal is a sinusoid with a frequency that changes continuously over a certain band Ω: ω1 6 ω 6 ω2 over an iteration with time period 0 6 t 6M : r(t) = A cos(ω1t+ (ω2 − ω1)t2/(2M)) (5.60) This signal has the same crest factor as a pure sinusoid, i.e. √ 2, and it gives good control over the excited frequency bandwidth. Figure 5.30 shows that the spectrum of the chirp signal covering the frequency band 0 6 ω 6 0.5. ω2 should be chosen based on the size of the designed closed-loop bandwidth. C2r = maxt u2(t) limN→∞ ∑N t=1 u 2(t) (5.61) 5.5.3 Example A Chirp signal is applied on a simulated paper machine, which was explained in detail in Section 5.3.6. As can be seen from Figure 5.32 the result is compared 125 Chapter 5. Multivariable Windsurfer Adaptive Robust Control with the identification based on the step input. 5 10 15 20 25 30 35 40 45 50 −0.09 −0.08 −0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 0.01 CD (Space) Ex ci ta tio n Va lu e   True Response Identified Model Chirp Identified Model Step Figure 5.32: CD model identification by chirp and step-input. While step input creates a model with corrupted corners, the chirp input generates a better response with higher resolution. We observed that applying chirp signal as the input of the multivariable Wind- surfing gives better identification for each spatial frequency and consequently re- duces the CD model uncertainties. In Figure 5.33 the black line stands for the true model, the dash line corresponds to the model identified by chirp input and the doted line shows the model identified by step input. It is obvious that in most low frequencies the model identified by the chirp signal is closer to the true model than the model identified by the step input. 5.6 Conclusions The main section of this chapter has presented the framework and constructive technique for the system identification and iterative controller design on spatially distributed paper-making processes. The design technique incorporates the Spa- tial Frequency Decomposition method to convert MIMO system to the family of 126 Chapter 5. Multivariable Windsurfer Adaptive Robust Control 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −110 −100 −90 −80 −70 −60 −50 −40 −30 Normalized Frequency  (×pi rad/sample) Po w er /fr eq ue nc y (d B/ ra d/s am ple ) Power Spectral Density Estimate via Periodogram   Ideal Plant Chirp Input Step Input Figure 5.33: Spectrum of the model identified by chirp input versus step input. independent SISO systems. Then each separated, spatial-frequency classical IMC- based Windsurfing Adaptive Robust Control is used to do the system identification and control design. This approach is insensitive to model uncertainties at low fre- quencies. In conclusion we would like to emphasize that this approach allows a significant reduction in spatial and temporal model uncertainties with increased ro- bust performance without losing robust stability of the multivariable closed-loop system. Section 5.3 has presented the modified H∞-based Windsurfing on CD pro- cesses. This iterative method is used to create a system identification and control design package. In contrast to IMC-based Windsurfing the controller designed by H∞-based Windsurfing is more practical and cautious. In Section 5.4 of this Chapter, a novel technique for the improvement of up- dating CD controllers and increasing their closed-loop bandwidth, based on the ν-gap method has also been developed. The homotopy structure of ν-gap makes it a suitable method to decrease the number of iterations, before we reach the last possible designed controller. A closed-loop simulation example has also been pre- 127 Chapter 5. Multivariable Windsurfer Adaptive Robust Control sented and has shown the performance of our method compared to the H∞-based multivariable Windsurfing. Finally, suitable input signals on MD and CD directions have also been de- signed for the iterations of controller design. The advantage of choosing a sinusoid on CD is that the spatial frequencies are excited one at the time. Our input sig- nal is compared to the bump test which has high gain on all frequencies of the spectrum. Moreover, in an iteration of system identification or control design, the remaining spatial frequencies remain unexcited, making our method less intrusive. We recommend a chirp signal on MD. The advantage of applying a chirp signal on a system controlled by a Windsurfing method is that we can excite all of the designed closed-loop bandwidth with almost the same power. This is a noteable improvement when compared to the regular bump test that has moderate DC gain and very weak spectrum. 5.7 Notation i iteration number of the system identification j iteration number of the controller design ν spatial frequency number 1 6 ν 6 n ω dynamical frequency number 1 6 ω 6 n n number of actuators G0 interaction matrix Ĝ0 circulant extension of the interaction matrix G̃0 decomposed interaction matrix Kj,i jth controller for the ith model K̃(z) controller matrix K̂(z) multivariable controller (reconstructed) λj,i bandwidth of the controller 128 Chapter 6 Summary and Future Work This thesis has concentrated on the problems involved with modeling, alignment, system identification and design of feedback control design for dynamical systems of spatially distributed and controlled by an array of actuators. Paper machine is an industrial example of such a process and is the central role of this thesis. This dissertation consists of two parts: The first part focuses on the analysis and cross- directional alignment or mapping problem of the paper produced by these indus- trial paper machines. The second part has concentrated on the modeling, analysis, system identification and the design of an iterative cross-directional controller for these machines. 6.1 Summary In Chapter 1 the structure of the broad range of industrial paper machines is ex- plained. Depending on the type of actuator, the width of one actuator’s spatial response varies between two ‘actuator-width’ to about one third of the full-width of the paper sheet and also the actuator’s cross directional response can have Gaus- sian, Mexican-hat or bimodal shape. It is illustrated in Chapter 2, that paper machine CD control systems belong to a set of large-scale, multivariable, spatially-distributed, uncertain and ill-conditioned control systems. Separability of Machine Directional (MD) from the Cross Direc- tional (CD) models is a common assumption and has been used in this work. In Chapter 3 a review of the state-of-the-art on CD mapping of the paper ma- chine’s actuators to their corresponding scanned databoxes was presented. It is well-accepted that variations in the shrinkage of paper are one of the main causes of misalignment. Therefore, a vast academic work have been done on the iden- tification of shrinkage. Based on different approaches to mapping, the alignment 129 Chapter 6. Summary and Future Work methods have been classified into the following categories: open-loop or closed loop and intrusive or nonintrusive. Also new insight into the actuator-mapping by applying image processing based techniques was discussed. The cancellation of the correlation between input disturbance and the actuator setting is presented as recent developments on the closed-loop and non-invasive alignment methods. As illustrated in a simulation example in Chapter 3 these methods are based on re- strictive assumptions. They are also known as batch algorithms, which means quite number of samples are needed before the method produces satisfactory results. Considering that shrinkage is a nonlinear phenomenon happening in the drying section, it is assumed that the paper will develop new dimensions due to shrinkage. Two batch algorithms were discussed in the early sections of Chapter 4 and their outcomes and disadvantages are discussed. Then a new tensor model is designed for the CD process by having the number of actuators, time, scanning databoxes and the new dimension of the sheet before shrinkage in Chapter 4. Next the Paral- lel Factor decomposition (PARAFAC) method was used in an alternating fashion to decompose the tensor model. We demonstrated that a method of actuator re- sponse’s system identification and mapping diagram can successfully be achieved using our newly proposed approach. This work constitutes the first known application of Windsurfing Adaptive Robust Control to a multivariable industrial process. The iterative Windsurfing method consists of iterations of an IMC based controller design and Hansen’s indi- rect parametric closed-loop system identification. At each iteration the bandwidth of the closed-loop system is increased. The examples of Section 5.2 have shown that this method reduces model uncertainties by exciting higher bandwidth. It also brings some degree of robustness. Modification of the IMC-based control design of multivariable Windsurfing with H∞ technique is presented in Chapter 5.3. Relevant control engineering cri- teria (closed-loop stability, performance and robustness) are taken into account by the H∞-based Windsurfing method. Closed-loop performance objectives were set in the design of the weighting functions of the controller. Two examples demon- strated that application of our modified Windsurfing method has improved both the rejection of disturbance and process response to a reference input. Most of the model requirements mentioned in Section 5.2, have thus been removed with this 130 Chapter 6. Summary and Future Work approach. The constraints of our current Windsurfing method was discussed in Chapter 5.4. In classical and H∞-based Windsurfing method the bandwidth of the con- troller increases by arbitrary and constant steps with each iteration. This may cause unnecessary iterations of controller design at low frequencies where the process gain is high. To solve this issue a new enhancement of the Windsurfing based on the ν-gap has been proposed. Depending on the number of CD actuators, this tech- nique reduces the number of on-line iterations while maintaining close-to-optimal robust performance and stability. The ν-gap metric can be easily calculated directly from frequency responses, which means that the distance between two controllers can be estimated directly from frequency response measurements. It provides a homotopy of controllers that stabilize the same model of the plant. We provided an example demonstrated that the ν-gap technique decreases the number of iterations. An input signal is also designed to be used instead of a step-input at the itera- tions of control design and system identification. Hansen’s framework of indirect parametric closed-loop system-identification provides us with the ability to identify the closed-loop system through an open-loop procedure. To excite all the frequen- cies of the closed-loop bandwidth at each iteration and for every spatial frequency, a Chirp signal is chosen as an input. In order to decrease the amount of perturba- tions and not to excite the remaining spatial frequencies, while we were designing a new controller to handle the specific spatial frequency, a sinusoidal signal was chosen as the input. One example in Chapter 5.4 demonstrated that the designed input signal is significantly more efficient than the conventional bump tests (step response) and still perturbing the system, but with much less power. 6.2 Future Work Some of the possible research directions directly related to the work presented in this thesis are outlined below: The Alternating Least Square (ALS) methodology is currently accepted as be- ing valid in PARAFAC decomposition. Obviously, a more efficient iterative tech- nique would allow for even better alignment and system identification. In par- ticular, the use of better optimized methods in conjunction with the tensor model 131 Chapter 6. Summary and Future Work plus PARAFAC decomposition should result in better CD mapping and modeling. Many algorithms have been proposed to fit PARAFAC models [51]: • Penalty Diagonalization Error • Gauss-Newton • Preconditioned Conjugate Gradients • Levenberg-Marquardt • Direct Trilinear Decomposition • Alternating Least Squares • Pseudo-Alternating Least Squares • Alternating Trilinear Decomposition • Self-Weighted Alternating TriLinear Decomposition • Alternating Slice-wise Decomposition In this thesis a three-way model of PARAFAC was applied to the distributed system of paper-making. PARAFAC can be written as an ‘n-way model’ for multi- variable systems. Therefore, its n-way modeling can be applied to paper machines, where the cross-coupling of different array of actuators is considered as a consid- erable disturbance. Theoretical issues surrounding the Windsurfing are yet to be addressed. The controller design of classical Windsurfing contains the parameter . The value of  must be made sufficiently small to satisfy robust stability. Although the modi- fied H∞ Windsurfing method has improved this matter, no work has been done to determine an optimal value for . In this dissertation the performance of the Windsurfing method on CD pro- cesses is proven through different simulations. The proposed method should be applied to a real paper machine and some mill trials should be done to justify our simulation results. 132 Chapter 6. Summary and Future Work Finally, the Windsurfing approach to iterative adaptive robust control can be combined with non-causal transfer function of CD processes. As it becomes more accepted to control CD systems in steady-states conditions, Windsurfing can be used to increase the bandwidth of the spatial frequencies. 133 Bibliography [1] J. Adamy. Adaptation of cross-direction basis-weight control in paper ma- chines using fuzzy decision logic. International Journal of Approximate Reasoning; Elsevier Science Inc., 16(1):25–42, 1997. [2] B. D. O. Anderson. Windsurfing approach to iterative control design, chap- ter 7, in Iterative Identification and Control: Advances in Theory and Appli- cations, pages 142–166. Springer-Verlag, Germany, 2002. [3] B. D. O. Anderson, X. Bombois, M. Gevers, and C. Kulcsár. Caution in iter- ative modeling and control design. In Adaptive systems in control and signal processing, pages 13–19, Workshop, Glasgow, Royaume, August 1998. [4] B. D. O. Anderson and M. Gevers. Fundamental Problems in Adaptive Con- trol, chapter in Perspectives in Control, Theory and Applications. Springer- Verlag New York, Inc., 1998. [5] B. D. O. Anderson and R. L. Kosut. Adaptive Robust Control: On-Line Learning. In Proceedings of the 30th Conference on Decision and Control, pages 297–298, 1991. [6] B. D. O. Anderson and Y. Liu. Controller reduction: concepts and ap- proaches. IEEE Transactions on Automatic Control, 34:802–812, 1989. [7] J. U. Backstrom, C. Gheorghe, G. E. Stewart, and R. N. Vyse. Constrained model predictive control for cross-directional multi-array processes. Control Systems Conference, 102(5):33–36, 2000. [8] B. Bamieh, F. Paganimi, and M. Dahleh. Distributed control of spatially invariant systems. IEEE Transactions on Automatic Control, 47(7):1091– 1107, 2002. 134 Bibliography [9] J. Berggren, A. Zehnpfund, D. Conway, J. Vasel, and S. C. Chen. Multiple- input and multiple-output CD control applications. In Proceedings of the TAPPI papermakers and PIMA international leadership conference, page 18p, Jacksonville, FL, 2007. [10] J. Berggren, A. Zehnpfund, D. Conway, J. Vasel, and S. C. Chen. Multiple- input and multiple-output CD control applications - practical aspects. In Proceedings of the TAPPI papermakers and PIMA international leadership conference, page 14p, Jacksonville, FL, 2007. [11] W. L. Bialkowski. Newsprint uniformity and process control - the untapped competitive edge. In In Proceedings of 76th CPPA Annual Meeting, Techni- cal Section, pages 255–269, Montreal, Canada, 1990. [12] J. Borch. Optical and Appearance Properties, volume 1, chapter 4, in Hand- book of Physical Testing of Paper, pages 95–148. Wiley Interscience, New York, 2nd edition, 2002. [13] J. Borch, R. E. Mark, and M. B. Lyne. Handbook of physical testing of paper, volume 2. CRC Press, 2nd edition, 2001. [14] R. D. Braatz. Identification, estimation and control of sheet and film pro- cesses. In Proceedings of the 13th World IFAC Congress, volume 7, pages 319–324, San Francisco, California, June 1996. [15] C. E. Brandon. Pulp and paper chemistry and chemical technology, vol- ume 3, chapter in Properties of Paper, pages 1715–1972. Wiley Interscience, New York, third edition, 1972. [16] W. Brecht and Wanka. Neue beiträge zur kenntnis der dimensionsstabilität von paieren. Des papier, 21(7):354–360, 1967. [17] D. P. Brewster. CD controls have unique problems: they need unique trou- bleshooting tools. In TAPPI Proceedings: Joint Conference Process Con- trol, Electrical & Information Conference, pages 215–228, Vancouver, BC, Canada, 1998. 135 Bibliography [18] T. S. Brinsmead and B. D. O. Anderson. Homotopy for the ν-gap. Commu- nications in Information and Systems, 1(4):333–380, December 2001. [19] M. Cantoni. Gap metric performance bounds for linear feedback systems. In Proceedings of the 38th Conference on Decision and Control, pages 4505– 4510, Phoenix, AZ, USA, December 1999. [20] J. D. Carroll and J. Chang. Analysis of individual differences in multi- dimensional scaling via an n-way generalization of Eckart-Young decom- position. Psychometrika, 35(3):283–319, 1970. [21] J.D. Carroll, S. Pruzansky, and J.B. Kruskal. Candelinc: A general approach to multidimensional analysis of many-ways arrays with linear constraints on parameters. Psychometrika, 45(3), 1980. [22] S. C. Chen. Actuation cell response and mapping determinations for web forming machines. Technical Report 5122963, United States Patent, 1992. [23] S. C. Chen. Full-width sheet property estimation from scanning measure- ments. In Proceedings of the Control Systems Conference, Whistler, Canada, 1992. [24] S. C. Chen. Analysis of two-dimensional sheet variations from a non- scanning measurement. In Preprint of IFAC Workshop on Adaptive Systems in Control and Signal Processing, pages 369–374, Glasgow, Scotland, UK, 1998. [25] S. C. Chen and R. Subbarayan. Cross-machine direction (CD) response modeling with two-dimensional sheet variation measurements. In Proceed- ings of the American Control Conference, pages 3082–3086, San Diego, California, 1999. [26] S. C. Chen and R. Subbarayan. Identifying temporal and spatial responses of cross machine actuators for sheet-forming processes. IEE Proceedings Control Theory Application, 149:975–995, 2002. 136 Bibliography [27] S. C. Chen, R. Subbarayan, K. Kristinsson, and R. Snyder. Paper machine applications with fullsheet imaging measurement. In Proceedings of the Control Systems Conference, pages 330–337, Porvoo, Finland, 1998. [28] S. C. Chen and P. Q. Tran. Automated optimization of cross machine di- rection profile control performance for sheet making processes. Technical Report US006564117B1, United States Patent, 2000. [29] T. Chen and B. Francis. Optimal sampled-data control systems. Springer- Verlag London Ltd., 1996. [30] Z. P. Chen, H. L. Wu, J. H. Jiang, Y. Li, and R. Q. Yu. A novel trilinear decomposition algorithm for second-order linear calibration. Chemometrics and Intelligent Laboratory Systems, 52:75–86, 2000. [31] A. Cichocki and S. Amari. Adaptive Blind Signal and Image Processing. John Wiley and Sons LTD, 2002. [32] R. P. A. Constantino, S. J. I’Anson, and W. W. Sampson. The Effect of Machine Conditions and Furnish Properties on Paper CD Shrinkage Profile. In Fundamental Research Symposium, Cambridge, volume 13, pages 283– 306, September 2005. [33] B. Cordons. Process modelling for control. Springer-Verlag London Ltd., 2005. [34] K. W. Corscadden and S. R. Duncan. Second order optimization methods for use with an artificial neural network estimator used in a plastic extrusion process. Transactions of the Inst of Chemical Engineers, Rugby, England, 74(A6):689–694, 1996. [35] K. W. Corscadden and S. R. Duncan. Convergence analysis of a re- duced model estimator applied to an industrial plastic film extrusion text bed. In Proceedings of the 11th IFAC Symposium (SYSID), pages 171–176, Fukuoka, Japan, July 1997. 137 Bibliography [36] K. W. Corscadden and S. R. Duncan. Reduced-order estimator for closed- loop online estimation of cross-directional parameters in a plastics extru- sion process. IEE Proceeding Control Theory Application, 144(6):549–557, 1997. [37] H. Corte. Perception of the optical properties of paper, chapter in The Fun- damental Properties of Paper Related to its Uses. British Paper and Board Industry Federation, 1976. [38] P. J. Davis. Circulant matrices (Pure & Applied Mathematics). New York: Wiley, 1979. [39] A. L. F. de Almeida, G. Favier, and J. C. M. Mota. Tensor decomposition and applications to wireless communication systems, chapter in Telecom- munications: Advances and trends in transmission, networking and appli- cationsking and applications, pages 57–82. University of Fortaleza Press, 2006. [40] A. L. F. de Almeida, G. Favier, and J. C. M. Mota. PARAFAC-based unified tensor modeling for wireless communication systems with application to blind multiuser equalization. Signal Processing, 87:337–351, 2007. [41] A. Dehghani, A. Lanzon, and B. D. O. Anderson. An H∞ algorithm for the windsurfer approach to adaptive robust control. International Journal of Adaptive Control and Signal Processing, 18:607–628, 2004. [42] A. Dehghani, A. Lanzon, and B. D. O. Anderson. The windsurfer approach to iterative control and identification using anH∞ algorithm. In IFAC work- shop on adaptation and learning in control and signal processing, and IFAC workshop on periodic control systems, pages 777–782, Yokohama, Japan, August 2004. [43] A. Dehghani, A. Lanzon, and B. D. O. Anderson. AnH∞ model referencing design utilizing a two degree of freedom controller scheme. In Proceedings of 44th IEEE conference on decision and control, and European control conference, pages 2302–2307, Seville, Spain, December 2005. 138 Bibliography [44] A. Dehghani, A. Lanzon, and B. D. O. Anderson. A newH∞ design method for high performance robust tracking control. In Proceedings of 16th IFAC world congress, Prague, Czech Republic, July 2005. [45] A. Dehghani, A. Lanzon, and B. D. O. Anderson. H∞ design to generalize internal model control. Automatica, 42:1959–1968, 2006. [46] A. Dehghani, A. Lanzon, and B. D. O. Anderson. A two-degree-of-freedom H∞ control design method for robust model matching. International Jour- nal of Robust and Nonlinear Control, 16:467–483, May 2006. [47] J. C. Doyle and G. Stein. Multivariable feedback design: Concepts for a classical/modern synthesis. IEEE Transactions on Automatic Control, 26(1):4–16, February 1981. [48] S. R. Duncan. Cross-directional control of web forming processes. PhD thesis, Imperial College of Science Technology and Medicine, University of London, London, UK, 1989. [49] S. R. Duncan. Dynamic Modeling of Cross-Directional Actuators: Implica- tions for Control. In Proceedings of the Control Systems Conference, pages 19–22, May 1996. [50] S. R. Duncan, J. M. Allwood, W. P. Heath, and K. W. Corscadden. Dynamic modeling of cross-directional actuators: Implications for control. IEEE Transactions on Control Systems Technology, 8(4):667–675, July 2000. [51] N. M. Faber, R. Bro, and P. K. Hopke. Recent developments in CANDE- COMP/PARAFAC algorithms: a critical review. Chemometrics and Intelli- gent Laboratory Systems, 65(1):119–137, January 2002. [52] J. Fan. Model Predictive Control for Multiple Cross-Directional Processes: Analysis, Tuning, and Implementation. PhD thesis, University of British Columbia, September 2003. [53] J. Fan, G. E. Stewart, and G. A. Dumont. Two-dimensional frequency analysis for unconstrained model predictive control of cross-directional pro- cesses. Automatica, 40:1891–1903, 2004. 139 Bibliography [54] F. Farahmand, G. A. Dumont, M. S. Davies, and P. D. Loewen. Online identification and alignment of MIMO cross directional controlled processes using second order statistics. In Proceedings of the IEEE International Con- ference on Control and Automation (ICCA), pages 989 – 994, Guangzhou, China, May-June 2007. [55] F. Farahmand, G. A. Dumont, and P. D. Loewen. The windsurfer approach to adaptive robust identification and control in cross- directional processes. In The IEEE Advanced Process Control Applications for Industry Workshop, Vancouver, BC, Canada, May 2007. [56] F. Farahmand, G. A. Dumont, P. D. Loewen, and M. S. Davies. Iterative adaptive robust control of multivariable CD processes. International Journal of Adaptive Control and Signal Processing, (Submitted), 2009. [57] F. Farahmand, G. A. Dumont, P. D. Loewen, and M. S. Davies. Tensor- based blind alignment of MIMO CD processes. Journal of Process Control, 19(5):732–742, May 2009. [58] F. Farahmand, B. Gopaluni, M. S. Davies, P. D. Loewen, and G. A. Dumont. Autonomous alignment of CD control on paper machines. In Proceedings of the Control Systems Conference, pages 209–214, Tampere, Finland, June 2006. [59] A.P. Featherstone and R.D. Braatz. Input design for large-scale sheet and film processes. Industrial & Engineering Chemistry Research, 37(2):449– 454, 1998. [60] A.P. Featherstone, J.G. VanAntwerp, and R.D. Braatz. Identification and Control of Sheet and Film Processes. Springer-Verlag London Ltd., 2000. [61] C. Fu and S. Nuyan. CD control sensitivity analysis. In Proceedings of the American Control Conference, pages 3087–3091, San Diego, 1999. [62] C. Fu, S. Nuyan, and S. Bale. CD control performance assessment. In Proceedings of the Control Systems Conference, pages 289–297, Porvoo, Finland, September 1998. 140 Bibliography [63] W. Gallay. Stability of dimensions and form of paper (part 1). Journal of Tappi, 56(11):54, 1973. [64] S. Gendron and J. Hamel. On paper streaks caused by mapping misalign- ment of the cross-direction control system. In Proceedings of the Control Systems Conference, pages 117–122, Stockholm, Sweden, 2002. Swedish Pulp and Paper Research Institute. [65] M. Gevers. A decade of progress in iterative process control design: from theory to practice. Journal of Process Control, 12:519–531, 2002. [66] J. Ghofraniha. Cross-directional response modelling, identification and con- trol of the dry weight profile on paper machins. PhD thesis, The University of British Columbia, Department of Electrical and Computer Engineering, Vancouver, Canada, 1997. [67] J. Ghofraniha, M. S. Davies, and G. A. Dumont. Cross direction response identification and control of paper machine using continous wavelet trans- form. In Proceedings of the American Control Conference., volume 3, pages 1483–1487, Albuquerque, New Mexico, June 1997. [68] J. Gigac. Using optical methods in the assessment of pulp and paper. Bul- letin, 32(2):23, 1996. [69] B. Gopaluni, M. S. Davies, P. D. Loewen, and G. A. Dumont. An online non-intrusive method for alignment between actuators and their response centers on a paper machine. ISA Transactions, Elsevier, 47:241–246, 2008. [70] B. Gopaluni, F. Farahmand, M. S. Davies, P. D. Loewen, and G. A. Dumont. An online non-intrusive method for alignment between actuators and their response centers on a paper machine. In Proceedings of the Control Systems Conference, pages 203–208, Tampere, Finland, June 2006. [71] D. M. Gorinevsky and C. Gheorghe. Identification tool for cross-driectional process. IEEE Transactions on Control Systems Technology, 11(5):629–640, September 2003. 141 Bibliography [72] D. M. Gorinevsky and E. M. Heaven. Automated identification of actuator mapping in Cross-Directional control of paper machine. Technical report, United States Patent PIN# 6,086,237, 18th October 1996. [73] D. M. Gorinevsky and E. M. Heaven. Automated identification of actuator mapping in cross-directional control of paper machine. In Proceedings of the American Control Conference, pages 3400–3404, Albuquerque, New Mexico, June 1997. [74] D. M. Gorinevsky and E. M. Heaven. Automated identification of web shrinkage and alignment parameters in sheet making machinery using a modeled actuator response profile. United States Patent US006086237A, Meaurex Devron Inc., July 2000. [75] D. M. Gorinevsky, E. M. Heaven, C. Hagart-Alexander, M. Kean, and S. Morgan. New algorithm for intelligent identification of paper alignment and nonlinear shrinkage. In Proceedings of the Control Systems Conference, pages 41–47, Halifax, Nova Scotia, 1996. [76] D. M. Gorinevsky, E. M. Heaven, C. Hagart-Alexander, M. Kean, and S. Morgan. New algorithm for intelligent identification of paper alignment and nonlinear shrinkage. Pulp and Paper Canada, 98(6):76–81, 1997. [77] D. M. Gorinevsky, E. M. Heaven, C. Lynch, and C. Hagart-Alexander. Au- tomated identification of web shrinkage and alignment in paper and film machines. In IEEE International Conference on Systems, Man and Cyber- netics, volume 4, pages 3358–3362, October 1995. [78] D. M. Gorinevsky, R. N. Vyse, and E. M. Heaven. Performance analy- sis of cross-direction process control using multivariable and spectral mod- els. IEEE Transactions on Control Systems Technology, 8(4):589–600, July 2000. [79] A. R. Guesalaga, F. M. Cabezas, and C. Gottschalk. Measuring physical properties of paper with an on-line visual sensor. Meas. Sci. Technol, 8:574– 580, 1997. 142 Bibliography [80] A. R. Guesalaga, A. D. Foessel, H. W. Kropholler, and F. Rodriguez. On- line measurement of paper shrinkage for monitoring & quality control. Pulp and Paper Canada, 95(7):36–40, 1994. [81] A. R. Guesalaga, F. Rodriguez, and H. W. Kropholler. Method for recog- nition of shrinkage using image processing. Celulosa y Papel (Celul. Pap. (Chile)), 8(3):10, 1992. [82] B. C. Haggman and W. L. Bialkowski. Performance of common feedback regulators for first-order and deadtime systems. In Control Systems, pages 75–83, Whistler, Canada, September 1992. [83] F. R. Hansen. A Fractional Representation Approach to Closed-Loop Sys- tem Identification and Experiment Design. PhD thesis, Stanford University, 1989. [84] F. R. Hansen, G. F. Frankin, and R. L. Kosut. Closed-loop identification via the fractional representation: experiment design. In Proceedings of the American Control Conference, pages 1422–1427, Pittsburgh, PA, June 1989. [85] R. A. Harshman. Foundations of the PARAFAC Procedure: model and conditions for an “explanatory” multi-mode factor analysis. UCLA Working Papers in Phonetics, 16(1):1–84, 1970. [86] R. A. Harshman. Determination and proof of minimum uniqueness condi- tions for parafac. UCLA Working Papers in Phonetics, 22:111–117, 1972. [87] B. Haznedar and Y. Arkun. Single and multiple property CD control of sheet forming processes via reduced order infinite horizon MPC algorithm. Journal of Process Control, 12(1):175–192, 2002. [88] G. X. He. Automatic cross-directional control zone alignment for sheetmak- ing systems. US Patent 5400258, Meaurex Corporation, Menlo Park, CA, USA, September 1995. [89] E. M. Heaven, D. M. Gorinevsky, M. Kean, and C. Sung. Integrated tool for intelligent identification of CD process alignment, shrinkage and dynamics. Pulp and Paper Canada, 99(2):40–44, 1998. 143 Bibliography [90] E. M. Heaven, D. M. Gorinevsky, C. Sung, and M. Kean. Multivariable identification of paper machine cross-direction process using an automated tool. In Dynamic modeling Control Applications for Industry Workshop; IEEE Industry Application Society, pages 56–61, May 1997. [91] E. M. Heaven, I. M. Jonsson, T. M. Kean, M. A. Manness, and R. N. Vyse. Recent advances in cross-machine profile control. IEEE Control Systems Magazine, pages 36–46, 1994. [92] S. M. Hoole, S. J. I’Anson, and R.W. Hoyland. CD shrinkage profiles on a commercial paper machine. Technical report, Paper Science Research Forum, 2000. [93] S. M. Hoole, S. J. I’Anson, M. Ora, T. N. Ashworth, D. Brigg, B. Phillips, and R.W. Hoyland. CD shrinkage profiles of paper: Experiments on a com- mercial paper machine. Paper Technology, 40(10):63–70, 1999. [94] M. Hovd, R. D. Braatz, and S. Skogestad. Optimality of SVD controllers. Technical report, Norwegian University of Science and Technology, Decem- ber 1996. [95] M. Hovd and S. Skogestad. Control of symmetrically interconnected plants. Automatica, 30(6):957–973, 1994. [96] M. A. Hubbe, J. J. Pawlak, and A. A. Koukoulas. Paper’s appearance: a review. NCSU BioResources Journal, 3(2):627–665, 2008. [97] R. S. Hunter and R. W. Harold. The measurement of appearance. John Wiley and Sons LTD, New York, 2nd edition, 1987. [98] S. J. I’Anson, R. P. A. Constantino, and S. M. Hoole. Estimation of the profile of cross-machine shrinkage of paper. Institute of Physics, 19(1):1– 11, 2008. [99] S. J. I’Anson and S. Kropholler. Enhancing visibility of wire mark by image analysis. Journal of Pulp & Paper Sciences, 17(1):J22–J26, 1991. 144 Bibliography [100] S. J. I’Anson and W. W. Sampson. The effect of raw material and paper machine variables on CD shrinkage profile of paper. Technical Report 41, Paper Sciences, 2001. [101] S. J. I’Anson and W. W. Sampson. Estimating paper machine CD shrinkage profiles from headbox actuator data. In Paper Technology, volume 45, pages 25–29, 2004. [102] T. Jiang and N. Sidiropoulos. Kruskal’s permutation lemma and the iden- tification of CANDECOMP/PARAFAC and bilinear models with constant modulus constraints. IEEE Transactions on Signal Processing, 52(9):2625– 2636, 2004. [103] A. P. Kaestner and C. M. Nilsson. Estimating the relative shrinkage pro- file of newsprint. Society of Photo-Optical Instrumentation Engineers, 42(2):1467–1475, May 2003. [104] S. Karimzadeh. Online detection of picketing and estimation of cross- direction and machine-direction variations using discrete cosine transform. Master’s thesis, University of British Columbia, Vancouver, BC, Canada, December 2008. [105] F. Kayihan. Practical estimation of CD control performance. In Proceed- ings of the American Control Conference, volume 5, page 3097, San Diego, California, June 1999. [106] H. A. L. Kiers and W. P. Krijnen. An efficient algorithm for PARFAC of three-way data with large numbers of observation units. Computational Psy- chometrics, 56(1):147–152, 1991. [107] A. P. Kjaer. Modelling, sensing and identification of web-forming processes. PhD thesis, Control Systems Centre, The University of Manchester Institute of Science and Technology (UMIST), 1994. [108] A. P. Kjaer, W. P. Heath, and P. E. Wellstead. Identification of cross- directional behaviour in web production: techniques and experience. Jour- nal of Control Engineering Practice, 3(1):21–29, 1995. 145 Bibliography [109] J. Kniivila. Optimization of machine direction tension profile using multi- ple CD actuator systems in a papermaking process. In Proceedings of the Control Systems Conference, Stockholm, Sweden, June 2002. [110] R. Kosut. Iterative Adaptive Control: Windsurfing with Confidence, chapter Model Identification and Adaptive Control: From Windsurfing to Telecom- munications. Springer Verlag, 2001. [111] K. Kristinsson and G.A. Dumont. Cross directional control on paper ma- chines using gram polynomials. Automatica, 32(4):533–548, 1996. [112] J. B. Kruskal. Three-way arrays: Rank and uniqueness of trilinear decom- positions, with application to arithmetic complexity and statistics. Linear Algebra and itsApplications, 18:95 – 138, 1977. [113] L. D. Lathauwer and J. Castaing. Tensor-based techniques for the blind separation of DS-CDMA signals. Signal Processing, Elsevier, 87:322–336, 2007. [114] D. Laughlin, M. Morari, and R. D. Braatz. Robust performance of cross-directional control systems for web forming processes. Automatica, 29(6):1395–1410, 1993. [115] W. S. Lee, B. D. O. Anderson, R. L. Kosut, and I. M. Y. Mareels. A New Approach to Adaptive Robust Control. International Journal of Adaptive Control and Signal Processing, 7:183–211, 1993. [116] W. S. Lee, B. D. O. Anderson, R. L. Kosut, and I. M. Y. Mareels. On Robust Performance Improvement through The Windsurfer Approach to Adaptive Robust Control. In Proceedings of the 32nd Conference on Decision and Control, pages 2821–2827, December 1993. [117] W. S. Lee, B. D. O. Anderson, I. M. Y. Mareels, and R. L. Kosut. On Some Practical Issues in System Identification for the Windsurfer Approach to Adaptive Robust Control. In 10th IFAC Symp. on System Identification, 1994. 146 Bibliography [118] W. S. Lee, B. D. O. Anderson, I. M. Y. Mareels, and R. L. Kosut. On some key issues in the windsurfer approach to adaptive robust control. Automat- ica, 31(11):1649–1636, 1995. [119] E. Lindem. Paper surface properties versus stress during drying. In Pro- ceedings of the Paper Physics Conference, pages 327–339, 1991. [120] X. Liu and N. Sidiropoulos. Cramer-Rao bounds for low-rank decomposi- tion of multidimentional arrays. IEEE Transactions on Signal Processing, 49(9):2074–2086, 2001. [121] L. Ljung. System Identification: Theory for the user. Prentice Hall, NJ, 1999. [122] T. Mast, K. Starr, P. Tran, and S. C. Chen. New Optimization of CD Control for Global and Localized Profile Performance. In TAPPI Spring Technical Conference & Trade Fair, 650 Ackerman Road Columbus OH 43202, 2003. ABB Inc. [123] E. C. Di Mauro, K. R. Wadhams, and P. E. Wellstead. On- and off-line measurement of paper shrinkage. In Proceedings of the IEE International Conference on Control, volume 2, pages 1229–1234, Stevenage, UK, 1994. [124] D. McFarlane and K. Glover. A loop shaping design procedure using H∞ synthesis. IEEE Transactions on Automatic Control, 37(6):759–769, 1992. [125] D. A. McFarlin. Control of cross-machine sheet properties on paper ma- chines. In Process Control Symposium, pages 49–54, 1983. [126] T. Metsälä and J. Shakespeare. Automatic identification of mapping and response for paper machine cross directional control. In Control Systems, pages 298–306, Porvoo, Finland, September 1998. [127] T. Metsälä and J. Shakespeare. Method and apparatus for identifying map- ping of paper machine actuator. US Patent US7128808B2, Metso Automa- tion Oy, Helsinki, Finland, October 31, 2006. 147 Bibliography [128] S. Mijanovic. Paper Machine Cross-Directional Control Near Spatial Do- main Boundaries. PhD thesis, University of British Columbia, March 2004. [129] D. Moore. ABB Control CD solutions with LV control perform very well. The resulting profile uniformity is consistent and stable. ABB Inc., 2008. [130] M. Morari and E. Zafiriou. Robust Process Control. Prentice Hall, Engle- wood Cliffs, NJ, 1989. [131] G. C. Myers. Effect of restraint before and during drying on edgewise com- pressive properties of handsheets. Journal of Tappi, 50:97–100, 1967. [132] A. J. Niemi, J. Berndtson, and S. Karine. Towards faster control of paper machine by dry line measurement. In International CD Symposium, pages 214–219, 1997. [133] C. M. Nilsson, L. Malmqvist, and J. Carlsson. Red fluorescence sensor for noncontact on-line measurements in paper production. Journal of Optical Engineering, 40:1674–1681, 2001. [134] L. S. Nordman. Laboratory investigations into the dimensional stability of paper. Journal of Tappi, 41:23–30, 1958. [135] S. Nuyan. Requirements for the cross-machine direction basis weight con- trol of the paper web. In Proceedings of the American Control Conference, pages 1416–1422, Seattle, WA, June 1986. [136] S. Nuyan and J. Shakespeare. Robustness and stability in CD control. In Control Systems, pages 193–196, Victoria, Canada, May 2000. [137] G. Obinata and B. D. O. Anderson. Model Reduction for Control System Design. Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2001. [138] N. Pauler. Paper optics. Lorentzen and Wettre Corp., Communications Dept., Kista, Sweden, 2002. [139] A. J. Paulraj and C. B. Papadias. Space-time processing for wireless com- munications. IEEE Signal Processing Magazine, 14:49–83, 1997. 148 Bibliography [140] M. Phan, J. Juang, L. Horta, and R. Longman. System identification from closed-loop data with known output feedback dynamics. Journal of Guid- ance, Control, and Dynamics, 17:661–669, 1994. [141] B. R. Philips, S. J. I’Anson, and S. M. Hoole. CD shrinkage profiles of paper-curve fitting and quantitative analysis. Appita Journal, 55(3):235– 243, 2002. [142] S. J. Popson and D. D. Malthouse. Measurement and control of the optical properties of paper. Technidyne Corp., 2nd edition, 1996. [143] H. Prasst. Differential shrinkage of paper in CD. In Proceedings of 2nd European Research Symposium on Image Analysis for Pulp and Paper Pro- duction, pages 95–111, Darmstadt, Germany, September 1993. [144] W. Reinelt. A ν-gap factsheet-with applications to model validation. Tech- nical report, Division of Automatic Control, Department of Electrical Engi- neering, Linköping University, 581 83 Linköping, Sweden, May 2001. [145] C. P. Robert and G. Casella. Monte carlo statistical methods. Springer texts in statistics. Springer-Verlag New York, Inc., 2nd edition, 1999. [146] O. J. Rojas, G. C. Goodwin, and G. V. Johnson. Spatial frequency an- tiwindup strategy for cross-directional control problems. IEE Proceeding Control Theory & Applications, 149(5):414–422, 2003. [147] R. T. Ross and S. Keugans. Component resolution using multilinear models. Methods in Enzymology, 246:679–700, 1995. [148] K. S. Rucker. Dryer for paper machine. US Patent 5729913, Meaurex Cor- poration, Pasco, Washington, August 1998. [149] W. Rudin. Fourier analysis on groups. New York: Interscience-Wiley, 1962. [150] E. L. Russell and R. D. Braatz. The average-case identifiability and control- lability of large scale systems. Journal of Process Control, 12(7):823–829, 2002. 149 Bibliography [151] D. R. Saffer and F. J. Doyle. Closed-loop identification with MPC for an industrial scale CD-control problem. IEE Proceeding Control Theory & Applications, 149(5):448–456, September 2002. [152] D. R. Saffer, F. J. Doyle, A. Rigopoulos, and P. Wisnewski. MPC study for a dual headbox CD control problem. Pulp and Paper Canada, 102(12):97– 101, 2000. [153] G. Schmidt. Optical properties of paper. Technical report, Verlag Dr. Martin Saendig GmbH., D-6229 Walluf 1/Wiesbaden, Germany, 1976. [154] W. E. Scott, J. C. Abbot, and S. Trosset. Properties of paper: An Introduc- tion. TAPPI Press, 1995. [155] W. E. Scott, L. R. Dearth, and B. D. Jordan. Optical properties of paper, volume 9, chapter 6, in Mill Control and Control Systems: Quality and Testing, Environmental, Corrosion, Electrical, pages 152–191. Pulp and Paper Manufacture, 1992. [156] J. Shakespeare, J. Pajunen, V. Nieminen, and T. Metsala. Robust optimal control of profiles using multiple CD actuator systems. Pulp and Paper Canada, 104(11):261–264, 2003. [157] J. A. Shands and J. M. Genco. Cross-Machine variation of paper curl on a twin-wire machine. Journal of Tappi, 7(9):165–169, 1988. [158] N. Sidiropoulos and R. Bro. On the uniqueness of multiliner decomposition of n-way arrays. Journal of Chemometrics, 14:229–239, 2000. [159] N. Sidiropoulos, R. Bro, and G. Giannakis. Parallel factor analysis in sensor array processing. IEEE Transactions on Signal Processing, 48(8):2377– 2388, 2000. [160] N. Sidiropoulos, G. Giannakis, and R. Bro. Blind PARAFAC receivers for DS-CDMA systems. IEEE Transactions on Signal Processing, 48(3):810– 823, 2000. 150 Bibliography [161] S. Skogestad, M. Morari, and J. C. Doyle. Robust control of ill-conditioned plants: high-purity distillation. IEEE Transactions on Automatic Control, 33(12):1092–1105, December 1988. [162] S. Skogestad and I. Postlethwaite. Multivariable feedback control: Analysis and design. New York: Wiley, 1996. [163] A. Smilde, R. Bro, and P. Geladi. Multi-way analysis. Applications in the chemical sciences. Wiley, Chichester, UK, 2004. [164] G. E. Stewart. Two Dimensional Loop Shaping Controller Design for Pa- per Machine Cross-Directional Processes. PhD thesis, University of British Columbia, August 2000. [165] G. E. Stewart. Reverse bumptest for closed-loop identification of CD con- troller alignment. US Patent US 2007/0039705 A1, Honeywell International Inc., 2007. [166] G. E. Stewart, J. U. Backstrom, P. Baker, C. Gheorghe, and R. N. Vyse. Controllability in cross-directional processes. Practical rules for analysis and design. Pulp and Paper Canada, 103(8), 2002. [167] G. E. Stewart, P. Baker, D. M. Gorinevsky, and G. A. Dumont. An ex- perimental demonstration of recent results for spatially distributed control systems. In Proceedings of the American Control Conference, pages 2216– 2221, Arlington, VA, USA, 2001. [168] G. E. Stewart, D. M. Gorinevsky, and G. A. Dumont. Feedback Controller Design for a Spatially-Distributed System: The Paper Machine Problem. IEEE Transactions on Control Systems Technology, 11(5):612–628, 2003. [169] G. E. Stewart, D. M. Gorinevsky, and G. A. Dumont. Two-dimensional loop shaping. Automatica, 39(5):779–792, 2003. [170] J. H. Stock and M. W. Watson. Introduction to Econometrics. Addison- Wesley, 2nd edition, 2007. 151 Bibliography [171] S. Talwar and M. Viberg. Blind separation of synchronous co-channel digital signals using an antenna array- Part I: algorithms. IEEE Transactions on Signal Processing, 44:1184–1197, 1996. [172] A. Taylor and S. Duncan. Actuator mapping and stability in real-life cross- directional control systems. In Proceedings of the Control Systems Confer- ence, pages 197–202, Tampere, Finland, June 2006. [173] A. Taylor and S. Duncan. Bayesian methods for system identification in cross-directional control. In Proceedings of the Control Systems Conference, pages 191–196, Tampere, Finland, June 2006. [174] J. F. M. ten Berge and N. Sidiropoulos. On uniqueness in CANDE- COMP/PARAFAC. Psychometrika, 67:399–409, 2002. [175] G. Tomasl. PARAFAC and missing values. Chemometrics and Intelligent Laboratory Systems, 75(2):163–180, February 2004. [176] P. Q. Tran, K. Starr, and T. Mast. Method and apparatus for controlling cross-machine direction (CD) controller settings to improve CD control per- formance in a web making machine. US Patent 7300548 B2, ABB Inc., Filed on Nov. 27, 2007. [177] A. J. van der Veen. A subspace approach for blind space-time signal pro- cessing. IEEE Transactions on Signal Processing, 45(1):173–190, 1997. [178] A. J. van der Veen. Algebraic methods for deterministic blind beam- forming. IEEE Proceeding, 86:1987–2008, 1998. [179] J. G. VanAntwerp, A. P. Featherstone, R. D. Braatz, and B. A. Ogunnaike. Cross-directional control of sheet and film processes. Automatica, 43:191– 211, 2007. [180] S. M. Veres and D. S. Wall. Synergy and Duality of Identification and Con- trol. Taylor and Francis Inc., 2000. [181] P. H. Viitaharju and K. J. Niskanen. Dried-in shrinkage profiles of paper web. Journal of Tappi, 76(8):129–134, 1993. 152 Bibliography [182] K. Villforth and S. Schabel. Clustermethoden für die analyse von sieb- mustern. Papier Technik Mensch, 4(10), 2005. [183] G. Vinnicombe. Frequency domain uncertainty and the graph topology. IEEE Transactions on Automatic Control, 38(9):1371–1383, 1993. [184] G. Vinnicombe. Measuring the robustness of feedback systems. PhD thesis, University of Cambridge, Cambridge CB2 1PZ, UK, 1993. [185] G. Vinnicombe. A ν-gap distance for uncertain and nonlinear systems. In Proceedings of the 38th Conference on Decision and Control, pages 2557– 2562, Phoenix, AZ, USA, 1999. [186] G. Vinnicombe. Uncertainty and feedback-H∞ loop-shaping and the ν-gap metric. Imperial Collage Press, 2001. [187] R. Vyse, B. Zarbi, U. Dos Santos, J. Rumel, and T. Steele. Trends to- wards finer cross direction moisture profile control. Pulp and Paper Canada, 102(11):T318–322, November 2001. [188] K. R. Wadhams, S. J. I’Anson, D. M. James, and H. W. Kropholler. The measurement of differential CD shrinkage. Paper Technology, 32(1):36–38, 1991. [189] T. Wahlström, K. Adolfsson, S. Östlund, and C. Fellers. Numerical mod- elling of the cross direction shrinkage profile in a dryer section, a first ap- proach. In TAPPI International Paper Physics Conference, pages 517–531, San Diego, 1999. [190] T. Wahlström and J. O. Lif. Dryer section simulator for laboratory investi- gations of shrinkage profile. In TAPPI International Paper Physics Confer- ence, pages 169–174, Victoria, Canada, 2003. [191] T. Wahlström and P. Makela. Predictions of anisotroic multiply board prop- erties based on isotropic ply properties and drying restraints. In Fundamen- tal Research Symposium, Cambridge, volume 13, pages 241–281, Septem- ber 2005. 153 Bibliography [192] G. L. Wedel. No-Draw drying restraint. In Tappi Engineering Conference, pages 275–281, 1988. [193] R. G. Jr. Wilhelm and M. Fjeld. Control algorithms for cross directional control: the state of the art. In Preprints of 5th IFAC/IMEKO Conference on Instrumentation and Automation in the Paper, Rubber, Plastics and Poly- merization Industries, pages 139–150, Antwerp, Belgium, 1983. [194] K. Zhou. Essentials of robust control. Prentice Hall, NJ, 1998. [195] K. Zhou, J. C. Doyle, and K. Glover. Robust and optimal control. Engle- wood Cliffs, NJ: Prentice-Hall, 1996. 154 Appendix A Alignment Programs This program generates the simulated data for a given dimension of the paper ma- chine. This routine first generates the closed loop data and then fits them in a tensor. Finally the generated tensor is decomposed to find the shrinkage model and the actuators’ response. The code can be found at the following address: http : //www.ece.ubc.ca/ ∼ fazelf/programs/parafac.zip 155 Appendix B Windsurfing Programs B.1 IMC-based Windsurfing This program is the main routine of the the Classical Iterative Windsurfer Adaptive Robust Control that simulates a paper machine and does system identification and control design. The code can be found at the following address: http : //www.ece.ubc.ca/ ∼ fazelf/programs/imcwindsurfing.zip B.2 H∞ Windsurfing This program applies the H∞-based Windsurfing on a simulated paper machine. The code is available at: http : //www.ece.ubc.ca/ ∼ fazelf/programs/hinfwindsurfing.zip B.3 ν-gap This program applies the enhanced Windsurfing, which is based on ν-gap to a simulated CD process. The code is available at: http : //www.ece.ubc.ca/ ∼ fazelf/programs/nugapwindsurfing.zip 156

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0070915/manifest

Comment

Related Items