MODEL PREDICTIVE CONTROL FOR MULTIPLE CROSS-DIRECTIONAL PROCESSES: ANALYSIS, T U N I N G , AND IMPLEMENTATION By Junqiang Fan B . Eng. (Electrical Engineering), Zhejiang University, 1988 M . Eng. (Electronics Engineering), Shenyang University of Technology, 1991 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF DOCTOR OF PHILOSOPHY " in T H E F A C U L T Y OF G R A D U A T E STUDIES DEPARTMENT OF ELECTRICAL A N D COMPUTER ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH C O L U M B I A September 2003 © Junqiang Fan,, 2003 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Electrical and Computer Engineering The University of British Columbia 2356 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract In this thesis, practical techniques for analyzing, tuning and implementing an industrial model predictive controller (MPC) for paper machine cross-directional (CD) processes are developed. These techniques include model reduction, performance and robust stability analysis of a linear closed-loop system, and optimal prediction of steady-state performance for constrained M P C . Paper machine cross-directional processes are large-scale two-dimensional systems. In order to make the online computational load feasible in real time, it is necessary to reduce the dimensions of the model, especially for multiple-array systems. Due to the very large dimension and ill-conditioning of the process, the plant model can be effectively reduced by wavelet matrices based on a modified wavelet packet method. The reduced model captures the controllable spatial components of CD processes. After model reduction, M P C based on the reduced model is implemented in the wavelet domain. Meanwhile, the constraints can be exactly represented in the wavelet domain. The main benefit of implementing M P C in the reduced domain is that the on-line computational time for solving the quadratic programming (QP) problem is greatly reduced compared to the implementation of M P C based on the original model while good control performance is preserved. This method is applicable for multiple-array systems. For large-scale spatially-distributed systems such as CD processes, the process model, the additive structured uncertainty, and the linear portion of the M P C controller are approximated as linear, spatially-invariant, and time-invariant. The transfer function analysis will hold as long as the disturbances are small enough magnitude that the M P C does not hit constraints. To analyze the relevant closed-loop transfer matrices, the novel concept of rectangular circulant matrices (RCMs) is proposed. RCMs can be diagonalized by complex Fourier matrices, allowing analysis in terms of a family of single-input singleoutput (SISO) transfer functions across the spatial frequencies. Familiar concepts from control engineering such as bandwidth and stability margin are extended into the twodimensional frequency domain, providing intuitive measures of closed-loop performance and robustness. Consistent criteria are given for the analysis of the closed-loop effect of the industrial C D M P C tuning weights based on standard robust control theory. Properly tuning the magnitude of weights, choosing the structure of weights, and considering the structure of model uncertainty in the M P C design stage can greatly improve the performance. This method can be used to design robust C D M P C for multiple-array systems. In order to assess the steady-state performance of constrained CD M P C , the state of the art requires to run closed-loop simulations. Due to the large-scale characteristic of C D processes, especially for multiple-array systems, it is very time-consuming and ii inconvenient for use as a practical tuning tool. However, fast and correct prediction of the steady-state performance is necessary for tuning the C D M P C , especially for the cases with active constraints. A one-step static optimizer is proposed to predict the steady-state closed-loop performance of CD M P C . Two examples are given for demonstrating that the static optimizer is significantly more efficient (up to two orders of magnitude) than the conventional closed-loop simulation method while'reliably and accurately predicting the steady-state closed-loop performance. iii Table of Contents Abstract ii List of Tables viii List of Figures ix Acknowledgements xv Chapter 1 1 Introduction 1.1 Paper Machine Cross-Directional Processes 1 1.2 Process Characteristics 6 1.2.1 Single-array System 6 1.2.2 Multiple-Array System 7 1.2.3 Ill-Conditioned System 8 1.2.4 Model Uncertainty 9 1.2.5 Coordination of Actuator Arrays 1.3 Industrial C D Model Predictive Control 12 1.4 Literature Overview 12 1.4.1 Spatial Frequency and Rectangular Plant 12 1.4.2 C D Model Predictive Control (MPC) 14 1.4.3 Robust C D Control 15 1.4.4 Steady-State Closed-Loop Performance with Constrained CD-MPC 1.5 Chapter 2 11 17 Aims and Contributions of the Work 17 1.5.1 Motivation and Objective 17 1.5.2 Contributions 18 1.5.3 Thesis Overview . : 18 Problem Statement 20 2.1 20 Process Model 2.1.1 Single-Array Process Model 2.1.2. Multiple-Array Process Model 20 24 2.2 Industrial Cross-Directional M P C 25 2.3 Model Reduction 27 2.4 Robust Stability and Performance 27 2.5 Steady-State Performance Prediction 28 iv Chapter 3 Two-Dimensional Frequencies in C D Control 29 3.1 3.2 3.3 3.4 Closed-Loop Transfer Functions 29 Rectangular Circulant Matrices 32 Frequency Domain Interpretation 34 Approximate Rectangular Circulant Matrices 38 3.4.1 Spatial Boundary Effects 39 3.4.2 Input and Output Dimensions Related by a Non-Integer Factor 40 3.5 Two-Dimensional Frequency Analysis 42 3.5.1 Open-Loop Transfer Functions in the Two-Dimensional Frequency Domain 42 3.5.2 Structured Uncertainty in the Two-Dimensional Frequency Domain 44 3.5.3 Closed-Loop Transfer Functions in the Two-Dimensional Frequency Domain 47 3.6 Summary 48 Chapter 4 M o d e l Reduction 50 4.1 Signal Decomposition and Reconstruction Using Discrete Wavelet Transform 50 4.2 Signal Decomposition and Reconstruction Using the Wavelet Matrices 53 4.3 Model Reduction 55 4.3.1 Modified Wavelet Packet Analysis 55 4.3.2 Model Reduction Using the Modified Wavelet Packets Method [21] 57 4.3.3 Quantification of the Model Reduction Requirements . 59 4.4 Practical Issues in Reduction Implementation 61 4.4.1 Distortion Correction 61 .4.4.2 Choice of Wavelets 62 4.4.3 Spatial Controllability and Model Reduction 64 4.4.4 Model Reduction Distortion Test 65 4.4.5 Comparison With Other Reduction Methods 66 4.5 MPC Implementation in the Wavelet Domain 68 4.5.1 Constraint Handling 68 4.5.2 Objective Function in the Wavelet Domain 69 4.6 Extension to Multiple-Array Systems 70 4.7 Examples 72 v 4.8 4.7.1 Basis Weight Profile Controlled by Slice L i p Actuators . 72 4.7.2 Three Properties Controlled by Four Actuator Arrays . 75 Summary 78 Chapter 5 Robustness and Performance 80 5.1 Closed-Loop Robust Stability for Single-Array Systems 5.1.1 80 Robust Stability Condition for the Unstructured Model Uncertainty 5.1.2 80 Robust Stability Condition for the Structured Model Uncertainty 81 5.2 Closed-Loop Performance Indices 82 5.3 Role of M P C Tuning Parameters 83 5.3.1 Structure of the Weighting Matrix Q 84 3 5.4 Robust Stability for Multiple-Array Systems 88 5.5 Examples 91 5.5.1 Unstructured Model Uncertainty Case 92 5.5.2 Structured Uncertainty vs Unstructured Uncertainty . . 94 5.5.3 Moisture Controlled by Steam B o x and Rewet Shower Actuators 5.6 96 Summary 99 Chapter 6 Steady-State Performance Prediction .6,1 100 C D - M P C Implementation for Multiple-Array Systems 100 6.1.1 Objective Function •. 100 6.1.2 State-Space Model 101 6.1.3 Output Predictions 102 6.1.4 Transfer to Standard Q P Problem 104 6.2 Static Optimizer 105 6.3 T w o Sensitivity Functions 107 6.4 6.5 6.6 6.3.1 Closed-loop Sensitivity Function for Unconstrained M P C 107 6.3.2 Sensitivity Function for the Static Optimizer 109 Design for Matching Optimization Problems 110 6.4.1 Approximate Solution for the Static Optimizer 110 6.4.2 Special Cases 112 Examples 113 6.5.1 Two Actuator Arrays and One Controlled Property 6.5.2 Four Actuator Arrays and Three Controlled Properties Summary . . 113 115 118 vi Chapter 7 Conclusions and Future Work 121 7.1 Summary 121 7.2 Future Work 123 Appendix A Closed-Loop Transfer Functions of Unconstrained M P C 124 Appendix B Wavelet Packet Analysis Tree 128 Appendix C Proof of R C M ' s Properties 130 C.1 Proof of R C M Property 1 130 C.2 Proof of R C M Property 3 130 C.3 Proof of R C M Property 4 130 C.4 Proof of R C M Property 2 134 C.5 Proof of RCM's Property 5 134 C. 6 Proof of RCM's Property 6 135 Appendix D Proof of Theorems 136 D. 1 Proof of Theorem 1 136 D.2 Proof of Theorem 2 137 Appendix E Calculation of the Eigenvalues of T u c iA 141 Appendix F Derivation of the Standard Q P Format 143 Appendix G Derivation of Equation (6.61) 144 Appendix H Relationship Between Appendix I G sum Solution of (6.95) and Gss 145 146 Bibliography 148 vii List o f Tables 4.1 The model reduction performance criteria for different wavelet filters ... 63 4.2 The model reduction performance criteria for the single-array system. . . . 75 4.3 Performance and computational time for the M P C using the full and reduced models 78 5.1 Robustness and performance: Q3 using different structures 88 5.2 The performance criteria for the conservative and aggressive controllers . . 92 5.3 RS and performance test results for structured and unstructured uncertainties (the single-array system) 5.4 96 RS and performance test results for structured and unstructured uncertainties (the multiple-array system) 6.1 99 Performance and computational time for the M P C using the closed-loop simulation and the static optimizer viii 120 List of Figures 1.1 Wide view of the paper machine showing typical positions of the various actuator arrays and scanning sensor(s). (Artwork courtesy of Honeywell) . . . . 1.2 Illustrating the scanning sensor's (scanner) zig-zag measuring path 1.3 (a) Steady-state response of basis weight to slice lip actuators; (b) Steady-state 2 3 response of moisture to slice lip actuators; (c) The slice lip actuator profile. The thick smooth solid lines in (a) and (b) are the modelled responses. Data was obtained from a Canadian newsprint machine 1.4 Two-dimensional feature of a C D process 1.5 Typical steady-state responses of basis weight to a slice lip actuator and a 4 7 dilution actuator for a A5gsm newsprint papermaking process 1.6 Singular values of a typical CD process model 1.7 The block diagram of Honeywell's industrial cross-directional model predictive controller 2.1 8 9 13 (a) Columns of the spatial response matrix Go in (2.4) identified by a basis weight process controlled by n = 64 slice lip actuators; (b) Corresponding slice lip actuators positions; (c) The singular values of Go (Data courtesy of Honey well) 21 2.2 The spatial interaction matrix G in (2.2) has a band diagonal structure. . . . 23 3.1 The closed-loop system block diagram with the unconstrained M P C in (2.26). 30 3.2 Equivalent transfer matrix representation of Figure 3.1 31 3.3 (a) Nonzero elements of a S C M N ; (b) Non-zeros of the diagonal elements 0 s of N ; (c) Nonzero elements of an R C M N ; (d) Non-zeros of the diagonal s r elements of N 35 r 3.4 Response of the closed-loop system to steady-state disturbances at spatial frequencies v = 1.42 cycles/metre and v = 2.84 cycles/metre, illustration of modelled process outputs \y\ (dotted) in (3.33) and approximated by \y\ in (3.35), and steady-state actuator set-points \u\ (dotted) in (3.34) and approximated by |w| in (3.36) 3.5 41 The two-dimensional frequency approximation giy?, z) in (3.42) of the nominal plant model G(z) in (2.1) 43 ix 3.6 The nominal spatial response g in (2.4) (solid+dotted line) and its correspondt ing one with uncertainty g ip in (2.10) (solid line): (a) the gain parameter uncertainty <5 = —0.25 in (2.13); (b) the attenuation parameter uncertainty 7 S = —0.25 in (2.13); (c) the width parameter uncertainty 8c = —0.25 in a (2.13); (d) the divergence parameter uncertainty Sp — 0.25 in (2.13); (e) the alignment parameter uncertainty 8 = —3 in (2.13); (f) all the parameters C uncertainty limits Si e [-0.25,0.25] and S <E [-0.6,0.6] in (2.13), (note 2 that only some extreme cases are shown for clear display) 3.7 The additive structured uncertainty \S(v™, e ™ in (3.44) for UJ \ l2 45 — 0 in the spa- tial frequency domain due to: (a) the gain parameter uncertainty S € 1 [—0.25,0.25] in (2.13); (b) the attenuation parameter uncertainty 5 a € [—0.25,0] in (2.13); (c) the width parameter uncertainty 8$ € [—0.25,0] in (2.13); (d) the divergence parameter uncertainty 8/3 € [0,0.25] in (2.13); (e) the alignment parameter uncertainty 8 6 [—3,3] in (2.13); (f) all the paC . rameters uncertainty limits 81 € [—0.25,0.25] and 8 <E [—0.6,0.6] in (2.13). 2 (note that the vertical axes have different scales for clear display) 46 3.8 One structured model uncertainty in the two-dimensional frequency domain. . 3.9 The structured model uncertainty in Figure 3.8 in the input singular vectors and dynamical frequency domain 47 48 3.10 The two-dimensional frequency approximation t d(Vj,z) in (3.29) of the sensi- y tivity function T d(z) in (3.6) with the nominal process model G(z) in (2.1) y and the controller parameters H = 8, H = l,Qi = Al, Q = 10~ /, and 4 p u 2 Q = 0.002/ in (2.26) 49 z 3.11 The two-dimensional frequency approximation t„d(^", z) in (3.30) of the closedloop transfer matrix T^z) in (3.4) with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H p = 1, Qi = Al, Q — 10~ /, 4 U 2 and Q = 0.002/ in (2.26) 49 3 4.1 (a)The scaling function <j>(t) of sym4; (b)The wavelet function ip(t) of sym4- • 4.2 Two-stage two-band analysis tree 52 4.3 Two-stage two-band synthesis tree 52 4.4 The decomposition and reconstruction filters of sym4 53 4.5 Spectra of the decomposition and reconstruction filters for syn%4 54 4.6 Two-stage wavelet packet analysis tree 56 4.7 A n example of the modified wavelet packet analysis 58 4.8 A n example for illustrating the singular values of 4.9 Illustration of the construction of the signals that define the model reduction, y in (2.1) and y in (4.18) WyW y (a) and W j W u (b). . 59 60 x 51 4.10 (a) The spectra of the wavelet filters: sym4 (dotted) and sym8 (dashed); (b) Illustration of the distortion difference in the reduced model G i n (4.25) r due to different wavelet filter lengths 64 4.11 (a) The spectra ofthe wavelet filters: coif2 (dotted) and sym6 (dashed); (b) Illustration of the distortion difference in the reduced model G in (4.25) r due to different wavelets 65 4.12 Singular values and distortions in a reduced model from dilution actuator response 66 4.13 (a) T h e probability of the distortion ratio rji in (4.35) using sym4; (b)The probability of the distortion ratio r\\ in (4.35) using sym8 67 4.14 Spectra of the approximation output profiles using different methods 68 4.15 M P C implementation in the wavelet domain 70 4.16 Basis weight response to a slice lip actuator and its power spectrum 73 4.17 The singular values and distortion indices of the reduced model 74 4.18 (a) The output profile before control, (b) The output profile after control using the reduced model, (c) Difference of two controlled output profiles using the original model and the reduced model, (d) Spatial power spectra of the output profiles (dotted: before control, solid plus 'circle': after control using the original model, solid: after control using the reduced model). Note that the x-axis of (d) is enlarged for clear display 76 4.19 Controller performance (standard deviation) and computational time for different reduced models 77 4.20 (a) Performance and (b) computational time for the multiple-array system M P C using the full model and the reduced model 5.1 79 A n example shows that the controller satisfies the R S condition (5.6) for structured uncertainties but violates the R S condition (5.1) for unstructured uncertainties 5.2 82 Contour plot illustrating the closed-loop bandwidth contour in the two-dimensional frequency domain from Definition 2 for the closed-loop transfer function t d(Vj,z) y in (3.29) with the nominal process model G(z) i n (2.1) and the controller parameters H — 8, H p = l , Q i = 41,Q = 1 0 / , and Q = -4 u 2 3 0.002/ in (2.26) (compare with the surface plot in Figure 3.10) 5.3 83 Illustration of the effect of changing Q on the gain of the closed-loop transfer 2 function t d(Vj,z) y i n (3.29) at spatial frequency v]j = 0 with respect to dynamical frequency with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H = 1, Qi — 41, and Q — 0.002/ i n (2.26). p U 3 xi 85 5.4 Contour plot illustrating the effect of changing Q in (2.26) on the closed2 loop bandwidth contour B (Uj,u) w in the two-dimensional frequency domain for the closed-loop transfer function i «f (z/V, z) in (3.29) with the nominal y process model G(z) in (2.1) and the controller parameters H — 8, H — p u l,Qi = AI, and Q = 0.0027 in (2.26). Note the greater influence on the 3 dynamical component of B (Vj,u) 85 w 5.5 (cf. Fig. 1 in [20]) Illustration of the effect of changing Q in (2.26) on the 3 gain of the closed-loop transfer function t i(Vj,z) in (3.29) at steady-state yc with respect to spatial frequency with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H p = l,Qi = 47, Q = 1 0 " / 4 U 2 in (2.26) 5.6 86 (cf. Fig. 2 in [20]) Illustration of the effect of changing Q in (2.26) on the 3 gain of the closed-loop transfer function t d(y^, z) in (3.30) at steady-state u with respect to spatial frequency with the nominal process model G{z) in (2.1) and the controller parameters H = 8, H = \,Q\ = 47,Q = 10~ / 4 p u 2 in (2.26) 5.7 86 Contour plot illustrating the effect of changing Q in (2.26) on the closed3 loop bandwidth contour B (Vj,u>) in the two-dimensional frequency domain w for the closed-loop transfer function t {Vj,z) in (3.29) with the nominal yd process model G(z) in (2.1) and the controller parameters H — 8, H p u = l,Qi = 47, and Q — 10~ / in (2.26). Note the greater influence on the 4 2 spatial component of B (U^,OJ) as compared with Figure 5.4 w 5.8 The spatial frequency representation |6 (f")| of the bending moment matrix m Bm m 5.9 87 (2.31) . 88 The sensitivity functions \t d(Vj, 1)| in the spatial frequency domain: Diagonal y Q vs Q using B 3 3 5.10 A n example of M T u d A in (5.10) 89 which shows the nonzero elements. 91 5.11 (a) and (b): the disturbance profile d(z) in (3.6) before control and,its spatial spectrum \d(z)\ = \F • d(z)\; (c) and (d): the basis weight profile y(z) in m (3.6) and its spatial spectrum \y(z) \ = \F -y(z) \ after using the conservative m controller; (e) and (f): the basis weight profile y(z) in (3.6) and its spatial spectrum \y(z) \ = \F -y(z) \ after using the aggressive controller (note that m (b), (d), (f) have different vertical scales) xii 93 5.12 (a) and (c): Illustration of t {v^(? ™) 2 ud in (3.28) and t {^,e ™) i2 yd in (3.27) in the two-dimensional frequency domain which satisfies the RS condition (5.1) for the unstructured uncertainty case; (b) and (d): Illustration of z) in (3.28) and t d{Vj, z) in (3.27) in the two-dimensional frequency tud(Vj, y domain which satisfies the RS condition (5.6) for the structured uncertainty case 95 5.13 Application of the C D - M P C controller designed for the structured uncertainties to constrained cases: (a). The output profile before control; (b) a class of output profiles after control; (c) a class of variations of the output profiles with time; (d) a class of input profiles after control; (e) the second order derivatives of the input profiles after control; (f) a class of variations of the input profiles with time (note that the y-axes are different in (a) and (b) for clear display 97 5.14 (a) and (c): Illustration of the singular values of T {z) in (5.15) in the dyud namical frequency and singular vectors domain and the sensitivity function t d{v j,z) in (3.27) in the two-dimensional frequency domain for the unV V structured uncertainty case; (b) and (d): Illustration of the singular values of T (z) in (5.15) in the dynamical frequency and singular vectors domain ud and the sensitivity function t d(Vj,z) y in (3.27) in the two-dimensional fre- quency domain for the structured uncertainty case (for the multiple-array system) 6.1 98 The block diagram of an unconstrained M P C diagram for a multiple-array C D system 6.2 108 (a) The moisture profiles before control (dashed) and after control (solid: static optimizer; dotted: closed-loop simulation); (b) The steam box actuator profiles (solid: static optimizer; dotted: closed-loop simulation); (c) The rewet shower profiles (solid: static optimizer; dotted: closed-loop simulation). 115 6.3 The basis weight profiles for the example in Section 6.5.2: static optimizer (solid+asterisk), closed-loop simulation (solid+circle), disturbance (dashed line), set-points (flat solid line) 6.4 117 The moisture profiles for the example in Section 6.5.2: static optimizer (solid+asterisk), closed-loop simulation (solid+circle), disturbance (dashed line), set-points (flat solid line) 6.5 118 The caliper profiles for the example in Section 6.5.2: static optimizer (solid+asterisk), closed-loop simulation (solid+circle), disturbance (dashed line), set-points (flat solid line) 119 xiii The actuator profiles for the example in Section 6.5.2 using the static optimizer (dashed plus point line) and closed-loop simulations (solid line) xiv 119 Acknowledgements I wish to thank my research supervisor Prof. G u y A . Dumont at U B C for his guidance and support during my postgraduate studies. I would also like to thank my industrial supervisor D r . Gregory E . Stewart at Honeywell Industry Solutions for his advice and friendship throughout the course of this project. The efforts of my committee members and readers, Drs. Stephen Duncan, Michael Davies, Peter Lawrence, V i k r a m Krishnamurthy, John Meech, W i l l i a m G . Dunford, E z r a Kwok, and Nicolas Jaeger are grateful acknowledged. I would like to thank my past and present members of control group at the P u l p and Paper Centre, Stevo Mijanovic, Shiro Ogawa, Dr. Michael Chong Ping, Manpreet Sidhu, Stephan Bibian, Setareh Aslani, D r . Leonardo Kammer, Dr. Mihai Huzmezan for their valuable discussions and friendship. The help of the Pulp and Paper Centre staff is greatly appreciated. In addition, Drs. Bruce Allison, Jan Maciejowski, and W i l l Heath have all provided technical input that shaped this work. I also would like to thank Honeywell Industry Solutions - Vancouver Operations for technical and financial support of the project. I am grateful to Johan Backstrom, Pengling He, Bob Vyse, Rob Spinner, Dan Stevens, Scott Morgan, Cristian Gheorghe, P a u l Baker, Amor Lahouaoula for their help. The financial support of the British Columbia Advanced Systems Institute, and the Natural Sciences and Engineering Research Council of Canada is greatly acknowledged. I would like to thank my mother, Caigui X u , my brothers and sister for their continuous support during the past 30 years. Finally, I would pay special tribute to my wife, Hongying (Mary) Zhang, and my daughter, Huiqiao (Melinda) Fan for their love and support during my studies at U B C . xv In Memory of My Father: Biquan Fan (1931-1983) I would like to dedicate this thesis to my father. Without his love and encouragement during my early education years, I would never have fulfilled my dream - to become an electrical engineer. xvi Chapter 1 Introduction 1.1 Paper Machine Cross-Directional Processes Papermaking is an important industry worldwide worth billions of dollars annually. The paper machine plays one of the most important roles in paper production. Its function is to continuously transform a slurry of water and cellulose fibers into a sheet of paper as efficiently as possible. Figure 1.1 illustrates a typical Fourdrinier style paper machine. The following is a brief description of its operation and more detailed information may be found in [9, 14, 31]. In general, a paper machine can be divided into four sections: the wet end section, the press section, the dryer section, and the post drying operations section. First, a dilute suspension of fibers, typically 0.3 ~ 0.6% consistency (i.e., 0.3 ~ 0.6% fibers, and 1 99.4 ~ 99.7% water), is pumped into the headbox. The headbox then evenly distributes the slurry over the cross direction on an endless wire screen which may be moving at speeds up to 120 km/h. Water in the slurry is removed by gravity and suction devices under the wire screen. The sheet of paper is formed and at this point (the end of the wet end section) it is 18 ~ 23% consistency. More water is squeezed out in the press section by counterrotating rollers to a consistency of 35 ~ 55%. Then, the sheet enters the dryer section and is dried by a series of steam heated dryer cans to 91 ~ 95% dryness (finished paper contains 5 ~ 9% water content). The dry sheet then enters the post drying operations section and passes a series of rotating rolls known as the calender stack for improving the surface smoothness and gloss, and controlling the sheet thickness (caliper). Finally, the dry paper is wound up on a reel as shown in Figure 1.1. The direction of sheet travel is known as the machine direction ( M D ) and the direction perpendicular to the sheet travel is known as the cross direction ( C D ) . In industry, paper machine control is split in two parts: M D control and C D control. The paper properties are measured by sensors mounted on a scanner that traverses back and forth across the moving paper sheet, taking between 15 — 45seconds to cross its width. In each pass across the moving paper sheet, the scanning sensor measures the sheet properties at up to 2000 evenly spaced locations, called measurement points or data boxes. The paper properties of particular concern to the papermaker are moisture (water 1 Consistency is defined as weight concentration of fibers in the suspension 1 Chapter 1. slice lip Introduction 2 machine direction steam box ' dilution actuator wet end section r e w e t showe press section cross direction dryer section scanning sensor induction heating post drying operations section Figure 1.1: Wide view of the paper machine showing typical positions of the various actuator arrays and scanning sensor(s). (Artwork courtesy of Honeywell) content of the sheet weight in %), basis weight or grammage (the sheet weight per unit area in grams per square metre gsm), caliper (the thickness of the sheet in //m), and gloss (a measure of the optical reflectance of the sheet measured in points). As illustrated in Figure 1.2, the combined effect of the moving paper sheet with the motion of the scanning sensor results in the sensor measuring a zig-zag path over the produced sheet. However the industrial practice and the actuator design are concerned with the sheet properties partitioned into pure MD component and pure C D profile, while the extraction of true C D and MD signals from scanning sensor measurements is currently unresolved (see the research on two-dimensional aliasing issues in [73, 91]). The industrial practice is to simply filter the scan profiles and create estimates of the CD and MD signals. Such filtering will generally slow down the observed dynamics of the system in an effort to produce a more reliable estimation. In order to control the time-varying M D component, the overall paper weight is controlled by manipulating the flow rate of pulp into the paper machine and the speed of the machine itself. The M D moisture content is controlled via the overall drying effect of the machine, typically actuated by manipulating the overall steam flow into the dryer section. M D control is outside of the scope of the current work, and no M D actuators are illustrated in Figure 1.1. Cross-directional control refers to the control of the paper profile in the CD. This is a finer-resolution control that is actuated by arrays of identically-constructed actuators distributed across the CD width of the paper machine. Figure 1.1 illustrates five of the most common types of CD actuator arrays in their typical positions on the paper machine Chapter 1. Introduction 3 Figure 1.2: Illustrating the scanning sensor's (scanner) zig-zag measuring path. including slice lip, dilution actuators, steam boxes, rewet shower and induction heating actuators. Slice lip actuators are traditionally designed for controlling C D variations in the basis weight profile. The pulp exits the headbox shown in Figure 1.1 through a narrow opening across the width of paper machine. The size of this opening can be changed locally by bending the flexible slice lip using a set of actuators (screws). Tightening one of the screws locally narrows the opening thus decreasing the amount of fibers at that location on to the wire and decreasing the basis weight. It is worth noting that the change in the basis weight profile by the slice lip will also bring about changes in the moisture and caliper profiles [3, 14]. Increasing the basis weight at one location in the C D will reduce drainage and thus locally increase moisture, and locally increase the thickness of the sheet due to increased fibre deposition in the absence of closed-loop moisture and caliper control. Figure 1.3 illustrates the steady-steady responses of basis weight and moisture to slice lip actuators with data from a Canadian newsprint paper machine. Figures 1.3a and 1.3b show the basis weight and moisture profiles including the process responses and the modelled responses respectively, and Figure 1.3c shows the actuator profile. Note that all of these figures are showing the deviations attributed to the bump test, which is normally used for identifying the process model. A typical set of slice lip actuators (screws) has n = 50 ~ 118 actuators and the actuator spacing is between 7.0cm and 20cm [14, 77]. Besides the common maximum/minimum and rate constraints, slice lip actuators are subject to bending constraints in order to prevent damaging the expensive slice lip due to excessive flexing and bending. Chapter 1. Introduction 4 200 400 600 800 1000 1200 1000 1200 CD locations [measurement points] 400 600 800 CD locations [measurement points] 200 (c) 100 I 0 -100 -200 10 15 20 ' 25 30 I 35 40 ' 45 50 55 Actuator numbers Figure 1.3: (a) Steady-state response of basis weight to slice lip actuators; (b) Steady-state response of moisture to slice lip actuators; (c) The slice lip actuator profile. The thick smooth solid lines in (a) and (b) are the modelled responses. Data was obtained from a Canadian newsprint machine. A newer type of actuator for basis weight control is the dilution actuator [90]. The basis weight profile can be changed by an array of dilution actuators distributed across the back of the headbox through locally altering the concentration of pulp fibers in the headbox. Increasing the flow of a dilution actuator reduces the local concentration of pulp fibers and thus decreases the local basis weight. A typical set of dilution actuators has n = 150 ~ 300 actuators with actuator spacing 3.5cm ~ 7.0cm. Detailed information about basis weight response to a dilution actuator may be found in [77, 90]. Again the adjustment in the basis weight profile by the dilution actuators will likely impact the moisture and caliper profiles although experimental data do not yet exist. Steam boxes utilize the principle of hot pressing to improve water removal by increasing the sheet temperature. The sheet temperature is increased by condensing steam on and in the sheet and thus the water in the sheet can be more easily removed. The number of actuators in an array of steam boxes is n = 55 ~ 171 and actuator spacing is between 7.5cm ~ 15cm. Analogous to slice lip actuators, steam boxes have an impact on both moisture and caliper [14]. Chapter 1. Introduction 5 Rewet shower actuators are used for preventing over drying of the sheet and correcting the dry steaks. These actuators atomize water into fine droplets and spray on the sheet in order to reduce the variations of the moisture profile. There are n = 50 ~ 120 actuators for a rewet shower actuator array with actuator spacing 5cm ~ 15cm. The calendar stack is located at the dry end of paper machine and is a series (stack) of counter-rotating rollers through which the paper sheet is fed in order to modify the caliper or thickness of the sheet. C D control of the caliper profile is actuated via an array of induction heating actuators that are mounted about 5mm from one of the rollers. Each actuator in the array is an independently controlled A C current source. Increasing the amplitude of the current has the effect of inducing heat-generating eddy currents in the roller. As the temperature increases, the diameter of the roller locally increases at that CD location, due to thermal expansion. This in turn increases the applied pressure on the paper sheet and reduces the resulting caliper. Induction heating actuators thus have a strong impact on the caliper profile and have also been observed to positively influence the sheet gloss due to the effect of increased surface temperature and applied pressure. A typical set of induction heating actuators has n = 100 actuators on average with actuator spacing 7.5cm. As seen above, although each actuator array is originally designed for controlling a single sheet property, strong interactions exist between multiple actuator arrays and multiple sheet properties. The estimated standard deviation of the paper properties is used as a quality standard in the sale of paper [1]. If a reel of paper is measured by the scanning sensor in Figure 1.1 at m locations in the C D and the reel consists of T scans, then the sheet measurements are given by y{k) e R m with discrete-time index k = 1, • • • , T, and the conventional paper quality index K is defined as (1.1) T m vA ) k m-T where y avg (1.2) is the overall mean of the two-dimensional array of measurements. Minimizing the the quality index K has been used in many C D control strategies. Chapter 1. Introduction 1.2 6 Process Characteristics The goal of this section is to describe the features of the C D process that make its control a challenging problem. 1;2.1 Single-array System In most C D control literature, a C D process only describes a single actuator array and single controlled property case, which is modelled as a cascade connection of a dynamic and a spatial response. In other words, the discrete-time model of the process has separable dynamic and spatial responses given by, y(z)=G(zHz), where y(z) € C ,u(z) m G{z) = G -h(z), 0 (1.3) e C " are the Z-transforms of the measurement profile and the actuator set-point profile. As discussed in Section 1.1, the scanning sensor measures the sheet at m < 2000 locations in the CD. The actuator array is composed of n < 300 individual identically-constructed actuators. Go € ] R m x n is the constant interaction matrix which describes the spatial response of the process, h(z) is the Z-transform of the scalar dynamic response of the process, which is usually considered as a first-order-plus-deadtime transfer function. The separability of the dynamics h(z) from the spatial response matrix G in (1.3) is a common simplification in C D control. Its implication is that the response 0 of the paper sheet to one of the actuators does not change in shape after changing the actuator's magnitude. Figure 1.4 shows a typical basis weight response to a slice lip which illustrates the separability of the two-dimensional (space and time) feature of a C D process. Depending on the type of arrays and of different controlled properties, the spatial responses could be very different. For example, the basis weight spatial response to a dilution actuator is usually much narrower than that to a slice lip actuator as shown in Figure 1.5. Similarly, the dynamic responses could vary from a very fast one to a very slow one. For example, the dynamic response of moisture to a rewet shower is very fast compared with the scan time (the rise time could be less than one scan time), and thus the process dynamics are dominated by the sensor filtering and the transport delay. The rise time of the dynamic response of caliper to an induction heating actuator could be more than 10 times the scan time while the dead time due to the transport delay is quite small as the scanning sensor is often installed just after the calender stack. Thus this rise time is the main factor of this process dynamics. Chapter 1. Introduction 7 OUTPUT RESPONSE PROFILE ACTUATOR INPUT PROFILE u(t) 0 CD [SPACE] TIME 0 CD[SPACE] 0 Figure 1.4: Two-dimensional feature of a C D process. 1.2.2 Multiple-Array System In order to address the interactions between the actuator arrays as qualitatively described in Section 1.1, a multiple-array C D process model has been proposed in [3, 4, 44]. It is a large scale model composed of single-array model blocks of the type in (1.3), Y(z) = G(z)U(z), (1.4) ut(z) Y(z) = , U(z) = (1.5) UN (Z) U Guhu(z) ... Gi hi {z) Nu Nu G(z) (1.6) GNylhNyl(z) -.. GNyN hNyN (z) u u where G^- and hij(z) are spatial interaction response matrix and dynamic response for each actuator array and each property, N and N are the numbers of actuator arrays and y u controlled sheet properties respectively. As introduced in Section 1.1, for a typical paper machine, there could be a total of as many as 600 individual actuators and 6000 profile measurements if all C D actuator arrays and all controlled properties are included in one control system. In other words, the dimensions for the multiple-array process model G(z) in (1.6) could be 6000 x 600. Chapter 1. Introduction 8 Figure 1.5: Typical steady-state responses of basis weight to a slice lip actuator and a dilution actuator for a A5gsm newsprint papermaking process. 1.2.3 111-Conditioned System The majority of industrial C D processes are severely ill-conditioned [38, 50, 69, 77]. In other words, the condition number of the spatial response matrix G 0 in (1.3) is large, |P « 1 T(Go) := where o(G ) 0 » and o(Go) denote the largest and smallest singular values of G M 0 respectively. Figure 1.6 shows the singular values of a typical C D process model for basis weight to slice lip actuators. Thus the gain of the C D process response is strongly dependent on the singular vector directions of the process. The smaller singular values corresponds to the weak or low-gain directions. The sign of the small gains of a C D process is usually not certain [29] and up to half the gains are vanishingly small even at steady-state u — 0 [77]. From Figure 1.6, it is evident that the controller needs to provide larger magnitudes of actuation to remove output disturbances that are aligned with these weakly controllable directions of the process (from the direction 40 to 64). A s stated in [74, 75], the small value (close to zero) of the smallest singular value Q_(GQ) is undesirable and may cause control problems, such as sensitivity to Chapter 1. Introduction 9 0.06 V - \ \ v.\ \ O b" in a \-. v D 13 > 3 DJ C .E co 0.02 - \ 1 ' 20 v. • 30 40 50 Singular vector direction j Figure 1.6: Singular values of a typical C D process model. uncertainty. If the controller has integral action, the closed-loop system could be unstable due to continuously erroneously increasing the magnitudes of the actuator actions in the low-gain directions which are not properly handled by the controller. A large magnitude actuation will more likely conflict with the physical constraints. Consideration of actuator limitations (constraints) will form an important part of this work. 1.2.4 Model Uncertainty Due to the process complexity and absence of physical modelling, the C D process model parameters are usually identified by input-output data [29, 37, 35]. The C D process spatial and dynamic responses are usually assumed as a constant transfer matrix and a linear low-order transfer function respectively. This assumption is only an approximation of the real industrial CD process. The presence of model uncertainty is thus inevitable in CD control. A n important class of process model uncertainties that plague the C D process can be characterized as global changes in the process model. In other words, deviations in the process model that influence the response of all actuators simultaneously. There are a variety of mechanisms for this and we present an overview of the main types here. • Sheet wander. The paper sheet may travel several hundred metres in the MD from the headbox to the scanning sensor (see Figure 1.1). The sheet can wander back and Chapter 1. Introduction 10 forth in the C D by several centimetres, corresponding to a change in the alignment between the inputs and outputs. • The base loading of an actuator array can alter the cross-directional spatial response shape. For example, increasing the average actuator setpoint base load to a steambox will alter the manner in which the steam leaves each individual actuator, thus resulting in a changed CD response shape. • Paper machine operators will sometimes speed up or slow down the paper machine for production reasons. Slowing the machine speed will have the effect of increasing the dwell time of the paper sheet travelling under each actuator array and altering the response of the sheet. • The overall slice opening is sometimes changed to improve the formation of the paper sheet. This action has a drastic impact on the spatial response of the weight profile to the slice lip actuators. To safely change this machine parameter requires having a C D control strategy that can handle the model change (typical strategies either use robustness or schedule the controller tuning for slice opening changes). • Many paper mills use the same set of controller tuning parameters for several grades of paper manufactured on a machine. Each grade of paper will respond differently to the actuators, as the weight, moisture and caliper targets, as well as the pulp stock properties, and machine speed may be significantly different. • A common form of global uncertainty is in terms of neglecting the multivariable characteristics of the paper machine (meaning that the off-diagonal blocks in (1.6) are neglected). For example, the slice lip illustrated in Figure 1.1 is typically used to control the basis weight profile of the paper sheet. However, it is well-known to papermakers that achieving a good moisture profile requires to first achieve a good weight profile. • Finally, the process models are identified from input-output data, and there are always unmodelled dynamics and features of the process that are not captured by the parameters of the process model. There exist other forms of non-global model uncertainty that will affect individual actuators such as: • Faults such as failure or partial failure of one or more actuators in an array. • Due to the physical construction of the paper machine, the paper sheet will typically respond differently to the C D actuators near the sheet edges that in the middle of the sheet. Chapter 1. Introduction 11 • As water is removed from the drying paper sheet, it may shrink in the C D resulting in the paper sheet being narrower at the scanning sensor than it was at the actuator array. It is true that no modelling technique can accurately capture all features of an industrial paper machine. Our structured model uncertainty described in Section 2.1.1 of Chapter 2 addresses the effects of the global uncertainty as described above. Conventional presentation of process model uncertainty is achieved by including "unknown" matrix perturbations A on the process model [75]. Model uncertainty can be either unstructured or structured [75]. Unstructured uncertainty is defined as a "full" complex perturbation A , where at each frequency any A(ju) satisfying o(A(ju)) < 1 is allowed. Structured uncertainty typically means that the perturbations can be rearranged into block diagonal form A = diag(Ai). In [28, 29, 82, 84], an additive unstructured model uncertainty based on uncertain singular values is used, while a structured model uncertainty description may be found in [33, 39, 40, 51, 52]. In [51, 52], the authors addressed the parameter uncertainty for both dynamic response and spatial response, however they did not explicitly address an important uncertainty - misalignment between the inputs 1 and the outputs in real-life C D processes. In [39, 40], the process model is described by rational transfer functions of spatial and dynamic Laplace variables and includes a structured uncertainty representation. In [33], a frequency method is used to show that the misaligned mapping between the actuators and the measurements results in a narrower achievable closed-loop bandwidth. 1.2.5 Coordination of Actuator Arrays As mentioned in Sections 1.1, one actuator array could affect more than one controlled property or one controlled property may be controlled by two actuator arrays. For example, slice lip actuators can affect basis weight, moisture and caliper. On industrial paper machines, the steam box (drying) and water spray (rewet) actuators are coordinated for controlling moisture [42, 72]. More recently, model predictive control for multiple-array systems in (1.6) has been proposed in [3, 4] for coordinating the actuator arrays and has proven to get much better control performance and energy savings than conventional C D control strategy (i.e., ignoring the interactions between the actuator arrays). Misalignment means that the actual CD locations of the outputs are not consistent with the predicted ones from the corresponding inputs. 1 Chapter 1. Introduction 1.3 12 Industrial C D Model Predictive Control A block diagram of Honeywell's industrial C D model predictive control ( M P C ) is shown in Figure 1.7. In Figure 1.7, the C D process block refers to the physical plant which needs to be controlled (including the actuator arrays, the sheeting forming process, and the scanning sensor(s), see Figure 1.1 for the various actuator arrays and scanning sensors). The M P C controller block refers to a computer on which the control software which is used for user interface and communication with other equipments has been installed. The real time Q P solver block is another computer which is used for calculating the optimal solutions and is always connected with the M P C controller computer. The'local area network ( L A N ) is used for communicating the computer with the actuators and the scanning sensors, because in a paper mill, they are often located far away. The model identification and tuning tool block represents the identification software [35, 36] and a software tool for tuning the performance and robustness of traditional feedback C D control [82]. The dashed line block refers to the software used to tune the C D - M P C parameters and predict the steady-state performance of the closed-loop system. The identification and tuning software is installed in a laptop computer which is connected with the M P C controller computer only when needed. Currently trial and error and closedloop simulations are mostly used for tuning the large-scale C D - M P C and predicting the steady-state closed-loop performance. The implementation procedure of the C D - M P C in the paper mill is as follows: first, the user turns off the closed-loop control of the process and uses the identification software to identify the process model. Then, based on the model, the user tunes the C D - M P C parameters and predicts the closed-loop performance until he/she is satisfied with the predicted closed-loop performance. This can be completed off-line. Finally, the user transfers the tuning parameters to the M P C controller for implementation. 1.4 1.4.1 Literature Overview Spatial Frequency and Rectangular P l a n t Spatial frequency methods have been used in C D control for some time. Because the actuators in an actuator array are identically constructed and evenly distributed across the sheet, the responses to these actuators can be considered as identical except for those near edges. If we have enough large number of actuators in an array, the C D process is almost spatially-invariant [18, 20]. In other words, the output profile can be obtained through convolution of the input profile and one single actuator's response [20, 33]. It is known that the convolution of two signals in time (space) domain can be transferred into Chapter 1. Introduction 13 Actuator setpoints Scanner sensor measurements LAN LAN MPC Controller Direct connection Real time QP solver Trial and error, Closed-loop simulations L A N connected when needed Model identification and tuning tool CD-MPC weights and closed-loop prediction Figure 1.7: The block diagram of Honeywell's industrial cross-directional model predictive controller. the production of two signals in the (spatial) frequency domain. Therefore, the output profile can be transferred from the spatial domain into the spatial frequency domain by the Fourier transform. The spatial frequency concept for web processes was thoroughly addressed by Duncan in [18]. The spatial bandwidth of CD control systems for web processes was investigated by Duncan and Bryant [20], reporting that the spatial bandwidth of the whole array of actuators depends on the spatial frequency response of a single actuator. The spatial frequency method brings an important advances for CD control by separating the controllable and uncontrollable components of a CD disturbance in an intuitive spatial frequency domain, because this method is directly associated with a physical variable (spatial frequency) unlike the singular value decomposition method, which is usually not linked with any intuitively understandable physical variables. A two-dimensional frequency domain loop shaping design technique was developed in [77, 83] and applied to CD controller design in [82]. This two-dimensional loop shaping method assumed that the spatial response matrix G in (1.3) is a (square) symmetric 0 Chapter 1. Introduction 14 circulant matrix [17], which can be diagonalized by the Fourier matrix. Thus, this method greatly simplified the controller design via decoupling the multi-input-multi-output C D system into a family of single-input-single-output transfer functions across the discrete spatial frequencies as in [6]. More recently, a spatial frequency anti-windup strategy for cross-directional control problems was provided by Goodwin et al in [34, 68]. Global closed-loop stability condition can be easily provided by this anti-windup spatial frequency method, while the steady-state performance attained with this method is close to the optimal Q P solution. In real industrial paper machine C D processes, the outputs generally significantly outnumber the inputs (by a factor between 3 ~ 10). These rectangular plants are usually transformed into square plants by first filtering the high resolution outputs with a low-pass anti-aliasing filter and then downsampling the outputs to the resolution of the inputs. The square C D model was used by a lot of C D control literature [46, 47, 51, 77, 82, 83]. However, squaring the process model before applying control techniques may be inappropriate for large-scale processes as pointed out in [69]. In addition, it is generally not possible to square the multiple-array process model G(z) in (1.6) when considering C D control with multiple actuator arrays and multiple controlled properties [3]. 1.4.2 C D M o d e l Predictive C o n t r o l ( M P C ) Model predictive control is a very successful control technique [63, 55] in process industry, especially in the petrochemical sector, due to its systematic capability for handling constraints on the inputs and outputs. However, it is less common in paper machine C D or other sheet forming process control. The main reasons include the large-scale and severe model uncertainty characteristics of C D processes as introduced in Section 1.2. There are mainly two methods for solving the large-scale problem of C D processes (up to 600 inputs and 6000 outputs). The first is to reduce the model size [20, 43, 44, 45, 50, 67, 71, 87]. A l l these model reduction methods use different basis functions to approximate the outputs and inputs profiles. In [20], Fourier transform is used to represent the controllable subspace in terms of spatial frequencies. In [43], spline functions are used to represent the sheet profiles. In [44, 67, 71], the principal component analysis method is used for model identification and reduction. In [45, 50], the G r a m polynomials or discrete Chebyshev coefficients are used for model reduction. In [87], the truncated singular values method is used for keeping robustly controllable components in the original model. The second method is to use a fast quadratic programming (QP) solver [3, 7, 13, 87] or using a linear programming (LP) solver [15, 16]. In [3], on-line computational time was shown to be greatly reduced by making use of the sparse structure of C D process models. In [7], an efficient Q P solver called QPSchur, which is basically an active set Chapter 1. Introduction 15 method using Schur complement algorithm, is developed for C D - M P C by re-organizing the constraints. In [13], the specific features of the system's state-space matrices are exploited to reduce the computational time for solving the Q P problem. In [87], the actuator constraint set is approximated by an ellipsoid and then an efficient algorithm is developed for solving the convex problem through scaling the size of the ellipsoid. In [15, 16], the linear programming method and modifications thereon have demonstrated their ability to solve large-scale control problems i n real time at the expense of some performance degradation relative to Q P formulation. Although there is significant literature about robust M P C for handling model uncertainty [5, 32, 48, 49, 53, 66, 70, 89], robust constrained M P C may be not implementable in real large-scale industrial processes due to intensive on-line computational requirements [8, 55]. The on-line computational demand of most of robust M P C algorithms grows exponentially with the problem size and thus is prohibitively high for C D processes [29]. 1.4.3 Robust C D Control Since Laughlin [52] in 1988 and Duncan [18] in 1989 analyzed the robustness of C D control systems, much has been published about robust C D control [10, 19, 51, 79, 82, 83, 88]. In [51], the authors used circulant matrix theory to develop methods for designing robust controllers based on the design of one single loop controller. The controller can be expressed as . . kAz)y{z), I u{z) = < I G , l y y ' k (z)y(z), 2 for decentralized control, ' for model-inverse-based control, (1.8) where ki(z) and k (z) are scalar transfer functions, G\ is the square spatial interaction 2 matrix. The decentralized controller's robustness is achieved by detuning its dynamics. One disadvantage is that forcing the controller k\(z) to be decentralized may introduce too much conservativeness because the well-controllable modes of the system have to be detuned. Another method presented i n [51] is model-inverse-based control. Both these techniques require well-conditioned interaction gain matrix Gi, while i n practice C D processes are usually ill-conditioned [20, 45, 50]. For plants where the singular vector directionality is independent of frequency (such is the case for the separable plant G(z) in (1.3)), the singular value decomposition (SVD) method can be used to decouple the system into nominally independent subsystems [46]. As evident from the model G(z) in (1.3), the C D process is such a plant because of the separation of the spatial response matrix Go and the dynamic response h(z). In [18, 19], singular value decomposition (SVD) of the system was used and the resulting set of n Chapter 1. Introduction 16 SISO systems was analyzed. The controller was designed as follows: u{z) K = K-k(z)y(z), (1.9) = {G G +rI )-'Gl eR T 0 Q n n x m , (1.10) where n and m are dimensions of the inputs and outputs respectively, K is the precompensator matrix. The pre-compensator matrix was detuned with a uniform penalty r in (1.10) so that the controller satisfied Nyquist stability requirements for all control modes based on limits on maximum and minimum singular values. This led to a degradation of the performance of the well-controllable modes of the system but to a lesser extent than the decentralized scheme in (1.8) of [51] (i.e. detuning the dynamics). In [10, 28, 88], the pseudo-SVD was used for robust C D control, which is defined as G{z) \(z) = G h(z) = UZV h(z) T 0 = UA(z)V , T '= T,-h(z) (1.11) (1.12) where U and V are real orthogonal matrix. The elements of the diagonal matrix A(z) are transfer functions. These diagonal elements Au(z) are referred to as pseudo-singular values [27]. The controller form was h(z) K{z) = V U. T (1.13) k (z) n In [28], a data-driven multivariable uncertainty identification algorithm was presented and a reduced-order pseudo-SVD controller based on confidence of each of the pseudo-singular values was proposed. This controller is robust because it does not try to control those directions whose corresponding pseudo-singular values' signs are not known with enough confidence. More recently, a robust loop shaping CD controller has been reported in [77, 80, 81, 82, 83]. The C D process is modelled as a linear, square, circulant symmetric system with a norm bounded additive unstructured model uncertainty. The control structure is as follows: u(t) = C(z)y{t) + Su{t-l), where C(z) is the 'control' part and S G R n x n (1.14) is a band diagonal matrix. If S = I, it is an integral controller and the closed-loop system could be unstable due to ill-conditioning and uncertain pseudo-singular values of the process model as mentioned above. Proper design of S can improve the robustness of the closed-loop control system. Details about Chapter 1. Introduction 17 how to design the blocks C(z) and S in (1.14) may be found in [56, 77, 82]. 1.4.4 Steady-State Closed-Loop Performance with Constrained CD-MPC The C D control application is particularly concerned with the appearance of the steadystate measurement profiles and actuator profiles. The main reason being that the CD process is typically so severely ill-conditioned that the familiar requirement of zero steadystate error cannot be imposed on the closed-loop [18, 20, 29, 45, 83, 92]. A significant number of the singular values of the plant model are vanishingly small and thus any component of the error profile which aligns with the corresponding singular vector directions must remain unattenuated. Model predictive controller (MPC) tuning is typically performed via simulations of the closed-loop system [12, 54, 76]. However, even with efficient M P C algorithms [7], it is still time-consuming and inconvenient to run such a large-scale simulation to evaluate the effect on the steady-state of any changes to the tuning parameters. It may require up to an hour for the simulation to evolve sufficiently to evaluate the steady-state. The difficulty inherent in this strategy then is in the development of a tuning tool that will be practical in an industrial environment. Rawlings et al [59, 60, 64, 65] proposed a target calculation method for steady-state of the inputs and outputs under an infinite horizon M P C structure. Unfortunately, their calculation method can not be used for finite horizon M P C used in the industrial strategy in hand [3, 4]. Wills and Heath [92] analyzed the steady-state behavior of CD control systems in a closed-loop for a class of CD controllers and gave a lower bound on output variations, but they did not explicitly provide a technique to directly calculate the steady-state input and output profiles with active constraints on actuators. 1.5 Aims and Contributions of the Work 1.5.1 Motivation and Objective As described in Section 1.4.2, M P C has advantages in constraint handling and actuator arrays coordination for CD processes. However, for an industrial C D - M P C such as described by Backstrom et al [3, 4, 7], there are some disadvantages including intensive computational time and the lack of guaranteed robustness in presence of the inevitable model uncertainty. Also, there are no currently consistent methods for designing optimization weights to achieve the desired performance and no convenient technique for predicting the constrained C D - M P C performance. Therefore, the main aim of this work is to solve Chapter 1. Introduction 18 the problems which exist in the industrial CD-MPC as outlined in Section 1.3. The detailed objectives include an effective and efficient model reduction method, an intuitive and consistent controller parameters analysis technique for performance and robustness of CD-MPC; and a fast tool for predicting the performance for multiple-array CD-MPC with actuator constraints. 1.5.2 Contributions The main contributions of this work are summarized as follows: 1. The development of a novel effective and efficient model reduction method for CD processes based on wavelet packets. This method is equally applicable to single-array and multiple-array systems. This technique was incorporated into the industrial CDMPC and greatly reduced the on-line computational time while maintaining close to optimal closed-loop control performance. 2. The new concept of rectangular circulant matrix is proposed and its properties are proven. An intuitive two-dimensional frequency technique based on rectangular circulant matrices is applied to CD processes including the open-loop plants, the structured model uncertainties, and the closed-loop transfer functions. The nonsquare nature of CD processes has rarely been preserved in industrial control design because of previous inability to perform two-dimensional frequency analysis. 3. The performance and robustness of the linear response of the CD-MPC is analyzed in the two-dimensional frequency domain. Properly tuning the CD-MPC controller parameters and considering the practical structured model uncertainty is shown to greatly improve the control performance while preserving robustness. 4. A practical performance-prediction tool based on the sensitivity function technique is developed for constrained CD-MPC. This prediction technique has reduced the computation time by up to two orders of magnitude when compared with the conventional method while accurately predicting closed-loop performance of the CD-MPC. 1.5.3 Thesis Overview The thesis is organized as follows: Chapter 2 first presents CD process models and the industrial CD-MPC structure. Then, the CD-MPC controller design requirements for model reduction, performance and robust stability, and constrained CD-MPC performance prediction are addressed. Chapter 1. Introduction 19 Chapter 3 introduces rectangular circulant matrices and their properties. These rectangular circulant matrices can be (block) digaonalized by the Fourier matrices. A twodimensional frequency technique based on rectangular circulant matrices is used to decouple the large MIMO transfer functions such as the open-loop plants into a family of SISO transfer functions. Chapter 4 presents a modified wavelet packet method for model reduction. Some practical and important issues such as correction of distortions are investigated. This method is extended to multiple-array systems. The implementation of the C D - M P C in the reduced domain is described. Two simulated examples based on industrial paper machine C D processes are given for demonstrating the effectiveness and efficiency of this method. Chapter 5 presents the robust stability for the unconstrained C D - M P C closed-loop system with structured and unstructured uncertainty. The robust stability condition for multiple-array systems is also derived. The roles of tuning parameters are described. Three examples show that properly tuning the parameters, choosing the structures of the optimization parameters, and considering the structured model uncertainties may greatly improve the closed-loop performance of the control system. Chapter 6 describes the development of a static optimizer for predicting the closedloop steady-state performance of the CD-MPC, especially for the constrained case. The method is based on the minimization of the difference between two sensitivity functions, one is for the static optimizer and another is for the conventional C D - M P C . Two examples are given for demonstrating the efficiency of the proposed method. Finally, conclusions are given and future research directions are discussed in Chapter 7. Chapter 2 Problem Statement In this chapter, the industrial CD process models, for a single-array system and a multiplearray system are introduced. Then, a particular industrial model predictive controller for these processes is described. Finally, the requirements for this work are specified. 2.1 2.1.1 Process M o d e l Single-Array Process Model The standard paper machine C D process model for a single-array system is used in this work and is given by y(z) = G(z)u(z)+d(z), G(z) = G -h(z), h(z)= 0 where y(z) e C , u(z) G C " , and d{z) e C m m . _ 1 Z d _ (2.1) (2.2) i ; (XZ are the Z-transforms of the measurement profile, the actuator set-point profile and the disturbances respectively. G € R ™ is mx 0 the interaction matrix which describes the spatial response of the process, h(z) is the Z-transform of the first-order with dead-time dynamic response of the process. The 1 dimension m of the measurements y is typically much higher than the dimension n of the actuator set-points u (usually 3n < m < lOn). In addition to the dynamics, each of the n actuators has an influence on the spatial response of the paper profile. As presented in [36, 37], the parameterized structure cos (2-3) models the effect of an actuator on the continuous paper profile, where g(x) is a continuous function, the spatial coordinate a; is a scalar real number, the gain parameter 7 , the attenuation parameter a, the width parameter £, and the divergence parameter 8 are It is very likely that the deadtime Td is not an integer number of sample times. For such a case, the transfer function h(z) may include a non-minimum phase zero in the numerator. The uncertainty of the process dynamics due to the presence of non-integer deadtime may be included in the unstructured uncertainty which will be defined in (2.7). 1 20 Chapter 2. Problem Statement 21 scalar parameters which are used for describing the spatial response shape. In practice, this continuous medium is measured by a scanning sensor at m evenly spaced positions across the sheet. As there are n actuators, the discrete process model can be given by G = [g g ,--. 0 u ,g }ER , (2.4) mxn 2 n where column p. e R represents the sampled spatial response to the i actuator given by , m th gi(X) = g(X - a), X = l,2,---,m (2.5) where the alignment parameter is the measurement point coordinate of the spatial response centre for, the i actuator (it is also known as the alignment parameter), X is the measurement point, g(X — C;) is the sampled version of g(x) in (2.3). For illustration purposes throughout this thesis, we will utilize an example model of a slice lip actuator array used to control the basis weight profile of a newsprint paper machine. The parameters of the nominal model G(z) in (2.2)-(2.3) are 7 = 0.0084, a = 1.56, i = 19.5, 3 = 0.13, a = 0.3679, T = 2, and d = 2.5 + 5(i - 1), i = 1, • • • , 64. The dimensions of the outputs and the inputs are m = 320, and n — 64. From equations th d 0.01 0.005 o (ft Ui - 0 A / \ \ 0 i \ -0.005 -0.01 50 100 150 f (a) 1 \ 200 V i 250 300 Scanner measurement locations [measurement points] (b) ^ o §1 •5 20 30 40 Actuator number C3_ eT 20 30 40 Singular vector direction j Figure 2.1: (a) Columns of the spatial response matrix Go in (2.4) identified by a basis weight process controlled by n = 64 slice lip actuators; (b) Corresponding slice lip actuators positions; (c) The singular values of Go (Data courtesy of Honeywell). (2.1)-(2.3), the nominal model G(z) can be constructed. Figures 2.1a and 2.1b illustrate Chapter 2. 22 Problem Statement the nominal steady-state response of the basis weight profile to slice lip actuators. The spatial interaction matrix G in (2.2) is a band diagonal matrix as shown in (2.6) when very small values of g in (2.5) are truncated. The non-zeros of Go in (2.2) are illustrated in Figure 2.2. This band diagonal and nearly spatially-invariant structure of the spatial interaction matrix Go in (2.6) is very important and extensively used in Chapters 3 and 5. 0 t 0i 96 pll • • 026 0 ••• 92 9s 010 • • 025 0 ••• 93 9A 99 • 024 0 ••• 9A 93 9s • • 023 0 ••• 95 92 97 • 022 0 ••• 96 0i 9e • • 021 026 97 92 9& • • 020 025 • • 921 016 • • 01 06 ••• 0 922 017 • • 02 05 ••• 0 923 018 • • 03 04 ••• 0 92A 019 • • 04 03 ••• 0 925 020 • • 05 02 ••• 0 926 021 • • 06 01 ••• 926 0 0 025 • • 010 05 ••• 0 0 026 • • 011 06 ••• 0 0 0 ••• 0 • • 0 Following the traditional approach used in robust control for representing model uncertainty [75, 93], we assume that the true process response belongs to a set of possible response models, described by an additive unstructured or structured perturbation on the nominal transfer matrix model G(z) in (2.2). The additive unstructured model uncertainty A (z) is given by Gl Gp^en.! := {G{z)+A {z)-.o(A {e ™)) Gl Gl a <l(u)}, (2.7) where 1(UJ) is a positive scalar function which limits the perturbed transfer matrix G (z) to a neighborhood of the nominal model G(z) in (2.1); o denotes the maximum singular value, z — e and the symbol UJ denotes the dynamic frequency. On the other hand, the effects of the global model uncertainty described in Section p l2mj Chapter 2. Problem Statement 23 0 20 40 60 Figure 2.2: The spatial interaction matrix G in (2.2) has a band diagonal structure. 0 1.2.4 of Chapter 1 can be addressed by the following additive structured model uncertainty, G {z) e n p 9i ( ) x P g (x) p P lp =9p{X = • • • ,g ]h(z)} 2p np = G{z) + A (z), G2 - Ci ), P <e = 7(1 ^ + <5), 8 ,5 8c,5 3 1 := {[g ,g , cos —-—• +e cos ^ Ci ~\~ 8c, Cip 7 2 cn l X — 1, 7 Op = a(l + 5 ), C = ^(1 + <y, /? = a e [Si,Si],S e ,m, i = P + S ), p [-5 m/n,8 m/n], c ••• P 2 2 1, • • • ,n, where <?(:r) denotes the parameterized structure g(x) in (2.3) with parameter uncertainty, g (X c ) is the sampled version of g (x) with alignment uncertainty, S ,8 ,8a,8^8j3 denote the parameters uncertainty, and 8 are parameters uncertainty limits. A matrix with a large condition number T(Go) in (1.7) is said to be ill-conditioned [75]. The majority of industrial CD processes are severely ill-conditioned [45, 50]. Indeed, because the minimum singular value a(Go) is very close to zero for most CD processes [20, 29, 45, 83], the ratio T(Go) in (1.7) between the maximum and minimum singular values could be very large. Figure 2.1c illustrates that about 20 smallest singular values of G(e ) in (2.2) are almost zero at steady-state OJ — 0. The severe ill-conditioning plus the model uncertainty A i(z) in (2.7) and A (z) in (2.8) make the CD control problem challenging. In fact, even the sign of small gains of the process model are usually not known [29] due to model uncertainty, which is rare in general control systems. If not properly handled, this steady-state sign uncertainty may destabilize the control system. p p — ip p c 2 t2mj G G2 1 Chapter 2. Problem Statement 24 A state-space model is much more convenient than a transfer matrix model for the forthcoming controller development, and Equation (2.1) can be rewritten in the following form, x(k + l) = Ax{k) + Bu(k), (2.14) y{k) = Cx(k)+d(k). (2.15) The generation of matrices A, B, C in (2.14) and (2 15) from the transfer matrix equations ? (2.1) and (2.2) is standard and may be found in Appendix D.2. A n augmented state-space model is used in the forthcoming model predictive control (MPC) framework in order to work in terms of changes to the actuator set-point profile Au(k) as suggested in [54, 58], x (k + l) = y(k) = a A Xa(k) a + (2.16) B Au(k), a (2.17) C x (k)+d(k), a a where x {k) = A 0 CA I a Ax(k) Au{k) = u(k) - u(k Cx{k) B B„ - CB 0 - 1), (2.18) / (2.19) To clarify notation, the operator A followed by a time-indexed variable denotes the first difference in time, such as Au(k),Ax(k) in (2.18), and the operator A with subscript is used to denote model uncertainty as A i{z) G 2.1.2 in (2.7) and AG2{Z) in (2.8). Multiple-Array Process Model In current industrial paper machine C D control systems, there are often one or two actuator arrays for controlling one of the paper properties (basis weight, moisture or caliper). However, one actuator array will affect more than one property. For example, slice lip actuators are originally designed for controlling basis weight, but they affect moisture and caliper too. Hence, it is necessary to coordinate all actuator arrays for controlling all paper properties [3, 4, 72]. In this work, we also use the following multiple-array C D process model: Y(z) = G(z)U(z) + D(z), (2.20) Chapter 2. Problem Statement 25 G h {z) n ... n G h (z) lNu 1Nu G(z) (2.21) GN \h y(z) y ... Ny G N y N u h N y N u (z) where Y(z) = 3/1(2), U{z) = ui(z), D(z) = di(z), (2.22) VNy(z) r , eC ^-^ ( ) x l , (z)eC f (2.23) n Uj d {z) (2.24) Ny z~ i Tdl Gij e R f hij{z) mxn where N y and N u ^ij 1 , 1= ,N , j - I,-- - ,N , y U (2.25) are the numbers of controlled properties Y(z) and actuator arrays U(z) respectively, G^ and hij(z) are the sub-plants' spatial interaction response matrices (constant transfer matrices) and dynamic responses (scalar transfer functions) respectively as defined i n (2.2), D(z) are the disturbances, T ij is the time-delay, a*,- is the open-loop d pole. The number of controlled properties 1 < iV < 5 and the number of the actuator y array 1 < N < 5. Note that a l l the controlled properties have the same dimension u m because they are usually measured by the same scanning system while the actuator arrays have different dimensions denoted as n i , • • • ,n . Nu The state-space model for the multiple-array system is a straightforward extension from the single-array system and will be discussed in Chapter 6. 2.2 Industrial Cross-Directional M P C It is well known that, due to its large-scale, the paper machine C D process can pose computational load problems for the on-line optimization used in model predictive control ( M P C ) . This thesis considers an industrial M P C controller for C D processes developed by making use of the sparse structure of G i n (2.2) [3] and choosing an efficient quadratic 0 programming (QP) solver [7]. In this industrial M P C , replacing y in (1.1) by the target paper profile y , the M P C problem is posed as follows, sp r minJ(fc) = H P min \y2(y(k+j)-y (k+j)) Q (y(k sp T + 1 H u - l Y, j=o j)-y (k+j))+ sp } [Aw(fc + j) Q2Au(k T + j) + u(k + j) Q u(k T 3 + j)] \ ) (2.26) Chapter 2. Problem Statement 26 subject to the augmented state-space equations (2.16) and (2.17) and the following constraints on the actuator array: Umin — -Udmax < li(A^) aV < max (2.28) Ali(fc) < Udmax, Uavgmax — U g{ki) ^ -M (2.27) Umax, B m U g , aV max (2.29) (2.30) • U{k) < Mmax, where y{k + i) G R : output predictions; m y G R : set-points for controlled properties; m sp H : prediction horizon; p H : control horizon; u Qi G R diagonal weighting matrix for controlled properties; m x m : Q G R ™ : diagonal weighting matrix for the changes in the actuator array; xn 2 Q €R n x n Qi,Q , 0,3 are positive semidefinite; 3 2 : weighting matrix for the actuator array; U in, U m E R : lower and upper actuator limits respectively; n max M ax G R " : bounds for bending limits; m Udmax G R : maximum change between successive control actions; n u vg(k) = a e = [1, • • • , 1] G R : average of one actuator array; n U, : maximum value allowed for the average of one actuator array; avgmax B m GR n x n : band diagonal bending moment matrix, which is given by -1 1 0 1 -2 1 0 1 -2 0 Br, (2.31) 0 -2 1 0 1 -2 1 0 1 -1 The constraint (2.30) is typically used only for slice lip actuators in Figure 1.1 as the actuators involve bending a flexible metal lip and need to limit bending to avoid permanently deforming the expensive lip. The remaining constraints (2.27)-(2.29) are typically required for all other actuator arrays. The constraints (2.27)-(2.30) are set up based on mechanical considerations of the paper machine. Chapter 2. Problem 2.3 27 Statement Model Reduction The MPC controller described in the previous section is designed for single-array systems but can be naturally extended to multiple-array systems. Due to the huge dimensionality in such systems (up to 600 inputs, 6000 outputs, 1800 hard constraints for a multiple-array system) and relatively short sampling interval (from 15 seconds to 30 seconds), solving a constrained QP problem such as (2.26)-(2.30) during the sampling interval is challenging. There are two ways to solve this problem: one is to find an efficient QP solver, another is to effectively reduce the model dimensions while maintaining good control performance. Both methods may be used together for solving some extreme cases, e.g. when having a very large number of actuators with very tight constraints on inputs. Because there is an efficient QP solver [7] for the above industrial CD-MPC, this thesis will focus on the model reduction method. The objectives of this part of the work include: 1. developing an effective and efficient model reduction method for a large, sparse, illconditioned band matrix such as G in (2.2). Detailed requirements for the model reduction are as follows: 0 • Dimensionality reduction; • Preservation of controllable components; • Elimination of non-robustly controllable components; • Reduction in the condition number of the spatial interaction matrix; • Sparsity of the reduced model; • Applicability to multiple-array systems. 2. developing an MPC for paper machine CD processes based on the reduced model. Detail requirements for this reduced model based MPC are as follows: • Fast enough for implementation; • Good closed-loop performance. 2.4 Robust Stability and Performance The output predictions y(k) are calculated based on the nominal process model in (2.1)(2.2). However, the true process response will differ from the nominal response and the model uncertainty such as A in (2.7) and A in (2.8) must be incorporated in an industrial control strategy. G 1 G 2 Chapter 2. Problem Statement 28 The design degrees of freedom available in the industrial model predictive controller are the weighting matrices Qi,Q , Q3, the prediction horizon H , and the control horizon 2 p H . These parameters should be designed such that: u • the closed-loop control system given by (2.26)-(2.30) and (2.1)—(2.2) is stable for all G (z) 6 I i i in (2.7) for the unstructured uncertainty or G (z) e U in (2.8) for the p p 2 structured uncertainty; • the paper quality index K in (1.1) is kept as small as possible. 2.5 Steady-State Performance Prediction The industrial M P C controller is tuned via the optimization weights Qi,Q , Qs and the 2 horizons H and H in (2.26). The standard tuning approach is to select some values for p Qi,Q2,Q3,H u p and H , then evaluate the expected closed-loop performance through offu line simulations. Through trial and error, values of the tuning parameters are determined for acceptable closed-loop performance. For conventional multivariable problems, the system is of small enough dimension to make this approach feasible and efficient [54]. However, to run a closed-loop simulation to steady-state for a multivariable CD control problem often takes several minutes (up to 60 minutes) due to the large-scale characteristics of G(z) in (2.2). This simulation would need to be run each time the user modified one of the tuning parameters in the performance index (2.26). In order to be used as a component of a practical tuning tool in an industrial environment, a prediction of the closed-loop's steady-state performance must be available quickly after the user modifies a parameter of the optimization problem in (2.26). A static optimizer is proposed to solve this problem in Chapter 6. The requirements for the static optimizer are: • significantly faster than the usual closed-loop simulation method; • optimally predicting the steady-state closed-loop performance; • applicable for multiple-array systems. Chapter 3 Two-Dimensional Frequencies in CD Control Circulant matrices and the Fourier matrices have been used for C D processes (square models) in [34, 47, 51, 77, 82, 83] in order to simplify the controller design. This previous literature motivates us to propose the novel concept of rectangular circulant matrix [23, 24], which can be applied to rectangular C D processes. This chapter first introduces closed-loop transfer functions with an unconstrained CDM P C given by (2.26) and (2.1)-(2.2). Then, the definition of a rectangular circulant matrix (RCM) and its properties are described in Section 3.2. The relationship between the R C M and two-dimensional frequency is described in Section 3.3. Then, the closed-loop transfer functions are proven to be approximate RCMs in Section 3.4. Finally, the openloop plant G(z) in (2.1), the structured uncertainty A (<2) G2 in (2.8), and the closed-loop transfer functions are analyzed in the two-dimensional frequency domain in Section 3.5. 3.1 Closed-Loop Transfer Functions Strictly speaking, the problem as described in Section 2.4 of Chapter 2 would require analyzing the robust stability of a closed-loop control system containing a nonlinear optimization. While such a solution would be beneficial [48, 49], it is quite challenging for CD processes due to their large-scale characteristics [8, 55]. In this section, we propose a technique based upon a linear approximation to the problem. We proceed as follows: 1. We first ignore constraints (2.27)-(2.30) such that the closed-loop system given by (2.26)-(2.30) and (2.1)-(2.2) is linear; 2. We then compute equivalent transfer matrices of the C D - M P C in (2.26); 3. We compute the closed-loop transfer matrices and analyze the performance and robustness of the system; 4. Finally, we re-introduce the constraints (2.27)-(2.30) for implementation and use the static optimizer, which will be presented in Chapter 6, to check the performance again. The subsequent development in Section 3.5.3 and Chapter 5 will justify the use of this approach in analyzing the expected performance of a linear closed-loop system. Numerous 29 Chapter 3. Two-Dimensional 30 Frequencies in CD Control simulations have shown that this technique is useful for evaluating the design of the optimization parameters for the problem in (2.26) even with active constraints. The model predictive control strategy that is used to optimize the performance index (2.26) can be split into two components: prediction and optimization. If the constraints (2.27) -(2.30) are neglected, then the system is linear and we may write down a closed-form solution for these operations [54]. Figure 3.1 illustrates the block diagram structure for such a system. In the figure, z~ is the unit-delay operator. Following common practice z l will also denote the .Z-transform complex variable and it should be clear which is applicable from the context of its use. Figure 3.1: The closed-loop system block diagram with the unconstrained M P C in (2.26). The M P C controller optimizing (2.26) is equivalent to the linear control law, Au(k) = K -y x (k) = L-y{k) a r s p + K -x (k) x + K ^-u(k-l), a (3.1) u + (I-L-C ){A x (k-l) a a + B Au(k-l)}, a a (3.2) where the constant matrices A , B , C are components of the augmented state-space a a a model in (2.19). The constant matrix L is the observer gain and its default value is L = x(k) [0%x , m I xm] T m in (2.14), 0 nxXm for the industrial M P C controller, here n is the size of states x denotes a n x m zeros matrix and I x mxm denotes a m x m identity Chapter 3. Two-Dimensional Frequencies in CD Control 31 matrix. The constant matrices K £ R , K £ R " ( " * + ) , and K _ £ M " are found from the optimization of the performance index ( 2 . 2 6 ) and are functions of the process model in ( 2 . 1 4 ) - ( 2 . 1 9 ) as well as the horizons H , H and the weighting matrices Qi, Q , Qz in ( 2 . 2 6 ) . The derivation of K , K , K _ may be found in Appendix A. The constant matrices and delay elements of Figure 3.1 can be re-arranged into the transfer matrices illustrated in Figure 3.2 with some tedious algebra. The constant prefilter matrix K £ R i preserved from Figure 3.1 and the complex transfer matrix Ki(z) € C and the constant transfer matrix K £ R may be found in Appendix A. N X M X r N X X u p r M x r x u N t u 2 x M s nXn n x m y d(z) u(z) ^ G (z) Ki(z) + p + u y(z) Figure 3.2: Equivalent transfer matrix representation of Figure 3.1. The linear closed-loop system in Figure 3.2 is robustly stable for all plants Gp(z) in (2.7) and (2.8) if it is nominally stable and p{T (z)A (z)) ud < 1, G (3.3) T {z) = K(z) • [/ - G(z) • K(z)]-\ (3.4) ud K(z) = Ki(z)K , (3.5) y where p denotes the spectral radius, A (z) is the unstructured uncertainty A (2;) in ( 2 . 7 ) or the structured uncertainty A (z) in ( 2 . 8 ) and the result ( 3 . 3 ) follows from the small gain theorem [75], the control sensitivity function T (z) £ C is the nominal closedloop transfer matrix connecting the actuator set-points u(z) to the output disturbance G G1 G2 ud nxm d(z). As CD control is largely a regulation problem with only occasional set-point changes, the performance criterion is a function of its disturbance rejection properties. The nominal transfer matrix linking the measured profile y(z) with the output disturbances d{z) is Chapter 3. Two-Dimensional Frequencies in CD Control 32 obtained from Figure 3.2 as, T {z) = yd where T (z) e C m x m yd [I-G(z)-K(z))-\ (3.6) is also known as the sensitivity function [75]. Note that T (z) and ud T (z) themselves should be stable in order to get the nominally stable closed-loop system. yd For stable plants like C D processes, properly choosing the control and the prediction horizon can guarantee that the transfer matrices T (z),T (z) ud yd are stable [30, 55]. The two transfer matrices T (z) in (3.4) and T (z) in (3.6) will be used to analyze yd ud the behavior of the linear closed-loop system illustrated in Figures 3.1 and 3.2. However, as T (z) e C n x m ud and T {z) E C mxm yd have very large dimension (n < 300,m < 2000), the study of their properties is difficult. A two-dimensional frequency domain analysis approach will be presented in the following sections for solving this problem. 3.2 Rectangular Circulant Matrices In this section, we first introduce circulant matrices, and then propose the new concept of a rectangular circulant matrix (RCM) and show the properties of RCMs. A circulant matrix is a special square matrix which is defined as follows [17, 41]: Co Cl c • • C„_i Cn-1 Co Cl • • Cn-2 Cn-1 Cl c 2 2 Co • C _2 n • Cn-3 c • • Co 3 where each row is a cyclic shift of the row above it. The properties of circulant matrices may be found in [17]. Here, we propose the definition of a rectangular circulant matrix (RCM): Definition 1 A matrix N € C mxn is a rectangular circulant matrix if it satisfies the following conditions: 1. If m — n • n, then column j = column (j — 1) circularly down shifted n a a times, j =2,3,--- , n , 2. If n = n • m, then row k = row (k a 2,3,-- - ,m, for a positive integer n . a 1) circularly right shifted n a times, k = Chapter 3. Two-Dimensional Frequencies in CD Control 33 Rectangular circulant matrices (RCMs) are analogous to square circulant matrices (SCMs) with respect to several properties: 1. If N N £C u m x are both R C M s , then N = Ni+ N is an R C M ; n 2 2. If Ni £ C 3 m i x n i and N £ C m i X n 2 2 2 are both R C M s , and either n = m (N is a S C M ) or n = rn (N has the same dimension as JV,), then N = N?N 2 2 l 2 £C n i X n 2 2 4 2 is an R C M ; 3. If JV! € C m x n is an R C M with m = n and JVi is invertible, then N = J V f is an 1 5 RCM; Properties (l)-(3) show that the class of R C M s is preserved under the matrix operations most commonly required for feedback control. A well known property of S C M s is their diagonalization with the Fourier m a t r i x [17], 1 F = — 1 1 1 Q- 1 o- -(n-1) 1 2 1 Q~ ~ {n -2(n-l) Q~ l) Q~ 2{n l) Q~ ( -^ 3 -(n-l)(n-l) £" (™- ) n 4 , (3-8) 1 where g = e ™ ,i = — 1. 2 A n analogous property holds for the R C M s defined above : (4) The matrix N £ C m x n is an R C M if and only if N = F NF" is a matrix of diagonal m blocks: [Ni, N = ••• , N ] na if m = T n -n a (3.9) if n = n -m Nn a a where H denotes conjugate transpose, each N k £ c m i n ^ '"^ m x m i n ^ '^ m is a diagonal matrix for k — 1, • • • , n . a For S C M s , the advantage of the diagonalization with the Fourier matrix F in (3.8) is n that the eigenvectors and singular vectors of all S C M s can be expressed in terms of the harmonic functions [47, 83]. While eigenvalues have no meaning for a rectangular matrix, an analogous result holds for singular values. 1 The Fourier matrix can be obtained by the Matlab command F = n fft(eye(n))/sqrt(n). Chapter 3. Two-Dimensional (5) If TV G C mxn 34 Frequencies in CD Control is an R C M , then N = F NF£ has the diagonal block structure of m (3.9), and a^N) = (N) 0j /£ = N(j + (k-l)n,j) \im = n -n a k=l (3.10) < ' £ N(jJ + if n — n • m (k-l)m) a k=l j = !,-•• ,min(m,n),fc= l,---,n . 0 Note that the singular values rjj(iV) in (3.10) are not necessarily sorted in descending order as usually presented. It is known [17, 77] that if JV e C is a SCM, then the diagonal elements of N = nxn F NF„ n occur in conjugate pairs. A n analogous property holds for the RCMs, (6) If N G C mxn is an R C M and n G C (q = max(m,n)) is a vector which represents q d the diagonal elements of N • • •, N iy AT: na Ni 0 0 N 0 in (3.9), 0 0 2 , n := d (3-11) diag(Af), 0 then, n (j d + 1) = n (q -j + d 1) , H j = 1, • • • , k, (3.12) where k = q/2 — 1 for even q and k = (q — l ) / 2 for odd Figure 3.3 shows nonzero elements of a SCM JV G C 3 2 0 x 3 2 0 S , an R C M i V e C r 3 2 0 x 6 4 and the nonzero diagonal elements of N and N . See Appendix C for the proofs of the above s r properties. 3.3 Frequency Domain Interpretation Given the spatial signals y G C m and u G C", then the formulas y = F -y m and u = F -u, n (3.13) Chapter 3. Two-Dimensional 35 Frequencies in CD Control 20 40 60 (c) 0 20 40 60 (d) Figure 3.3: (a) Nonzero elements of a SCM N ; (b) Non-zeros of the diagonal elements of N ; (c) Nonzero elements of an R C M N ; (d) Non-zeros of the diagonal elements of N . s a r where multiplication by matrices F m r and F , defined in (3.8), is equivalent to performing n the standard m-point and n-point discrete Fourier transforms respectively. The usual interpretation is that this is a decomposition of the (spatial) signals into spatial harmonic functions with m (respectively n) equally spaced frequencies between Q = 0 and = 2TT [62]. We will consider the case m = n • n in Definition 1 in some detail (the case n — n -m a a is a straightforward extension). If y = N • u, then y = N-u, with y and u in (3.13) and N = F NF". m (3.14) From (3.9) we find that for an input signal u of a single spatial frequency Cl (j), u (3.15) Chapter 3. Two-Dimensional Frequencies in CD Control 36 then the output V(j + (k- l)n) = ( " ^ 0, { j+ { ~ k 1 K j •^ ) k = ^ - ^ otherwise (3.16) 1 is composed of n spatial frequency components, a r (j ly + 27T / i — 1 \ [+ (fc-l) , (k-l)n) = — J A; = 1,2,.-. , n . 0 (3.17) Each of the discrete spatial frequencies 0 „ or Q may be converted into engineering units y with the usual formula I 2h( - )> 7T < < 27T ' 2lT n y where X denotes the sampling interval, and typical units for the spatial frequency v are For the C D control systems in general the input and output sampling in- cycles/metre. tervals X ^ X due to the fact that there are typically several times more measurements u y than actuators distributed across the same spatial domain. Two important spatial frequencies are the spatial Nyquist frequency of the measurements and the actuators defined by "» = 2^;, = W' (3-19) ( 3 y - 2 0 ) where X and X denote the sampling intervals of the outputs and the inputs respectively. y u Denote the j column of N e C th as % £ C , and define the Fourier transform of m x n m rij as n =F -n , j where the complex Fourier matrix F m nature of the columns of an R C M , then as / = 1^71 f ° ra n v m (3.21) j is defined in (3.8). Note that due to the circular —• for any j,k = !,••• ,n and is denoted — 3 = 1) • • • > , where / € M contains the magnitude of the spatial n m frequency components of the columns of N. A bandlimited / is one for which there exists an index j < floor(m/2) such that b /(?') = 0 for j b < j <m-j . b (3.22) ' J Chapter 3. Two-Dimensional Frequencies in CD 37 Control The discrete indices in (3.22) correspond to a physical spatial bandwidth v B 1 If N e C mxn Theorem limited to v B < is an from (3.18), RCM, y = N • u, and the spatial bandwidth in (3.23) is (the spatial Nyquist frequency of input u in (3.19)), then, I 0, otherwise where I m — n, j = floor(n/2) + 2, • • • ,n Proof: see Appendix C. This result indicates that when the column vectors of Af are bandlimited to the spatial Nyquist frequency of the inputs u , the matrix operator N can be replaced by a family N of scalar equations, one for each discrete spatial frequency Cl (j) with j = 1,2, • • • , n in u (3.15). Corollary 1 If N e C m x n is an RCM and v < v% [in (3.19), (3.23)), then B cr,(iV) = (3.26) |iV(z,j) where i = j for j = 1, • • • , floor (n/2) + 1, ori—j+m-nforj = floor(n/2) + 2, • • • ,n. Proof: From Theorem 1, the elements of N(i,j) are zeros except i—j and i = j +m — n, then (3.26) follows from (3.10). Remark: the singular vectors of Af in Corollary 1 can be given by the complex Fourier matrices F m and F defined in (3.8). n Based on the RCM's properties and Theorem 1, the following theorem is derived: Theorem 2 For the closed-loop system in Figure 3.1, if 1. G e R 0 m x n in (2.4) is an RCM, and 2. the columns g^ in (2.4) are bandlimited with v < ^ ^ (3-19) and (3.23), and n B 3. the weighting matrices Q e R y m x m , Q ,Qz € M 4- the observer gain matrix L = [0 eI] 2 T n x n in (2.26) are RCMs, and with e £ (0,1] in Figure 3.1, Chapter 3. Two-Dimensional Frequencies in CD Control 38 then, transfer matrices defining the open-loop plant model G(z) closed-loop transfer matrices T^z) GC n x m G C in (3.4) and T (z) G C mxm yd m x n in (2.2), the in (3.6) are also RCMs. Proof: See Appendix D. 3.4 Approximate Rectangular Circulant Matrices Theorem 2 shows that the R C M structure would be preserved for each of the transfer matrices important to the problem at hand. The optimization weighting matrices Qi, Q , Q3 2 in (2.26) are each a diagonal matrix and therefore satisfy Condition (3) in Theorem 2. 1 The industrial MPC has an observer gain matrix L satisfying Theorem 2 [3]. However, for practical CD control systems, the spatial interaction matrix G in (2.2) is not an R C M 0 according to Definition 1. There are two reasons for this. First, there is no guarantee that the number of measurements m in a scan is an integer multiple of the number of actuators n [29, 50]. Second, even assuming that m = n • n for integer n , there is still the problem a a that each column g^ of G in (2.4) is not equivalent to a circularly downshifted vector of 0 Pj-i- Each column g^ (or j — l , - - - ,n models the steady-state spatial response of the paper sheet to a unit step in the j th actuator. The spatial shift invariance holds for the actuators near the centre of the sheet, but does not hold for actuators near the edge of the sheet whose response is modelled as the truncated response to the last actuator in Figure 2.1 (also see (2.6)). In Figure 3.3(c), we know that all the column vectors are circularly shifted based on the RCM's definition, but in Figure 2.2, the column vectors are circularly shifted except the column vectors near the edges. Comparing Figure 2.2 with Figure 3.3(c), these two figures are exactly the same except for the upper right corner and the lower left corner. In other words, the interaction matrix Go in (2.2) is spatially invariant except for the edges. Before exploring the effect of this violation of Condition (1) in Theorem 2, it should be mentioned that those columns 7Jj in (2.4) that do satisfy the spatial shift invariance (of Definition 1), will typically also satisfy Condition (2) for the majority of practical C D control systems. Section 3.3 introduced the attractive properties of RCMs in terms of the ability to decouple across "spatial" frequencies. If we are to exploit two-dimensional loop shap- ing techniques [82, 83] for controller analysis and design then it is necessary to use a rectangular circulant approximation as the "true" system in (2.2). The weight matrix Q 3 could be equal to q^B^Bm (53 is a real scalar number) and then it is an approximate RCM, B is the bending moment matrix in (2.31). x m Chapter 3. Two-Dimensional 3.4.1 39 Frequencies in CD Control Spatial Boundary Effects The following experiment was designed for testing the difference between an R C M and its approximation arising from the spatial edge problem. In order to compare the "true" closed-loop transfer matrices to their R C M approximations, the following definitions are required: f {z) = F T (z)F%, (3.27) f {z) = F T {*)F% (3.28) yd ud where the Fourier matrices F m n yd ud and F are defined in (3.8), the integer n =m/n. Then m n a define {t (vf,z),--- ,tyd{v ,z)} = diag(f (z)), (3.29) { L K . 4 - , U C ^ ) } = DIAGCTU*)), (3.30) y yd where v and y m yd are spatial frequencies with engineering units in (3.18) corresponding to the normalized spatial frequency £l (j) in (3.17) and Q {j) in (3.15) respectively, 'dmp(T)' y u denotes the operation of extracting the diagonal elements of T , 'DIAG(T)' denotes the operation of getting the following elements of a rectangular matrix T £ C f (z)), ud n x m (such as 1 D I A G ( T ) = {T(l, 1), • • • ,T(k, k),T(k + l,k + l+m-n),--- ,T(n,m)}, where k = n/2 for even n or k — (n + l)/2 for odd n. If T (z) yd in (3.4) are true RCMs, then t (Vj,z) in (3.6) and T (z) ud and t d(Vj,z) would contain all of the frequency U yd information. (3.31) However, as mentioned above, the edge effects of the spatial interaction matrix Go in (2.4) and (2.6) lead to an approximation. Consider the closed-loop response of the system to a step disturbance d(z) £ C v m in Figure 3.1 of a single spatial frequency (i.e., the disturbance is a step function in the dynamical domain and a unit-magnitude y jd sine wave with a spatial frequency v in the spatial domain), d(z) = F d(z) y jd d {z) = { ^ ' 3 0, m f o r ^ It is a straightforward extension for the case T S c x = ^ for 1 < j < m,j mXn , y£j d , such as G{z) in (2.2). £ C , and m (3.32) Chapter 3. Two-Dimensional Frequencies in CD Control 40 the true response is given by y(z) = f {z)-d(z), (3.33) u(z) = f (z)-d(z), (3.34) yd ud and the R C M approximation response is obtained from (3.27) and (3.28), 1 < J < m,j ~ u K z ) .'-v = i £ Figure 3.4 compares the approximations \y(u^,z)\ in (3.35), . ^ = £ \u(Vj,z)\ ^jd . ( ,36) in (3.36) with the true responses \y(z)\ in (3.33) and \u(z)\ in (3.34) respectively at two spatial frequencies when they reach steady-state (i.e., z —> 1). The plots show a dominant response at the v = Vj frequency, the transfer matrices T (z) yd in (3.6) and T (z) ud in (3.4) will be approximated as a SCM and an R C M respectively, enabling the use of SISO transfer functions t (Vj,z) yd and t (Vj,z) ud in (3.29), (3.30). Please note that the spatial frequency of the disturbance in Figure 3.4f is very close to the spatial cut-off frequency v B in (3.23). It is clear that there is no big difference in the responses as shown in Figure 3.4d and Figure 3.4e for such a disturbance. 3.4.2 Input and Output Dimensions Related by a Non-Integer Factor In practical C D control it is generally not true that the number of measurements m is an integer multiple of the number of actuators n, in which case the process model matrix G in (2.2) will not be a rectangular circulant matrix by Definition 1. However a 0 technique analogous to that used to change the sampling frequency of a continuous signal in multirate signal processing [86] may be used to obtain a useful result by upsampling and downsampling in turn. First, create a new spatial interaction matrix G\ € R i m x n where mi = m • n by upsampling the continuous response g(x) in (2.4) as in (2.5) at the Chapter 3. Two-Dimensional Frequencies in CD 41 Control 0.08 (d) (a) 0.06 i^O.5 0.04 0.02 0 1 2 3 (e) I . ^ 1 2 3 4 •i 4 -, i-^ 1 2 3 1 2 I . 0.5 3 4 4 spatial frequency v [cycles/metre] spatial frequency v [cycles/metre] Figure 3.4: Response of the closed-loop system to steady-state disturbances at spatial frequencies v = 1.42 cycles/metre and v = 2.84 cycles/metre, illustration of modelled process outputs \y\ (dotted) in (3.33) and approximated by \y\ in (3.35), and steady-state actuator set-points \u\ (dotted) in (3.34) and approximated by \u\ in (3.36). corresponding high resolution. Then, create a second matrix 1 0 ••• 0 G 2 = 0 0 ••• 0 ••• 0 0 ••• 0 ••• 0 0 ••• 0 0 0 ••• 0 1 0 ••• 0 0 0 ••• 0 0 0 ••• 0 (3.37) 1 0 for downsampling from the high resolution m i to the true measurement resolution m . Then, the model matrix in (2.4) is given by, (3.38) Chapter 3. Two-Dimensional 42 Frequencies in CD Control The downsampling matrix G in (3.37) is an R C M by Definition 1, and we will approx2 imate the upsampled model matrix G\ G R miX ™ with an R C M through the use of periodic boundary conditions (as discussed in Section 3.4.1 above). However, G = G • Gi itself is 0 not an R C M due to the fact that the dimension of G G R m x m 2 i and d 2 GR ^ m xn do not satisfy property 2 of the R C M . From the RCM's property 4, both G = F G F" 2 m 2 1 and G\ = F GiF" are each a mi matrix of diagonal blocks as illustrated in (3.9), and further, as the column vectors of G\ satisfy the condition u B < of Theorem 1, ~ J ^0, ^3)\ 1=0, f o r i = j + fc . , otherwise G (3-39) where (similar to (3.25)), when j = 1, • • • , floor(n/2) +1 (3.40) when j = floor(n/2) + 2, • • • , n Therefore, G — G G\ £ C 0 2 m x n has the following characteristic: n c -\ ) ' Go(«,j)i -> for z = j + fc , . , otherwise (3.41) + where k is given in (3.25), but in this case m is not restricted to be an integer multiple of n. Comparing (3.41) with (3.24), G G C m x n 0 with a non-integer n = m/n will have a a properly analogous property to that of the R C M N in Theorem 1. 3.5 Two-Dimensional Frequency Analysis In this section, we will first analyze the open-loop system G(z) in (2.1) in the twodimensional frequency domain, then check the structured uncertainty A (z) G2 in (2.8) in the two-dimensional frequency domain, and finally examine the closed-loop system behavior of the control system in Figure 3.1 in terms of the two-dimensional frequency domain with t i(i'j,z) yc 3.5.1 in (3.29) and t d(i'j,z) u in (3.30). Open-Loop Transfer Functions in the Two-Dimensional Frequency Domain From the previous section, we know that the open-loop plant model G(z) in (2.1) of C D processes is an approximate R C M . Its corresponding representation approximation in the Chapter 3. Two-Dimensional Frequencies in CD 43 Control two-dimensional frequency domain can then be obtained by { S K , * ) , • • • ,9M,z)} = (3.42) BIAG(F G(z)F?), m similar to the approach used in (3.27)-(3.30). Figure 3.5 shows the nominal model's representation g(i/j, z) in (3.42) in the twodimensional frequency domain. The spatial Nyquist frequency of the inputs cycles/metre in (3.19) (X frequency = 0.033 cycles/second (here the sampling time T U>N u = O.llmetre = 4.54 in (3.19) in this case) and the dynamical Nyquist s = 15 seconds). Figure 3.5 illustrates that the process model g(u ,z) in (2.1), (2.2) and (3.42) has a large gain y for low spatial and dynamical frequencies, and small gain for high spatial and dynamical frequencies. This high-frequency gain roll-off is a familiar feature of physical processes and limits the controllability of disturbances d(z) in (2.1) to the low spatial and dynamical frequencies [20, 29, 82]. The ill-conditioning of Go in (2.2) can be interpreted using (3.26) in Corollary 1 as gain roll-off for high spatial frequencies [20, 29]. ° - T d y n a m i c a l Nyquist frequency 0 6 spatial Nyquist frequency 0.05- d y n a m i c a l frequency co [cycles/second] spatial frequency v [cycles/metre] Figure 3.5: The two-dimensional frequency approximation g(v?, z) in (3.42) of the nominal plant model G{z) in (2.1). Chapter 3. Two-Dimensional 3.5.2 Frequencies in CD Control 44 Structured Uncertainty in the Two-Dimensional Frequency Domain If the nominal spatial response matrix G in (2.4) is an R C M , then, all the possible models 0 G (z) p in (2.8) are RCMs. This is because changing parameters such as 7 , a in (2.3) and multiplication of the scalar function h{z) in (2.2) will not change the RCM's property. Based on one of the RCM's properties, i.e., two RCM's sum or subtraction is an R C M , all the additive structured uncertainties A (z) in (2.8) are RCMs. Then, from Theorem G2 1, like f (z) ud in (3.28), A (z) G2 A G2(i e C m x n can be "diagonalized" as follows, , ( z ) = F A (z)F» 3) m G2 J ^ °' *= I =0, otherwise, j + (3.43) where k = 0 for j = 1, • • • , n/2 and k — m — n for j = n/2 + 1, • • • , n for even n or k ~ 0 for j = 1, • • • , (n + l)/2 and k = m - n for j = (n + l)/2 + 1, • • • , n for odd n. The structured uncertainty A (z) has a two-dimensional frequency representation 8(Uj,z) G2 and its singular values can be obtained by {d(vlz),... ,8(v:,z)} a (A (z)) } G2 = DIAG(F A (z)F»), (3.44) = \6W,z)\, (3.45) m G2 j = l,---,n. For instance, if the nominal parameters in (2.3) of the nominal plant model G(z) in (2.1) have uncertainty limits 5i = 0.25 and 8 = 0.6 in (2.13), the additive structured 2 model uncertainty A (2;) can be calculated from (2.8). Figure 3.6 shows the nominal G2 steady-state spatial responses to a single actuator and the responses with parameter uncertainties in the spatial domain. Figure 3.7 illustrates the corresponding perturbations 5(Vj,z) in (3.44) in the spatial frequency domain. Note the figures only show the structured uncertainty at the dynamical steady-state case (OJ = 0) in the spatial frequency domain for clear display. In Figures 3.6a-e, only one parameter has uncertainty, while Figure 3.6f shows a class of response shapes combining uncertainty in all parameters. Similarly, Figures 3.7a-e illustrate the effect of uncertainty in one parameter as indicated in order to separate each parameter's contributions to the total structured model uncertainty. Figure 3.7f shows a class of structured model uncertainties which include all the parameters uncertainty in the spatial frequency domain. The dashed line may be considered as the maximum uncertainty limits across the spatial frequencies at steady-state. Figure 3.8 shows a realization of an additive structured perturbation £(i/", z) in the twodimensional frequency domain. Here only one G (z) is used which has different parameters p Chapter 3. Two-Dimensional Frequencies in CD Control x 10" 45 x 10" x 10" -20 x 10 ' x 10"' -20 -20 0 x 10 0 -20 20 CD locations [measurement points] 0 -20 20 CD locations [measurement points] 20 CD locations [measurement points] Figure 3.6: The nominal spatial response g in (2.4) (solid+dotted line) and its corresponding one with uncertainty g in (2.10) (solid line): (a) the gain parameter uncertainty 5 — —0.25 in (2.13); (b) the attenuation parameter uncertainty 6 = —0.25 in (2.13); (c) the width parameter uncertainty 6c = -0.25 in (2.13); (d) the divergence parameter uncertainty 8j3 — 0.25 in (2.13); (e) the alignment parameter uncertainty 6 = —3 in (2.13); (f) all the parameters uncertainty limits 8i € [—0.25,0.25] and 5 £ [—0.6,0.6] in (2.13), (note that only some extreme cases are shown for clear display). t ip 1 a C 2 than the nominal ones: 7 - 0.0084, a = 1.56, £ = 13, 8 = 0.13, a = 0.3679, q = 2, and d = 5.5 + 5(i - 1) with i = !,••• ,64 in (2.2)-(2.3). Generally speaking, it is always possible to diagonalize AG (Z) 2 in (2.8) at each dynam- ical frequency u using the singular value decomposition (SVD). However, it is difficult to directly relate singular vectors to physical variables, such as spatial frequencies. When the R C M structure AQ2(Z) in (2.8) is exploited, then its singular vectors can be given by the harmonic functions of the spatial variable [24, 25, 82, 83]. Figure 3.9 shows A (^)'s G2 representation in the input singular vectors vs dynamical frequency space. From Corollary 1 in (3.26), we know that the maximum singular value in Figure 3.9 is the same as the maximum gain in Figure 3.8. A casual analysis of Figure 3.7f indicates that uncertainty in the selected model param- Chapter 3. Two-Dimensional Frequencies in CD Control 46 x 10" (b) / \ * - 5 = -0.25 II a \\\| 1IAi = / il a 5 5 = 0.25 -0.05 (d) 5 = 0.05 0.005 0 2 v [cycles/metre] 4 0 2 2 4 4 v [cycles/metre] v [cycles/metre] Figure 3.7: The additive structured uncertainty \5(Vj,e \ in (3.44) for u — 0 in the spatial frequency domain due to: (a) the gain parameter uncertainty £ £ [—0.25,0.25] in (2.13); (b) the attenuation parameter uncertainty S £ [-0.25,0] in (2.13); (c) the width parameter uncertainty 5% £ [—0.25,0] in (2.13); (d) the divergence parameter uncertainty 6/3 £ [0,0.25] in (2.13); (e) the alignment parameter uncertainty S £ [—3,3] in (2.13); (f) all the parameters uncertainty limits Si € [—0.25,0.25] and S £ [—0.6,0.6] in (2.13). (note that the vertical axes have different scales for clear display) l2nuJ 7 a c 2 eters has its strongest effect at the mid-range spatial frequencies. While it is not unusual to find systems with low uncertainty at low frequencies, these plots should not be interpreted to indicate that high spatial frequencies are easily controllable. On the contrary, the relatively low model uncertainty simply indicates that we are relatively certain that the process model gain is small at high spatial frequencies. This important characteristic is very useful for robust C D control design. It helps to remove unnecessary conservativeness of the robust CD controller based on the maximum singular values of unstructured model uncertainty [23, 77]. Chapter 3. Two-Dimensional Frequencies in CD Control 17 v [cycles/metre] Figure 3.8: One structured model uncertainty in the two-dimensional frequency domain. 3.5.3 Closed-Loop Transfer Functions in the Two-Dimensional Frequency Domain Figures 3.10 and 3.11 show the two-dimensional frequency representations t d(Vj,z) y (3.29) and t d(Vj, u in z) in (3.30) of the closed-loop transfer matrices T (z) and T (z) in yd ud (3.6) and (3.4) respectively. Here again the nominal model G(z) in (2.1) of Section 2.1 is used, and horizons H p = 8 and H u = 1, and weights Qi = 41, Q = 10~ /, and 4 2 Qz = 0.002/ in (2.26). The flat surface in Figure 3.10 will be used for the definition of two-dimensional frequency bandwidth in Chapter 5. The low-gain part in Figure 3.10 means that the disturbances can be rejected by this controller in the low spatial and dynamical frequency range. If the disturbances belong to high spatial frequencies (such as over 3 cycles/metre) and high dynamical frequencies (such as over 10~ 2 cycles/second), they will not be rejected because of high gains in the corresponding frequency range. Note that if the disturbances belong to high dynamical frequencies (such as 0.05 cycles/second, those disturbances will be amplified because of the waterbed effects [75]. From Figure 3.11, the controller responds more aggressively to the lower spatial frequency disturbances (such as 0 ~ 3 cycles/metre) than to the higher spatial frequency disturbances (such as 3 ~ 4 cycles/metre). In other words, the controller tries to remove the lower spatial Chapter 3. Two-Dimensional Frequencies in CD Control 48 Figure 3.9: The structured model uncertainty in Figure 3.8 in the input singular vectors and dynamical frequency domain. frequency disturbances and leave the higher spatial frequency disturbances untouched. 3.6 Summary In this chapter, a rectangular circulant matrix (RCM) was defined and its properties were discussed in detail. If the column vectors of an R C M are bandlimited to the spatial Nyquist frequency of the inputs u^, its singular values can be easily obtained by pre- and post-multiplication of the complex Fourier matrices. The R C M structure can be preserved for the closed-loop transfer matrices in a linear feedback system if the open-loop plant and the controller are each an R C M and both of their column vectors are bandlimited to v^. The process model G(z), the structured uncertainty Ac2([z), and the closed-loop transfer functions T (z) and T {z) were approximated as RCMs and then decoupled into the yd ud two-dimensional frequency domain. This allows the large-scale MIMO closed-loop transfer matrices to be analyzed in terms of a family SISO transfer functions, one for each spatial frequency. The two-dimensional analysis method provides a better insight into the ill-conditioning characteristic of the process model and the role of the structured model uncertainty. It also provides theoretical basis for analyzing the closed-loop system's performance and robust stability, which will be explored in Chapter 5. Chapter 3. Two-Dimensional Frequencies in CD Control 49 Figure 3.10: The two-dimensional frequency approximation t (v}f, z) in (3.29) of the sensitivity function T (z) in (3.6) with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H = 1, Q = 41, Q = 10~ I, and Q = 0.002/ in (2.26). yd yd p 6CK U x 2 4 3 "" Figure 3.11: The two-dimensional frequency approximation t^iy^z) in (3.30) of the closed-loop transfer matrix T^z) in (3.4) with the nominal process model G(z) in (2.1) and the controller parameters H = 8, H = 1, Qi = 41, Q = 10^ /, and Q = 0.002/ in (2.26). 4 p u 2 3 Chapter 4 Model Reduction The problem statement of Section 2.3 requires that an effective and efficient model reduction method and a reduced-model based MPC be designed for a large-scale CD process to satisfy reduction and implementation specifications. The motivation for using wavelet packets for model reduction comes from well-known signal compression theory and previous application of wavelets in paper machine data analysis [61]. The reduced model using the wavelets method maintains the band diagonal and sparse structure features of the original CD process model, which are required by the CD-MPC QP solver [3, 7]. This chapter includes a summary of the main properties of wavelets, wavelet filters, and wavelet packets which may be found in [11, 85]. In Section 4.2, it is demonstrated that a given signal can be decomposed into and reconstructed from wavelet coefficients using wavelet matrices. In Section 4.3, a modified wavelet packets method is designed for model reduction and an example is shown that how efficiently reduce large-scale signals or models. In Sections 4.4, some practical issues such as distortion correction, the relationship between the spatial controllability and the reduced dimensions are discussed and summarized. In Section 4.5, implementation of MPC in the wavelet domain is investigated. In Section 4.6, the method is extended to multiple-array systems. Finally, Section 4.7 provides two examples for demonstrating the method's effectiveness and efficiency. 4.1 Signal Decomposition and Reconstruction Using Discrete Wavelet Transform Wavelet transform is an alternative to the Fourier transform. Instead of transforming a pure "time" description into a pure "frequency" description, wavelet transform provides a good compromise — a "time-frequency" description. Note that "space-spatial frequency" replaces "time-frequency" in CD processes. A given function f(t) € L (R) can be decomposed as [11]: 2 k = k j<J < fit), h,k > hAt) +E E < /(*)>> M*) k k j<J 50 (-) 4 2 Chapter 4. Model Reduction where Cj(k), dj(k) 51 are defined as the scaling coefficients and the wavelet coefficients respec- tively, the symbol ' < , > ' represents the inner product, the scale level J and the shifting variable k are integers, (j)j k t and tpj^ are the scaling functions and the wavelet functions, which are obtained by dilation and translation of the scaling and mother wavelet functions <j> and ip respectively. Figure 4.1 shows a scaling function and a wavelet function of the wavelets-symlets [57]. (a) _ 0 . 5 I 1 0 1 1 1 1 2 3 1 4 i 5 6 7 (b) 0 1 2 Figure 4.1: (a)The scaling function 3 4 5 6 7 of sym4; (b)The wavelet function <f>(t) ip(t) of syn%4- In most applications, one never deals directly with the scaling functions or the wavelet functions. Only the low-pass filter coefficients ho(n) and the high-pass filter coefficients hi{n) need to be considered [11, 85] through the following relationships: cj)(t) = ^ n h (n)V2(f)(2t 0 - n), ^(t) = 5 ^ / i ( n ) ^ ( 2 t - n ) . 1 (4.3) (4.4) n Finally, the coefficients cj(k) and dj(k) in (4.1) are obtained from the following equations: Cj+i(k) = Y h {l Q } - 2k)Cj{k), (4.5) Chapter 4. Model Reduction 52 dj+i(k) = ^2hi(lf (4.6) - 2k)Cj(k), where j > 0 (for j = 0, c {k) is the same as the original signal /(£)), and // represents the length of the filter. The above equations (4.5) and (4.6) can be easily implemented by the two-band analysis tree in Figure 4.2, where the symbol '| 2' (downsampling) throws away every second data point. 0 h,(-n) 12 ho(-n) 12 h.(-n) 1 2 •*j+2 ho(-n) 12 * j+2 C:'j+l C Figure 4.2: Two-stage two-band analysis tree. Similarly, the original signal f(t) can be perfectly reconstructed from the coefficients through the two-band synthesis tree in Figure 4.3, where the symbol '| 2' (upsampling) inserts one zero between every two samples in order to get the same length of the original signal at scale level 0, and g (n) = h {n), gi(n) = hi(n). Figure 4.4 and Figure 4.5 0 •*j+2 "j+2 f 2 ! 2 0 f 2 &(n) t 2 go(n) &(n) go(n) Figure 4.3: Two-stage two-band synthesis tree. illustrate the filters and their frequency responses. For orthogonal wavelets, the relationship between the low-pass filter coefficients and the high-pass filter coefficients is given by, /ii(fc) = (-l)* /io(//-fc), +1 (4-7) where k = 0,1,2, • • • // -1, // is the length of the filter. In addition, the low-pass filter Chapter 4. Model Reduction 53 Lowpass filter for decomposition Highpass filter for decomposition Lowpass filter for reconstruction Highpass filter for reconstruction Figure 4.4: The decomposition and reconstruction filters of sym4- ho(n) and the high-pass filter hi(n) are double-shift orthogonal, J2h (n)h (n-2k) = 6(h), (4.8) Y^ho(n)h {n-2k) = 0, (4.9) ^2h {n)h (n-2k) = 5(k), 0 0 v l l (4.10) where the Kronecker delta function 6(k) is 8(k) = l I 0, h 4.2 k ° otherwise = (4.11) Signal Decomposition and Reconstruction Using the Wavelet Matrices From Section 4.1, it is known that signal decomposition and reconstruction can be implemented by low-pass and high-pass filters. It is convenient to express the relationship between the signal and the coefficients (scaling and wavelet) in matrix form. One-stage 54 Chapter 4. Model Reduction The l o w p a s s filter for d e c o m p o s i t i o n The h i g h p a s s filter for d e c o m p o s i t i o n The l o w p a s s filter for reconstruction The h i g h p a s s filter for reconstruction 0.1 0.2 0.3 0.4 N o r m a l i z e d frequency 0.5 0.1 0.2 0.3 0.4 N o r m a l i z e d frequency 0.5 Figure 4.5: Spectra of the decomposition and reconstruction filters for sym4. two-band analysis can be expressed as r 0 0 r r ri r 0 0 r r n 0 0 0 0 0 0 0 0 • • r 0. 0 • • 0 0 3 Cj+l(l) 2 3 0 2 0 0 0 • • 0 0 • • 0 r ' •• 0 0 0 0 0 0 0 0 0 0 0 r ri r 0 r r c (2) j+l c i{m) j+ 0 0 3 2 3 0 2 (4.12) <Wi) Si So 0 0 <W2) S3 s Sl So 0 0 0 0 s s 0 0 0 0 0 0 2 3 2 Si 0 • • 0 0 • • 0 So • • 0 0 0 0 0 0 0 0 0 0 0 0 0 • • s 0 • • 0 s Si So s s2 d i(m) _ j+ 0 0 3 2 0 3 Chapter 4. Model Reduction where m = floor( ~ ), N+l f 1 55 the filter's length is assumed to // = 4, and r 1 2 n h\(ri), the large matrix in (4.12) is called the wavelet matrix W . T = h (n),s 0 n — For orthogonal wavelets, signal reconstruction can be easily implemented using the following formula, c,-+i(l) c (2) j+l c,-(l) c,-(2) = W c i(m) j+ (4.13) d (l) j+l Ci(N) d {2) j+l dj+i(m) It is worth noting that the above equations (4.12) and (4.13) only apply to one-stage analysis and synthesis. As for multi-resolution analysis and synthesis, the dimensions of wavelet matrices should be adjusted for different scale levels. For example, the original signal's length (at scale level 0) is 1000 and the filter length / / is 4, thus at scale level 1, the length of the scaling coefficients or the wavelet coefficients is 501, then the wavelet matrix between the scaling (or wavelet) coefficients and the original signal Wi is 501 x 1000, but at scale level 2, the wavelet matrix W between level 1 and level 2 is changed to 252 x 501. 2 Therefore, the wavelet matrix between the original signal and the coefficients at level 2 is W eR 3 4.3 4.3.1 2 5 2 x l 0 0 0 = W xW . 2 1 Model Reduction Modified Wavelet Packet Analysis Not all components of the measurement profiles (basis weight, moisture, and caliper) of paper machine CD processes can be controlled due to the limited spatial bandwidth [20]. These profiles can be split into lower and higher spectral parts [50]. The higher spectral part is considered to be uncontrollable. In CD control systems, there are two important spatial frequencies: one is the Spatial Nyquist frequency of the inputs v% in (3.19), another is the spatial cut-off frequency v which depends on the spatial response of a single actuator c and is usually quite smaller than the above v^. In [78], the authors provided some practical rules for estimating this frequency v in paper machine C D processes. If g x c ma represents the maximum gain of the spatial frequency response of a single C D actuator and g Vc 1 floor(X) is a Matlab function which represents the operation of getting the integer part of X. Chapter 4. Model Reduction 56 represents the gain at the cut-off spatial frequency, then (4.14) 9v = r • 9r c where r is a scalar number between 0.05 ~ 0.2. This cut-off spatial frequency v is smaller c than the spatial bandwidth v B in (3.23). From Theorem 1, Corollary 1 in Section 3.3 of Chapter 3 and (D.ll), it is known that this cut-off frequency is closely related to the cut-off singular value ji, 1V (4.15) C yJmjrC where m and n are dimensions of the outputs and inputs respectively. This cut-off frequency indicates how much dimensional reduction can be achieved without losing any controllable components. If the cut-off frequency v is known, the lower frequency part can be obtained by c decomposing the original signal. However, this part can not normally be exactly represented by the scaling coefficients obtained by the above standard discrete wavelet analysis tree. This is where wavelet packet analysis becomes useful. A two-stage standard wavelet packet analysis tree is shown in Figure 4.6, where S is the original signal, cal and cdl are h,(-n) h,(-n) 12 cdo^ ho(-n) 12 •cda. h,(-n) 12 cacl. h (-n) 12 -caa 12 12 0 0 Figure 4.6: Two-stage wavelet packet analysis tree. scale level l's coefficients, and caa2, cad2, cda2, and cdd2 are scale level 2's coefficients. Please note that those coefficients represent different frequency bandwidths but are not necessarily sorted from low to high [57]. Please see Appendix B for the correct frequency sequence. Here a modified wavelet packet method is introduced as an attractive alternative to Chapter 4. Model Reduction 57 The length of the original signal is assumed as N. standard wavelet packet analysis. First, a proper scale level L is selected using a standard wavelet analysis tree. At this level L , the approximation profile generated from the scaling coefficients with k = CL(1C) 1,2, • • • ,n,L has a spatial frequency bandwidth VB in (3.23) which is just smaller than the cut-off frequency v in (4.14). This can be justified by checking the bandwidth of the c approximation profile at scale level L — l, where the bandwidth should be bigger than v . c Second, only the wavelet coefficients di,(k) at scale level L are continuously decomposed to another proper scale level P (or several scale levels). At this P level, the spectral part contributed by the scaling coefficients cp(k) with k = 1,2, • • • ,np should be close to the spectral difference between VB and v . Finally, the reduced dimension signal includes the c above scaling coefficients CL and Cp, and its length is UL + np - C N. UL and rip can be computed by (4.16) fc=i where X is the scale level number. Compared to standard wavelet packet analysis, the modified wavelet packet method only decomposes the interesting scale levels and wavelet coefficients. An example of the modified wavelet packet analysis tree is shown in Figure 4.7, where n L = 134 and rip = 38, ho and hi are the low-pass and high-pass filter coefficients, Cj(k) and dj(k) are the scaling and the wavelet coefficients at scale level j. The original signal (1024 data numbers) is decomposed to 172 scaling and wavelet coefficients which come from scale level 3 and scale level 5. 4.3.2 Model Reduction Using the Modified Wavelet Packets Method [21] If y (z) and y(z) represent the controllable part of the measurement profile y(z) in (2.1) r in the wavelet domain and spatial domain respectively, we have the following relationship: y (z) (4.17) r = W y(z), y(zj = < y ( « ) = W^W y(z), y r (4.18) tf where W is a p x m (p <C m) matrix and y(z) is also the controllable approximation y profile of y(z). Generally speaking, WyW y works like a low-pass filter, for spatial frequency components of y{z) below v c for spatial frequency components of y{z) over v c (4.19) Chapter 4. Model Reduction 58 S i g n a l S=c(0), l e n g t h Vs length=8 \ cM d/k) Cjft) h^(-n)/ Level 1 \ k=515 * d,(k) \K(-n) Level 2 k=261 \ \ Level 3 \ k = @ ) ho(-nj/ \h,(-n) c (k) 4 h^-n)/ Level 4 k=70 d (k) 4 \ h ( - n ) i Level 5 k=C38 d (k) 5 Figure 4.7: An example of the modified wavelet packet analysis. where v is the cut-off spatial frequency in (4.14). This is due to the fact that different c scale wavelet coefficients have different frequency bandwidth [85] and the wavelet matrix W is composed of several scale level wavelet coefficients. y Similarly, the actuator profile can be decomposed and reconstructed by u {z) r u{z) = W u(z), (4.20) u = W^u (z) = WjW u{z), r (4.21) u where W is a q x n matrix and u(z) is the approximation profile of u(z) in (2.1). u W£W U removes high spatial components from the actuator profile u(z). Therefore, the reduced model can be obtained from (2.1), (4.17) and (4.21) by y (z) r = W G h(z)u(z) +Wd (4.22) W G h(z)u(z) +Wd (4.23) G {z)u {z) r y y = G (z) r dr(z) where G {z) G C r pxq 0 0 r r y y + d (z), (4.24) = W G h(z)WZ, (4.25) = Wyd(z), (4.26) y 0 is the reduced model in discrete wavelet domain, d (z) G C is a vector r Chapter 4. Model Reduction 59 of the disturbance in the wavelet domain. We also.can construct the approximation model G(z) from the reduced model G (z) in (4.25) as follows: r G(z) = W^G {z)W r = W W G(z)W^W . y y u Figure 4.8 shows the singular values of WfW (W e R y M 8 2 x 1 1 0 0 (4.27) T u y 8 6 x U 0 0 ) and W^W U (W u e ) for the model reduction example in Section 4.7.1. The discontinuity in Figures 4.8a and 4.8b is related to the cut-off singular 7; in (4.15). Note that the singular values and the eigenvalues of W W y y or W^W U are same. (a) 400 600 800 Singular vector directions 40 60 Singular vector directions Figure 4.8: A n example for illustrating the singular values of W^Wy (a) and W^W U (b). G(z) is useful for evaluating the distortions during the model reduction process, which will be discussed in Section 4.4. Figure 4.9 illustrates the construction of the signals that define the model reduction. 4.3.3 Quantification of the Model Reduction Requirements In general, a-system's controllability is discussed in terms of the output signal y(z) (e.g., the cut-off frequency in (4.14)). For C D processes, this cut-off frequency is related to the system's singular values [24, 78]. If a singular value is large, then the corresponding singular vector is easily controllable. Small singular values correspond to non-robustly controllable modes. Chapter 4. Model Reduction y 60 G W W GW W T T y y u u u Figure 4.9: Illustration of the construction of the signals that define the model reduction, y in (2.1) and y in (4.18). The singular value decomposition of the steady-state process model G ss G h(l) 0 (i.e. G ( l ) = in (2.2)) is: = U • E •V . G The input singular vectors of G (4.28) H ss are given by the columns of V in (4.28) denoted by Vj. ss Define the set of controllable input singular vectors as, IL := [Vj : oAGss) > li], (4.29) then the set of non-robustly controllable singular vectors will be defined by, U h := {VJ : CjiG,,) < -yi}, (4.30) where Oj is the singular value corresponding to the singular input vector Vj. The selection for -y depends on the process spatial response. Detailed information can t be found in [78] for a possible way to select its value. It is necessary to quantify the reduction performance for the subsequent reduction algorithm development in Chapter 4. 1. Dimensionality reduction. It is required that p. < m and q < n, The numbers m,n,p, and q are defined in (2.1) and (4.25). 2. Preservation of controllable components. The reduced model is required to preserve the controllable components of the system. This may be quantified by the following equation, ^ \ , ~ ° ^ <w s V h u eU s t (4.31) II s ||2 where y , y , and u are the reconstructed output profile in (4.18), the output profile U s s s Chapter 4. Model Reduction 61 in (4.18), and the input profile in (4.21) at steady-state respectively , 11/ is the set of controllable input singular vectors defined in (4.29), and wi is the tolerance which can be used as the distortion index. A smaller value of w is preferred as it t corresponds to a smaller distortion of the controllable components. 3. Elimination of non-robustly controllable components. The reduced model is required to eliminate the non-robustly controllable components of the system. This may be quantified by the following equation, u eU j j - ^ < w, II U II2 h s (4.32) h s where n „ is the set of non-robustly controllable singular vectors defined in (4.30) and Wh is the tolerance which can be used as the index of eliminating the nonrobustly controllable components. A smaller Wh means a better elimination of the non-robustly controllable components. 4. Reduction of the condition number. It is required that the reduced model's condition number T ( G ) be small. Well-conditioned systems are more easily conr trolled than ill-conditioned ones, because the non-robustly controllable directions do not exist in the well-conditioned systems. 5. Sparsity of the reduced model. For its use in the quadratic program solver the reduced matrix G G C pxq r is required to be a sparse, banded matrix. Sparseness of a matrix may be quantified by the ratio w of nonzero elements to the total number s of elements, 77 w, = — p-q <l, (4.33) where n is the number of nonzero elements in the reduced model G . z 4.4 4.4.1 r Practical Issues in Reduction Implementation Distortion Correction If the signal is of infinite length, the wavelet matrix is also infinite and no distortions exist. Because the signal is always finite and there are convolution operations between the signal and the wavelet filters during the decomposition and the reconstruction, there will be distortions in the first and last few rows of the wavelet matrix. The number of seriously distorted rows in the wavelet matrix depends on the length of the filter // and Chapter 4. Model Reduction 62 the scale level number X: /N />i n dr = floor{— + "V") - / («). Y. fe=i (-) „,|oor/-/V. \-^v If — 1. 4 34 where rirf is the number of the seriously distorted rows, and N is the original dimension. The following example illustrates how to remove these distorted rows in the reduced models. The model G, the output y(z) and the input u{z) in (2.1) are assumed to be 693 x 226,693 x 1 and 226 x 1 respectively; the filter's length l is 8. Decomposing the signal to scale level 1 is assumed. According to (4.16), the length of the reduced output y {z) is 350, that is, the distorted wavelet matrix W is 350 x 693. 4 rows have to be removed from W . Which rows need to be removed depends on what kind of wavelets is used. If sym4 (a member of the symlets family [57] shown in Figure 4.1) is used, the first row and the last three rows should be removed. Therefore, the wavelet matrix W is changed to 346 x 693. Similarly, the wavelet matrix W without the distorted rows is 113 x 226. Finally, the reduced model G without serious distortion is 346 x 113. r f r y y y u r 4.4.2 Choice of Wavelets Because the CD responses are usually symmetric except at the edges, it makes sense to choose wavelets • which are symmetric or nearly symmetric and have similar shapes as the CD response ones in order to effectively capture the features of the input and the output profiles. These wavelets include symlets (near symmetry), coiflets (near symmetry), and biorthogonal wavelets (symmetry) (refer to [57]).. However orthogonality is lost for biorthogonal wavelets. Therefore, symlets and coiflets are good wavelet families for CD process model reduction. One has to choose the length // of the low-pass filter and the high-pass filter in (4.7) used for model reduction. It depends on the desired sparsity of G in (4.25) and the allowable distortion between the original model G in (2.1) and the approximation model G in (4.27). If If is increased, G becomes more dense and the distortion between G and G will be decreased and vice versa. Our experience indicates that If between 8 to 16 is a good choice. As an example, consider the model of an industrial working paper machine which uses dilution actuators for controlling the basis weight. The original model at steady-state G in (2.2) is 693 x 226. The reduced model G (4.25) is 129 x 141. Figure 4.10 illustrates the relationship between the wavelet filter length and the distortion. Before we give detailed numbers about model reduction performance mentioned in Section 2.3, we may define the distortion ratios for preservation of controllable components and elimination of r r ss r Chapter 4. Model Reduction 63 non-robustly controllable components as E S^ -Vr( J^) (M Wi = = >lh — Wh(j) = E"i ( M ) > ^(G )' (4-35) s s (4.36) UG-G^-Vjh, -^n rr IC V 2^=min(j>,q) + l °jV^* ss) ' (4.38) WG-Vjh, where p, q are the dimensions of the reduced model G in (4.25) and n is the actuator number, G is in (4.27), the input singular vector Vj is the j column vector of V in (4.28). Smaller rji means smaller distortion introduced from the reduction process while smaller r] denotes larger elimination of non-robustly controllable components. These two ratios also apply for the original model. If this is the case, rji and r) are 0 and 1 respectively. The model reduction performance criteria for different wavelet filter lengths are shown in Table 4.1. From the numbers in the table, there is trade-off among the model reduction requirements, i.e., it is impossible to let all the indices such as wi in (4.31), Wh in (4.32) , and the sparsity index w in (4.33) be small. Note in Figures 4.10 and 4.11, and Table 4.1, only the values at steady-state (i.e., z = 1) are shown. r t h h h s Table 4.1: The model reduction performance criteria for different wavelet filters Wavelets Filter Condition Distortion Distortion Sparsity Model name length number T index r\i index t]h index w € 57.5 0 1 0.0236 eR sym4 8 3.36 0.1002 0.7076 0.1986 sym8 eR * 16 3.14 0.0449 0.3190 0.4210 coif2 12 3.34 eR * 0.0965 0.6781 0.2781 sym6 £ R129X141 12 3.24 0.0651 0.4922 0.2777 sym8 eR * 16 6.03 0.0457 0.2852 0.4413 s G G G G m x 1 4 1 rl 129 141 r2 m 141 150 155 r3 G rA G r5 For symlets and coiflets with same filter length (for example, sym6 and coif2, the filter length is 12), G in (4.27) generated by sym6 has smaller distortions than by coif2. Figure 4.11 shows the distortion difference due to different wavelets. The quantified indices may be found in Table 4.1. From this table, sym6 or sym8 and the reduced dimensions between 129 x 141 and 150 x 155 are good choices for this plant model. Chapter 4. Model Reduction 1.5 I 64 I I I ''''. '''., x> I I I (a) X sym4 1 3 1 ••] x cd \'-.. \ 0.5 - sym8 \ / N y / - •., ^ 0.05 - i i i i i i 0.1 0.15 0.2 0.25 0.3 0.35 '''•''<., 1 _ x 0.4 0.45 0.5 Normalized frequency 0.12 (b) 0.1 0)0.08 sym4: o(G )- c0.06 M CD sym&. o(G ) sym4: w ta sym4: iv sym8: w r2 So.04 ; 0.02 20 40 s y m g ; 60 80 w 100 120 140 160 180 200 220 Input singular direction v. Figure 4.10: (a) The spectra of the wavelet filters: sym4 (dotted) and sym8 (dashed); (b) Illustration of the distortion difference in the reduced model G in (4.25) due to different wavelet filter lengths. r 4.4.3 Spatial Controllability and Model Reduction Duncan and Bryant (1997) [20] proposed that the spatial bandwidth of an actuator array in web forming C D processes depends on single actuator's spatial response. Stewart et al. (2001) [78] further gave some practical rules about how to choose the cut-off spatial frequency v in (4.14) in C D processes. Only those spatial frequency components below the c cut-off spatial frequency v can be controlled. In fact, this frequency u is closely related c to the singular value c in (4.29), because the process model G(z) is an approximate R C M and its singular values can be obtained from (3.26) by multiplication of Fourier matrices. It is well known that most C D control systems are ill-conditioned or singular [50, 77]. In addition, due to model inaccuracies, some smaller singular values, which are related to high spatial frequencies, can not be reliably identified [27]. Therefore, the input directions corresponding to these smaller unreliable singular values can be discarded, that is, the input singular vectors belong to 11^ in (4.30) are not robustly-controllable. Generally speaking, choosing -ji in (4.29) between 0.1 ~ 0.2a(G ) is typically good ss for most C D processes [78], which implies that the reduced model's condition number Chapter 4. Model Reduction 0 65 50 100 150 200 Input singular vector v Figure 4.11: (a) The spectra of the wavelet filters: coif 2 (dotted) and sym6 (dashed); (b) Illustration of the distortion difference in the reduced model G in (4.25) due to different wavelets. r will be smaller than 10. We may use j to decide what numbers are good for the reduced dimensionsp and q of G in (4.25). Figure 4.12 shows the singular values and the distortion indices of the reduced model G e R . It is evident that the closed-loop bandwidth with this reduced model G is larger than one with the reduced model G 2 in Table 4.1 as G keeps 21 more input singular values than G . However, the condition number of G is larger than that of G . Therefore, there is trade-off between the reduced dimensions and the condition number (i.e., performance and robustness). t r 150x155 r5 r r5 r5 r2 r5 r2 4.4.4 M o d e l R e d u c t i o n D i s t o r t i o n Test We tested a large number of possible CD process models in order to verify the effectiveness of the model reduction method described above. The measurements rn in (2.1) ranged from 300 to 2000 by step 200, the actuator numbers n in (2.1) ranged from 40 to 300 by step 65, the attenuation parameter a in (2.3) was changed from 2 to 7 by step 1, the width parameter £ in (2.3) was changed from 2 to 10 by step 1, and the divergence parameter (3 in (2.3) was changed from 0 to 0.45 by step 0.15. The default controllable singular value 7* in (4.29) was set as Q.la(G ) for all these models. If the model's condition number is less ss Chapter 4. Model Reduction 66 0.12 0 50 100 150 Input sigular vectors v. Figure 4.12: Singular values and distortions in a reduced model from dilution actuator response. than 10, then the model will not be reduced. There are a total of 7704 effective models (needed to be reduced for the above variable parameters. We calculated the distortion ratio rji in (4.35) for each model using two different length wavelet filters. The statistical results for this ratio are shown in Figure 4.13. From the figure, the distortion ratios using the short wavelet filter are usually bigger than those using the long wavelet filters. 1 4.4.5 C o m p a r i s o n W i t h Other R e d u c t i o n M e t h o d s In industry, the measurement profile (high-resolution) is usually first filtered by an antialiasing filter (such as Hamming window) and then down-sampled (using an averaging method) to low-resolution profile (same resolution as the actuator profile), i.e. transfer the process model from a rectangle matrix to a square one. Here we calculate the distortions in the output profile using the proposed wavelet packets method and the commonly used averaging method. Figure 4.14 shows the distortions using these two methods for the previous same dilution actuator response model example in Section 4.4.2. From this figure, it is clear that the wavelet method is much better than the averaging method. These model are ill-conditioned because their condition numbers are more than 10. 67 Chapter 4. Model Reduction 0.25 0.25 distortion ratio r^ Figure 4.13: (a) The probability of the distortion ratio r)i in (4.35) using sym4; (b)The probability of the distortion ratio rji in (4.35) using sym8. Another advantage is that the condition number of the reduced model using the wavelet method (6.03 as shown in Table 4.1) is much smaller than that of the square model (which is 68 for this CD process). Compared with Gram polynomials in [50], this wavelet method can directly build relationship between the reduced dimension and the cut-off frequency. The wavelet method has also the localization capability, i.e. the wavelet coefficients are related to the spatial response locations. Other orthogonal bases including Gram polynomials and singular vectors (column vectors in U and V in (4.28)) lack this capability because those bases span over the entire sheet width. The sparsity feature of CD process models will be kept in the reduced model by using the wavelet method but lost if using other reduction methods. This method can be used for multiple-array systems, but the singular value decomposition method can not be used for such systems. Detailed information will be given in Section 4.6. Chapter 4. Model Reduction 68 v^=14.2 v =9.5 dotted: original solid:approximation u s i n g sym8 0.8 dashed+point: approximation using HamrAing+averaging <D 0.6 cd 0.4 0.2 11' ti my 6 8 10 12 Spatial f r e q u e n c y [cycles/metre] Figure 4.14: Spectra of the approximation output profiles using different methods 4.5 M P C Implementation in the Wavelet Domain 4.5.1 Constraint Handling The inequality constraints in (2.27 - 2.30) can be simply expressed by [22] AxAu(fc) < bi -Ci«(fc- 1), (4.39) where A bi, and Ci can be easily obtained after some algebra. Note that the above constraints are in the spatial domain, and need to be transferred to the wavelet domain. However, that transfer to the wavelet domain can be avoided simply by replacing u(k) in (4.39) with ii(fc).in (4.21), that is, inequality (4.39) is replaced by u AiW^Au (k) r <bi- CiW^u (k r - 1). (4.40) This inequality will guarantee that the constraints in the spatial domain will not be violated while the minimization of the upcoming objective function (4.41) is in the wavelet domain. The reason is that the optimal solution u (k) in the wavelet domain will be r Chapter 4. Model Reduction 69 projected back to the spatial domain by using the inverse wavelet transform (i.e., u(k) = W£u {k) in (4.21)) and then the solution u(k) is implemented in the spatial domain. Please note that it is u(k), not u(k), that is used for implementation. r 4.5.2 Objective Function in the Wavelet Domain The objective function for CD-MPC in the wavelet domain is: min 3 {k) r — min Au (k) Au (k r r Hu-1 ]T [Au {k+j) Q Au {k T r 2r + j) +u {k + j) Q u {k r T r 3r + j)) \ r (4.41) subject to (4.24) and (4.40), where y (k) in (4.24) and y (k) are the reduced output profile and the reduced desired profile in the wavelet domain respectively. y (k) can be gotten by W x y (k) analogous to y {k) in (4.17). Q , Q and Q have the same structures as Qi, Q , and Q in (2.26) respectively, but have the smaller dimensions. Due to the large size of the plant and simple dynamical characteristics, H is restricted up to 3, and H is chosen between 10 and 20, which correspond to typical settling times for paper machine CD processes. If we substitute y (k) with W y(k) from (4.17) and u (k) with W u{k) from (4.20) into (4.41), the objective function J(A;) is the same as the objective function 3{k) in (2.26) except for different weighting matrices. For the objective function 3 (k), the equivalent weighting matrices are WyQirW , W^Q W , and W^Q W . As discussed in Section 4.3.2, for controllable components of the output profile y below the cut-off frequency v in (4.14) (related to the cut-off singular value 7( in (4.29)), W W will not significantly change their magnitude. For uncontrollable components of the output profile y above the cut-off frequency v , W W will remove them. In other words, the equivalent weighting matrix W Q\ W is very close to the weight matrix Qi for the objective function 3(k) in (2.26) in terms of its effect on controllable components (refer to Figure 4.8a). Similarly, other two equivalent weighting matrices W^Q2 W and W^Q W are also very close to the weight matrices Q and Q respectively in the controllable subspace of G(z) (i.e., for controllable singular vectors Vj € Hi in (4.29)). Therefore, the optimal solution Au in (4.41) based on the reduced model is very close to the optimal solution Au in (2.26) based on the full model in the controllable subspace. Examples will be presented in Section 4.7. When the optimal solution of (4.41) Au (k) is calculated, then the reduced actuator profile u {k) is equal to Au {k)+u {k—1). Finally, the corresponding actuator profile u{k) in the spatial domain can be obtained by (4.21). Figure 4.15 shows the MPC controller r spr spr y sp T 2 ir 2r 3r 3 u p r y r u r r y 2r u 3r u c y c y T y y y y T 2 U 3r u 3 r r r r r Chapter 4. Model Reduction 70 implementation diagram in the wavelet domain. constraints sp vv, W„ • 'MPC vv, W T u — • plant " * u sensor Figure 4.15: MPC implementation in the wavelet domain Because the reduced model G (z) in (4.25) removes high spatial frequency components (i.e., small singular values), it inherently increases robustness of the MPC controller. The controller also will rarely result in a picketing phenomenon, which happens often in practical CD control systems, because the actuators are not permitted to operate in high spatial frequency range after transforming control actions from the wavelet domain to the spatial domain. The robustness and performance analysis method, which will be presented in Chapter 5, can be used for the wavelet-based MPC when W^QirWy, W^Q W , and W^Q2, W are considered as the equivalent weight matrices for the wavelet-based MPC. r 2r r 4.6 u u Extension to Multiple-Array Systems As mentioned in Section 1.2.5 and Section 2.1.2, it may be necessary to coordinate all actuator arrays for controlling all the properties in CD processes. It is known that the computational time of solving the QP problem such as (2.26) could exponentially increase with the variables of the problem (such as inputs and constraints) [87]. Due to several times more actuators, measurements, and constraints in a multiple-array system than in a single-array system, the on-line computational time for solving the QP problem is so long that the MPC controller may not be implementable in reality. Due to severe illconditioning both in single-array systems [20, 29, 50, 83] and in multiple-array systems [4], the model reduction method used for single-array systems can be used for multi-array systems too. The singular value decomposition method [87] cannot be used for multiple-array systems. The reason is that there is not such a constant spatial matrix in a multiple-array Chapter 4. 71 Model Reduction system model G(z) in (2.21) due to different sub-plants' dynamics (i.e., h (z) in (2.21)). However, the model reduction method based on the modified wavelet packets described above can be used for multiple-array systems. Here is the procedure: itj 1. We first decide the reduced dimensions for each controlled property and each actuator array. • For a particular property y we should keep the maximum controllable components in the reduced property y . For example, the moisture profile is controlled by slice lip actuators, steam boxes and rewet showers as shown in Figure 1.1 and the spatial cut-off frequencies for each of these subplants (a row of the multiple-array system model G(z) in (2.21)) are v , v , and i / in (4.14). Then, let the maximum of the three frequencies as the spatial bandwidth for the moisture. Based on the this bandwidth, using the method for the singlearray systems described in Section 4.3, step by step, we can get the reduced dimension for the moisture. iy Ti ci c2 c3 • For a particular actuator array Uj, we should keep the maximum control capacity after the model reduction. For example, a slice lip actuator array can affect the basis weight, the moisture, and the caliper profile. If the reduced dimensions for the slice lip actuators for each of these subplants (a column of the multiple-array system model G(z) in (2.21)) are q with j = 1,2,3, then, the reduced dimension for the slice lip actuators is the maximum of q j. Xj X 2. We then build the wavelet matrices W i, ••• , W ^ for each controlled property and W i, • • • ,W N for each actuator array. Here N and N are the numbers of the controlled properties and the actuator arrays in (2.21) respectively. y u U y y y U u 3. Finally, we build the reduced model for a multiple-array system as follows: Y (z) r G (z)U (z) r r + (4.42) D (z), r w, Vn{z) (4.43) Y (Z) T VrNyiz) W 0 ••• W yNy ••• yl VN {Z) V 0 G (z) G(z) r o ••• G>ll(z) 0 W,yN y • • w, uN u G iN {z) r u (4.44) GrN l{z) v G NyN {z) T u Chapter 4. Model Reduction = G (z) rij 72 WyiGijWujhi^z), i= 1, • • • , N , j = 0 u {z) rl U (z) y 1, • • • , N , Ui(z) (4.46) = r 0 •• • UrN {z) u drl(z) 'Wyy (4.4 u w _ u {z) uNu 0 " Nu _ dy(z) ' (4.47) D (z) r 0 •• drN {z) v _ d y(z) N where G , G , h y Uj, d N , N are from (2.21), y € C with i = 1, • • • , N , u E O' with j = 1, • • • , N . Pi tj ih u u y Ti u ri y u After the reduced model G (z) in (4.42) is constructed, the constraints and the objective function for the multiple-array system in the wavelet domain can be handled in a similar fashion as the single-array system in Section 4.5. r 4.7 Examples In this section, two examples are given for demonstrating the effectiveness and efficiency of the wavelet-based reduced-order CD-MPC,firstfor a single-array system and then for a multiple-array system. 4.7.1 Basis Weight Profile Controlled by Slice L i p A c t u a t o r s In order to illustrate the effectiveness of the above wavelet-based reduced-order MPC, a real paper machine CD process model for 39 gsm paper grade is used in the following example. This paper machine has 110 slice lip actuators and the scanning sensor has 1100 measurements per scan (30 seconds per scan). The shape of the actuator response is shown in Figure 4.16 and is very typical for slice lip actuators on a light weight paper grade. The frequency power spectrum of the response is also shown, the Spatial Nyquist frequencies in (3.19) can be calculated, and the cut-off spatial frequency v in (4.14) is chosen as the frequency whose magnitude is 10% of the maximum gain (i.e., r = 0.1 in (4.14)) [78]. This frequency v in (4.14) can be considered a tuning parameter for the robustness of the controller. The design procedure for the wavelet-based reduced-order MPC is as follows: c c 1. W e first calculate the reduced dimensions from the cut-off frequency v c in (4.14). Chapter 4. Model Reduction 1 73 1 1 1 1 (a) - Actuator locations 1 u> 0 1 ' 1 "— 1 -0.5 —h-- - - + - / • • + • • \ •+ -J>-—• • -4— •—1 H h 1 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 C D locations [metres] 0 1 2 3 4 5 6 7 8 9 10 Spatial frequency [cycles/metre] Figure 4.16: Basis weight response to a slice lip actuator and its power spectrum Based on the cut-off frequency v , the output profile can be reduced from 1100 to 86. Similarly, according to the controllable singular vectors, the input profile can be reduced from 110 to 82. c 2. After choosing the wavelets, we then calculate the wavelet matrices W y and W u and obtain the reduced model G . r In this example, we choose sym8 as the wavelets for model reduction. According to the modified wavelet packets method in Section 4.3.1, we will keep scale level 4's 69 and scale level 6's 17 coefficients for the outputs in the wavelet domain, and scale level 2's 55 and scale level 3's 27 coefficients for the reduced inputs. 3. We then check the model reduction indices such as the distortion ratios rji in (4.35) and n in (4.37), condition number. h If they are good, then we go to the next step. Otherwise, we go back to step 1 by changing the cut-off spatial frequency v or step 2 by changing the wavelet filter length If. Figure 4.17 shows the singular values and the distortion indices wi in (4.36) and w in (4.38). Table 4.2 shows the reduction performance. From the indices in the table, the reduction performance is pretty good except the condition number of the reduction model, which is slightly large. This is due to the slightly large cut-off spatial frequency v (i.e., smaller r in in (4.14)) and the original model's c h c Chapter 4. Model Reduction 74 large condition number (it is 102.07). 0 20 40 60 80 100 120 Figure 4.17: The singular values and distortion indices of the reduced model 4. Finally, we apply the reduced model for the objective function in (4.41), calculate the controller performance K in (1.1) and computational time, and compare them with those using the original model. If a good compromise between the performance and the time is obtained, then we finish the whole design procedure. Otherwise, we return to step 1 and do it again until we achieve our requirements. Figure 4.18 shows the output profiles before and after control using the reduced model G and the original model G. The output profile before control shown in Figure 4.18a is a realistic disturbance generated from a software simulator of Honeywell. When the original model G e C ° is used in MPC, the paper quality index or control performance index K in (1.1) is decreased from 0.578 to 0.199, and the average computational time for each control action is 29.25 CPU seconds (on an Intel Pentium III 700MHZ CPU, 256Mb RAM PC). While the reduced model G is used, the control performance index K is 0.213, which means this reduced-order controller rejects less about 3.7% disturbance than the controller based on the original model, however the computational time is reduced about r 1100xll T Chapter 4. Model 75 Reduction Table 4.2: The model reduction performance criteria for the single-array system. Original Model Reduced Model Dimensions m x n(p x q) 1100 x 110 86 x 82 Condition number T in (1.7) 102.07 19.44 Distortion ratio rji in (4.35) 0 0.0455 Distortion ratio n in (4.37) 1 0.2299 Sparsity index w 0.0206 0.3451 Computational time 29.25 12.38 [CPU seconds] Performance index K in (1.1) 0.199 0.213 h s 57.7% (12.38 CPU seconds compared to the above 29.25 CPU seconds). Because there are some higher spatial frequency components than the cut-off frequency v in (4.14) in the output profile before control, these components can not be removed by any controller. c There exists trade-off between the controller performance and the model reduction. Figure 4.19 shows the performance index K in (1.1) and on-line computational time for various reduced models. From this figure, the reduced model dimensions 86 x 82 may be a good compromise between the performance and the on-line computational time. 4.7.2 Three Properties Controlled by Four A c t u a t o r A r r a y s This example considers the wavelets-based reduced-order CD-MPC for a working industrial paper machine. This paper machine is equipped with four actuator arrays and one scanning sensor. The actuator arrays are slice lip actuators ui, steam box u , rewet shower u , and induction heating actuators u . Their dimensions Uj in (2.23) and actuator spacings X in (3.19) are ni = 57, n = 65, n = 88, and n — 84, X = 114.3 mm, X = 94.2 mm, X = 69.6 mm, and X = 72.9 mm. Three controlled properties (basis weight t/i, moisture y , and caliper y ) are measured by the sensor. The dimension for each of the measurements is 1280. Due to huge memory consumption in saving the predictions of the outputs and the limited resource of our computer, we downsampled the output dimension to 320 in order to implement the MPC in our simulation test. This can be realized by changing the maximum of k in (2.5) from 1280 to 320. The parameterized structure g(x) in (2.3) for the spatial response model Gij in (2.21) has the following 2 3 4 uj u2 2 u3 3 u4 2 3 4 ul Chapter 4. Model Reduction 76 Figure 4.18: (a) The output profile before control, (b) The output profile after control using the reduced model, (c) Difference of two controlled output profiles using the original model and the reduced model, (d) Spatial power spectra of the output profiles (dotted: before control, solid plus 'circle': after control using the original model, solid: after control using the reduced model). Note that the a;-axis of (d) is enlarged for clear display. parameters, Ciij = 0.00647 0 0 0 0.00502 -0.02359 0.04872 0 0.00549 -0.00450 0.00794 -0.01758 (4.48) 1.2 0 0 0 1.2 7 7 0 7 7 7 7_ (4.49) " 40 0 0 0 " 51 53 50 0 94 26 50 _ . (4.50) " 0.3 0 0 0 " 0 0.3 0 0 0 0 0.15 0 (4.51) 1 Pij = > 1 4 Chapter 4. Model Reduction 77 Figure 4.19: Controller performance (standard deviation) and computational time for different reduced models. Please note that there is no interaction between the basis weight profile yi and steam boxes u ; the basis weight profile y\ and rewet shower u , the basis weight profile y± and induction heating actuators u , and the moisture profile y and induction heating actuators u . That is why there are no identification parameters in the equations (4.48)(4.51). Other parameters for the subplants' dynamics hij(z) in (2.21) are as follows, 2 3 4 2 4 0.3679 0 0 0 0.5057 0.7575 0.4895 0 0.9024 0.6739 0.3679 0.9374 (4.52) 2 0 00 2 10 0 2 10 2 (4.53) Then, we follow the procedure described in Section 4.6 and get the reduced model. Through trial and error, the reduced dimensions for the outputs are p = 100 with i = 1,2,3, and the reduced dimensions for the inputs are qi — 50, q = 57, q = 77, and q = 74. Therefore, the dimensions for the original models dj in (2.21) and the reduced t 2 4 3 Chapter 4. Model models G rij Reduction 78 in (4.44) are rrii x Uj = 320 320 320 _ " 100 Pi x qj = 100 100 (4.54) 57 65 88 84 x (4.55) 50 57 77 74 Finally, we obtain the results. Figure 4.20 shows the computational time vs the performance for the MPC using the full model and the reduced model. Table 4.3 gives the detailed numbers for these results. The profiles before control are generated from a software simulator of Honeywell which is used for simulating realistic disturbances of real CD processes. From the numbers in Table 4.3, it is worth sacrificing a few percents of performanceforalmost half computational time reduction. In addition, it may be not implementable for the multiple-array system MPC using the full model because the sampling time for each scan is 30 seconds in this case and smaller than the on-line computational time shown in Figure 4.20 or in Table 4.3. Table 4.3: Performance and computational time for the MPC using the full and reduced models Before Control K Performance y i Ky 2 Ky 3 Average Computational Time per scan [CPU second] 4.8 0.4741 0.7861 0.6635 - After Control Full Model Reduced Model Change 0.0875 0.0894 0.39% 0.1337 0.1607 4.14% 0.0370 0.0431 0.98% 32.55 17.18 -47.22% Summary In this chapter, a novel model reduction method based on wavelet packets has been presented for paper machine cross-directional processes. This model reduction method has been connected with CD process controllability. The reduced model can effectively capture the controllable spatial components and reject uncontrollable components. Some practical issues in implementing the reduction such as distortion correction and choice of wavelets Chapter 4. Model Reduction 79 100 ' 0 ' 50 ' 100 scans 1 150 1 200 10 1 i 0 50 i 100 i 150 I 200 scans Figure 4.20: (a) Performance and (b) computational time for the multiple-array system MPC using the full model and the reduced model. have been investigated. Properly choosing the wavelets and their filter lengths will reduce the distortions. A practical procedure for model reduction of multiple-array systems has been provided. The constraints in the wavelet domain can be handled by using an inverse wavelet transform, and this method can guarantee that the constraints will not be violated in the original spatial (physical) domain and then no picketing will appear in the actuator profile. The major benefit of implementing MPC in the reduced domain is that the online computational time can be greatly reduced compared to the implementation of MPC based on the original model while good control performance is preserved. This has been demonstrated by two examples from real paper machine CD processes. Another minor benefit is that the controller implemented in the reduced domain will have some degree of robustness because it does not attempt control of the non-robustly controllable directions, that have been removed from the reduced model. However, this method can not guarantee robust stability of the closed-loop system. This issue will be explored in the next chapter. Chapter 5 Robustness and Performance In Chapter 3 , the large MIMO transfer matrices T (z) in (3.4) and T (z) in (3.6) could be conveniently transferred into a family of SISO transfer functions t (v!f, z) in (3.30) and t {v j, z) in (3.29) across the spatial frequencies. In this chapter, these SISO transfer functions will be used for deriving the robust stability condition for both the unstructured uncertainty A G I ( Z ) in (2.7) and the structured uncertainty A G ( Z ) in (2.8) described in Section 2.1 of Chapter 2, and will also be used for evaluating the closed-loop system's performance. Tuning of the MPC controller's parameters in the two-dimensional frequency domain will be demonstrated in Section 5.3. The method is extended to multiple-array systems in Section 5.4. Finally, several examples will demonstrate that both proper consideration of the structured uncertainty and careful design of the tuning parameters may greatly remove unnecessary conservativeness from the traditional controller design. ud yd ud V yd 2 5.1 5.1.1 Closed-Loop Robust Stability for Single-Array Systems Robust Stability Condition for the Unstructured Model Uncertainty The standard robust stability (RS) condition in (3.3) can be rewritten for the unstructured model uncertainty and may be transferred into the two-dimensional frequency domain. If T (z) in (3.4) is an RCM, then, based on (3.26) in Corollary 1, the maximum singular value of T (z) is equal to the maximum value of \t (Uj,z)\ in (3.30). In terms of the approximation t {v^z) the robust stability condition (3.3) may be rewritten as ud ud ud ud H T ^ ^ A d ^ H o o < a{T {z))<j{&ai{z)) < 1 O max(\t W,e ™)\) i2 ud ud < 1{UJ) j W " , u. J (5.1) This allows us to define a robust stability margin for unstructured additive model uncertainty as n := 7P= rr, (5.2 maxsup(|t (^,e ^)|) i2 ud for i/V e - K , i/ , • • • , vl%} and u e [0,u ]. u 2 N 80 Chapter 5. Robustness and Performance 5.1.2 81 R o b u s t Stability C o n d i t i o n for the S t r u c t u r e d M o d e l U n certainty From (3.28) and (3.43), the RS condition in (3.3) for the structured uncertainty A (^) in (2.8) can be rewritten as G2 p(T (e ™)A {e ™)) ud i2 G2 = i2 where Uj e {^i,^,--- > K} spatial frequencies with engineering units in (3.18) and co € [0,cu ] with the dynamical Nyquist frequency u . In the above derivation, from (5.3) to (5.5) we use the fact that F and F are unitary Fourier matrices, and eigenvalues are invariant under similarity transformations, that is, T A and F~ TudAF have the same eigenvalues [75]. From (5.5) to (5.6), we used (3.28) and (3.43) to conclude that T (e ™) • A(e ™) is a diagonal matrix. Note that for structured uncertainties the RS condition (5.6) is less conservative as it takes directionality into account. The unstructured RS condition (5.1) is equivalent to a r e N N n m l ud n 1 ud i2 i2 max(|?/„(z/\e )|)max(|c5K,e ; 1 r f K )|)<1, Va; => (5.6) (5.7) which can be significantly worse than (5.6) since it "pairs-up" the two maxima. Figure 5.1 demonstrates the conservativeness of the controller based on unstructured uncertainties. In this figure, t (v, e ™) and 5 (v, e ) comes from Figures 3.11 and 3.8 respectively (the steady-state case UJ = 0 is shown for clear display). In this figure, as max(|t (i/, , 1)|) max(|£(i/£, 1)|) = 50.4043x0.0280 = 1.4129 > 1, it seems that the closedloop system is unstable based on the RS condition (5.7) for unstructured uncertainties. However, max(\t (v", e ™) • 5(u , e ™)|) = 0.8497 < 1, and that means the closed-loop i2 ud ud l2nuJ u U(i i2 ud u 3 i2 J system is stable for structured uncertainties. Another benefit from (5.6) is that even for the same dynamical frequency UJ, the maxima of tud(Uj, e ™) and 6{v ,e ™) usually do not happen in the same spatial frequency v. In other words, given 5(Vj,e ™), it is possible to adjust t (u^, e ™) for good robust i2 l j i2 i2 !Note that = F~ [17, 24]. l t2 ud 82 Chapter 5. Robustness and Performance Spatial frequency [cycles/metre] Figure 5.1: A n example shows that the controller satisfies the R S condition (5.6) for structured uncertainties but violates the R S condition (5.1) for unstructured uncertainties. margin and performance, especially through tuning the weighting matrix Q i n (2.26). 3 This will be explored in detail in Section 5.3.1. 5.2 Closed-Loop Performance Indices The two-dimensional frequency sensitivity function t {u -,z) i n (3.29) may be used for yd v calculating the control system's bandwidth. The bandwidth i n the two-dimensional frequency domain can be similarly defined as that in the dynamical frequency domain [75]: Definition 2 The closed-loop bandwidth B (v?j,uj) w {fj,uj) such that \t d(v?,e ™)\ of CD processes is the contour in < -j- for all {vV,u) enclosed by B (V}J,LU). i2 w y Figure 5.2 illustrates the two-dimensional frequency bandwidth contour B (uj,u) w as the appropriate contour of Figure 3.10. In order to effectively evaluate C D control performance, we may define a better C D performance index x based on the two-dimensional frequency bandwidth: Definition 3 Let cu and v s s spectively when \t d(Vj, e be the dynamical frequencies and the spatial frequencies re- )\ first crosses over 0.7071 from (0,0) to the dynamical Nyquist l2nui y frequency UJ and the spatial Nyquist frequency of the inputs v^, then the CD control pern 83 Chapter 5. Robustness and Performance > Sc a) 1 2 3 4 dynamical frequency co [cycles/second] 5 x 1CT Figure 5.2: Contour plot illustrating the closed-loop bandwidth contour in the twodimensional frequency domain from Definition 2 for the closed-loop transfer function t d(Vj,z) in (3.29) with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H = 1, Qy = AI, Q = 10" 7, and Q = 0.002/ in (2.26) (compare with the surface plot in Figure 3.10). y 4 p U 2 3 formance index is defined by S3" * = £!£„"'(0-7071 ' | V K , e i 2 " " ) I K j 4 , 0.707L <'> 5 8 where d v and d denote frequency resolutions. The physical meaning of this index is that the volume of truly removed disturbance by the controller in the two-dimensional frequency domain is compared with the volume of theoretically removable disturbance, in other words, it truly reflects how much disturbance rejected by the controller. Please refer to Figure 3.10 for better understanding. A larger x means better control performance, i.e., more disturbances are removed. v 5.3 u Role of MPC Tuning Parameters Typically the choice of the horizons H and H in (2.26) is based on the dynamics of the process [2, 55, 76] and is not explored further here. In practice, H > 1 is rarely used as it leads to heavier on-line computational burden for solving the QP problem in (2.26) [7]. The number of degrees of freedom in the selection of tuning parameters in (2.26) can p u u Chapter 5. 84 Robustness and Performance be reduced by noting that without loss of generality, Qi in (2.26) may be fixed and the closed-loop performance modified via Q and Qz [76]. Typically in MPC, the weight matrix Q is considered to affect the dynamics of the closed-loop system while Qz affects the slower behavior and steady-state performance. Figures 5.3 and 5.4 illustrate the effects of changing the weight Q while Figures 5.5-5.7 illustrate the effects of changing the weight Qz on the closed-loop system in terms of the closed-loop transfer functions t d{v j,z) in (3.29) and t (uj,z) in (3.30). Figure 5.3 illustrates that increasing the weight Q has the effect of depressing the peak of the sensitivity function t {Vj,z) in (3.29) while reducing its bandwidth (only shows the spatial frequency v — 0 case). Figure 5.4 confirms this using a contour plot of the two-dimensional bandwidth and illustrates that in this case, Q has a negligible effect on the spatial component of the bandwidth B (i>j,u) in Definition 2. Figure 5.5 illustrates that decreasing the weight Qz helps to remove more low spatial frequency disturbances (only shows the steady-state case u = 0) because the sensitivity function t d(Vj, z) in (3.29) is linked to the disturbance rejection as mentioned in Section 3.1 of Chapter 3. However, the closed-loop transfer function t d(Vj, z) in (3.28), linked to the robust stability margin n in (5.2), has a larger peak for smaller Q as shown in Figure 5.6 (which is very similar as Figure 2 in [20]). Figure 5.7 confirms that changing Qz has a significant impact on the spatial bandwidth while in this case Q has a little effect on the dynamical component of the two-dimensional bandwidth B (V?,UJ) in Definition 2. The above results are in qualitative agreement with the traditional intuition regarding weights in the MPC. The above two-dimensional frequency analysis technique is useful for giving a consistent criterion for the selection of the MPC weights based on standard robust control theory, as compared to an ad-hoc design. 2 2 2 V ud y 2 yd 2 w y U 3 3 W 5.3.1 Structure of the Weighting Matrix Q 3 If the weighting matrix Qz in (2.26) is a diagonal penalty function as it was for Figures 5.5 and 5.6, it uniformly penalizes the actuator profile u(z) from low spatial frequency to high spatial frequency. Then, a smaller Q gets good performance and bad robustness, and vice versa. Given a model with structured uncertainty like in Figure 3.8, it is not beneficial to uniformly penalize the actuator profile from low to high spatial frequency. Also from a practical view, the most important task for a CD controller is to remove low spatial and dynamical frequency disturbance as much as possible. Can we design a penalty function Qz so that \i d(uj, e ) \ in (3.30) has a peak in low spatial frequency range and no peak in middle or high spatial frequency range? This kind of Qz would give a controller with good compromise between performance and robustness. The band diagonal bending moment matrix B in (2.31) is a good candidate. 3 u M i2voJ 85 Chapter 5. Robustness and Performance 1.8 1 Ql 0 . L _ 1 . . 0.01 0.015 0.02 0.025 0.03 dynamical frequency co [cycles/second] 0.005 1 0.035 Figure 5.3: Illustration of the effect of changing Q on the gain of the closed-loop transfer function t d(i>j,z) in (3.29) at spatial frequency = 0 with respect to dynamical frequency with the nominal process model G(z) in (2.1) and the controller parameters H = %,H = \,Qi=U, and Q = 0.002/ in (2.26). 2 y p u 3 Q -0.011 Q-0.1I 2 / 7/ " h • £i. |t Jv ,e y i&0> )|<0.7071 • % ,S 1 '< is < : J 1 0f 1 it 2 0 -o.oou 2 / f i i i a i ii 3 4 5 6 7 8 dynamical frequency <o [cycles/second] 1 i 9 U 10 x 10~ 3 Figure 5.4: Contour plot illustrating the effect of changing Q in (2.26) on the closed-loop bandwidth contour B (v?,u) in the two-dimensional frequency domain for the closed-loop transfer function t {v j, z) in (3.29) with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H = l,Qi = 41, and Q = 0.002/ in (2.26). Note the greater influence on the dynamical component of B {V^,UJ). 2 w yd V p U 3 W 86 Chapter 5. Robustness and Performance spatial frequency v [cycles/metre] Figure 5.5: (cf. Fig. 1 in [20]) Illustration of the effect of changing Q in (2.26) on the gain of the closed-loop transfer function tyd{Vj,z) in (3.29) at steady-state with respect to spatial frequency with the nominal process model G(z) in (2.1) and the controller parameters H = 8,H = 1, Q = 41, Q = 10~ I in (2.26). 3 p U v A 2 Q =0 3 00021 _Q_ Q =0 0021 021 ^ - 0 p- 100f- 1.5 2.5 2 3 3.5 spatial frequency v [cycles/metre] Figure 5.6: (cf. Fig. 2 in [20]) Illustration of the effect of changing Q in (2.26) on the gain of the closed-loop transfer function t (u^,z) in (3.30) at steady-state with respect to spatial frequency with the nominal process model G(z) in (2.1) and the controller parameters H = 8, H = l,Qi = 41, Q = lO" / in (2.26). 3 ud 4 p u 2 87 Chapter 5. Robustness and Performance 4.5 F dynamical frequency co [cycles/second] x 10 3 Figure 5.7: Contour plot illustrating the effect of changing Q in (2.26) on the closed-loop bandwidth contour B (v?,cu) in the two-dimensional frequency domain for the closed-loop transfer function t d(Vj, z) in (3.29) with the nominal process model G{z) in (2.1) and the controller parameters H = 8,H — 1,Q = AI, and Q = 10~ / in (2.26). Note the greater influence on the spatial component of B (Uj,uo) as compared with Figure 5.4. 3 w y 4 p u X 2 w The bending moment matrix B in (2.31) includes the first order derivative near the edges and the second order in the middle. If let B (l, 1) = B (n,n) = —2, B (m, 1) = B (l,m) = 1, then B is an RCM. Its spatial frequency frequency representation m m m m m m {6miK)> • • • ,LnK)} = BlAG(F B F^). n (5.9) m Figure 5.8 shows 6 (i/") in the spatial frequency domain (the actuator number is 64). The bending moment matrix B looks like a high-pass spatial filter itself, but in the cost function (2.26), if Q = q xB x B, (5.10) m m T 3 3 m l where q E R is a tuning number, then, Q works like a low-pass filter to remove high spatial frequency from the actuator profile u. This Q has low gains in low spatial frequency and high gains in high spatial frequency, i.e., it lightly penalizes the low spatial frequency part of u and heavily penalizes the high spatial one of u. Table 5.1 and Figure 5.9 show that Q using B in (5.10) gets both better performance and better robustness than using diagonal penalty function. In Table 5.1, the performance index x = 0.1579 and the robustness margin index p(T„dA) = 0.8120 for Q using B in (5.10) compared to x = 0.1443 and p(T A) = 0.8529 using diagonal Q . In Figure 1 3 3 3 3 m 3 ud 3 m 88 Chapter 5. Robustness and Performance - / - - / )l / / / - / / / / / - 1 1 1 1 1 1 1 0.5 1 1.5 2 2.5 3 3.5 = ^ 0 / / / spatial frequency [cycles/metre] L 4 Figure 5.8: The spatial frequency representation |6 (^")| of the bending moment matrix B in (2.31). m m 5.9, Q = 0.002 x / for diagonal Q and Q = 0.0007 x B^x B in (5.10) using B , and other tuning parameters in (2.26) are H = 8, H — l,Qi = 47, Q = 10~ J. Please note that the figure only shows UJ = 0 case for clear display. Definitely, there are ways to design better high-pass spatial filters than B for better performance and robustness. The bandwidth of such filters can be adjusted while that of B can not. This issue may be considered as a future research objective. 3 3 3 m m 3 p u 2 m m Table 5.1: Robustness and performance: Q using different structures Using B Diagonal Index Q = 0.002/ Q = 0.0007B£ x B Robustness p(T A ) in (5.6) 0.8529 0.8120 Performance x (5-8) 0.1443 0.1579 3 m 3 ud 3 m G2 m 5.4 Robust Stability for Multiple-Array Systems In this section, we will derive the RS condition for multiple-array systems. In order to simplify the following development without losing generality, we will let the numbers of actuator arrays and controlled properties N = N = 2. The structured uncertainty for multiple-array systems can be naturally extended from y u Chapter 5. Robustness and Performance 89 spatial frequency v [cycles/metre] Figure 5.9: The sensitivity functions \t Q vs Q using B in (5.10). (vV, yd 3 3 1)| in the spatial frequency domain: Diagonal m the single-array systems as described in Section 2.1 of Chapter 2. For a two actuator arrays and two controlled properties CD system, the plant model can be written as GP(*) = G(z) = A(z) = G(z) +A(z) e c 2 Gn(z) G (z) G2i(z) G (z) AnW A m x ( n i + " 2 ) (5.11) ) 12 (5.12) 22 1 2 l » (5.13) A (z) &2i(z) 22 where each A - satisfies (2.8) but is permitted to have different uncertainty limits in (2.13). Note that the coupling interactions between two properties Gu{z) and G {z) and their model uncertainties A ( z ) and A i(2) are usually neglected in the traditional CD control strategy (i.e., one actuator array controls one property). For a multiple-array system, the control sensitivity function T (z) is defined in a fashion similar to that in (3.4) for a single-array system: i3 2y 12 2 ud T ^ O ) = -K(z) [I * 2m (5.14) + G(z)K( )]- , 1 2m 2 where K(z) G C ) is from [3] (and may be found in (6.64) of next chapter), and G(z) € C ( ) is in (5.12). Similar to the partition of the additive structured uncertainty A(z) in (5.13) for a multiple-array system, the transfer matrix T (z) can be ( n i + n 2 2mx x 2 m ni+n2 ud 90 Chapter 5. Robustness and Performance split into different sub-matrices (the z argument will be ignored in this section for clarity of notation), r Rn • ud £ (£( l+ 2)x2m n Rn n (5.15) i?2i R22 where each of sub-transfer-matrices GC with i = 1,2 in (5.15) is an RCM. Now let us derive the RS condition in (3.3) for multiple-array systems. The calculation of the eigenvalues of T A. niXm ud \{T A) M = A ud M a b M (5.16) M c d M F, ni 0-. F H b M M c •n-2 (5.17) F H d n.2 Ma M M M c M a (5.18) b d where 2 22 (5.19) 22 22 (5.20) #i A , = i2 iA + i? A , M = RuAu + RuAn, M = R iA + M = F M F£, _M = F M F%, (5.21) M = F M F», M = F M F^. (5.22) a c a c 2 ni n2 n a c M E22A21, b d = R b M d n 2 ni A u + 12 b n2 d The complex Fourier matrices F and F have dimensions n\ xni and n xn respectively. Because all R - and A^- are RCMs, based on the RCM's properties [24], then M G C and M G C are each an RCM, and M and M are diagonal matrices. Similarly, because all R^ and A j are RCMs and their column vectors are bandlimited to the spatial Nyquist frequencies of the inputs v%, although M G C" and M G C are not RCMs due to non-integer n /ni (assuming without loss of generality that n > ni), M and M have a similar "diagonal" structure [24] as T {z) in (3.28), ni n2 2 2 niXni io a n2Xn2 d a d t lXn2 n2Xni b c 2 2 c b ud j = i + k, (5.23) otherwise, j=i + k, (5.24) otherwise, where k = 0 for i — 1, • • • , ni/2 and k = n — ni for i = n-i/2 + 1, • • • ,ni for even ni or k — 0 for i = 1, • • • , (ni + l)/2 and k = n - ni for i = (ni + l)/2 + 1, • • • , ni for odd ni. Now the matrix in (5.18) is highly structured since M and M are both diagonal 2 2 a d 91 Chapter 5. Robustness and Performance matrices, and M and M are composed of diagonal blocks as shown in (5.23) and (5.24) respectively as illustrated in Figure 5.10. b c 0 20 40 60 80 100 120 Figure 5.10: An example of T A which shows the nonzero elements. ud The RS condition (3.3) for a multiple-array system with structured uncertainties can be rewritten as p(T A) = max ud 3 Ma Mb M M c < 1 (5.25) d for all A given in (5.13). Appendix E exploits the block diagonal structure as shown in Figure 5.10 to derive an analytical formula for the eigenvalues. By contrast, the RS-condition (3.3) for a multiple-array system with unstructured uncertainties is a(T )a(A) < 1 & a(T ) < ud 5.5 ud (5.26) Examples In this section, we will demonstrate the effectiveness of the proposed method through three examples. The first is for a single-array system with unstructured uncertainties, the second is for a single-array system with structured uncertainties, and the third is for a multiple-array system with structured uncertainties. 92 Chapter 5. Robustness and Performance Table 5.2: The performance criteria for the conservative and aggressive controllers Conservative Qi = 1(T /, Q = 5 x 10" / 3 3 3 Stability margin r] in (5.2) F(QL) M Model uncertainty tolerance" Spatial bandwidth B (v?,0) [cycles/metre] Dynamical bandwidth B (0,UJ) [cycles/second] Performance index i n (5.8) Paper quality index K in (1.1) of the outputs y{k) in Figure 5.11 Q Aggressive = 10" /, Q = 9 x 10" / 4 2 4 3 0.0310 0.0134 0.0306 0.0130 81% 35% 2.51 2.79 w 4.65 x K T 4.71 x 10" 3 3 w 0.1332 0.2192 0.9370 0.9012 X Y 4 "IIGOJOHOO =0.0382. K,d = 1.2984, here d denotes the disturbance profile before control. b 5.5.1 U n s t r u c t u r e d M o d e l Uncertainty Case In this section, we will show how the closed-loop bandwidth B {V^,UJ) W in Definition 2, the performance index x hi (5.8), and the robust stability margin 77 in (5.2) can be used for practical tuning of the cross-directional model predictive controller (2.26). In this example, the nominal model G(z) in (2.1) with the parameters in Section 2.1 of Chapter 2 is used. Its two-dimensional approximation <?(^J, z) is illustrated in Figure 3.5. The model uncertainty in this case is assumed to be constant across dynamical and spatial frequencies with l{yj) = 0.0115, which is 30% of the (in this case ||G(f)||oo norm of the nominal model G{z) = 0.0382). Through trial and error, both a conservative and an aggressive controller were determined. For the conservative case, the tuning parameters in (2.26) are the horizons H = 8, p H = 1, the weights Qi = Al, Q = 10~ I, and Q = 5 x 10~ /. For the aggressive case, 3 u 2 3 3 only the weights Q — 10~ I and Q = 9 x 10~ / while the other parameters are un4 4 2 3 changed. The performance criteria B (VJ,OJ) w and n developed in this paper are found in Table 1. From the table, the approximate stability margin n is a close match to the true stability margin l/a(T ) ud for this example. Table 1 also shows that the aggressive con- troller has larger bandwidth but smaller model uncertainty tolerance than the conservative controller. The simulation uses a steady-state disturbance d(z) in (3.6) as shown in Figure 5.11a and target profile y sp = 45 gsm in (2.26). This disturbance profile d(z) comes from a real working paper machine. Using (3.6), Figure 5.11c and 5.lid show the output profile Chapter 5. Robustness and Performance y(z) and its spatial spectrum 2 \y(z)\ = \F 93 m 4 6 CD locations [metres] • y(z)\ at steady-state u 1 2 = 0 after using the 3 4 v [cycles/metre] Figure 5.11: (a) and (b): the disturbance profile d{z) in (3.6) before control and its spatial spectrum \d(z)\ = \F • d(z)\; (c) and (d): the basis weight profile y{z) in (3.6) and its spatial spectrum \y(z)\ — \F • y(z)\ after using the conservative controller; (e) and (f): the basis weight profile y(z) in (3.6) and its spatial spectrum \y{z)\ = \F • y{z)\ after using the aggressive controller (note that (b), (d), (f) have different vertical scales). m m m conservative controller while Figures 5.lie and 5.11f display the analogous plots after using the aggressive controller respectively. Note that the vertical axes in Figures 5.11b, 5.lid, and 5.1 If are displayed on different scales in order to clearly show the difference in the performance achieved using the aggressive and the conservative controllers. From the numbers of the estimated steady-state standard deviation in Table 1, the aggressive controller removes about 9.9% more disturbances than does the conservative controller (the open-loop variation was k = 1.2984). An improvement in performance is to be expected as the spatial bandwidth B (i>y, 0) and the performance index x were increased d w 94 Chapter 5. Robustness and Performance - from 2.51 cycles/metre to 2.79 cycles/metre and from 0.1332 to 0.2192 respectively. Note that only the spatial components lower than the actuator Nyquist spatial frequency v'ff = 4.54 cycles/metre are shown in this Figure and those higher than this frequency are not controllable. 5.5.2 Structured Uncertainty vs Unstructured Uncertainty In this example, the nominal model G(z) in (2.1) with the parameters in Section 2.1 of Chapter 2 is used again. This time we assume the parameter uncertainty limits Si = 0.25 and S = 0.6 in (2.13) respectively. Then, a class of perturbed models G (z) can be obtained by (2.8). Therefore, the maximum singular value of the additive uncertainty A (z) in (2.8) for all dynamical frequencies is 0.0307 (this can be obtained by performing a grid-search, i.e., let each parameter such as 7 in (2.13) run across its range by small steps). This number is much larger than 30% of the maximum singular value of the nominal model G(z) as we assumed in the previous section for the unstructured uncertainties. If we were to use the RS condition for unstructured uncertainty in (5.1) the we would require 2 p G2 max(|i (^,e^)|) < = 32.6042. ud (5.27) To satisfy the RS condition for unstructured uncertainty requires a controller K(z) such that a(T ) = a(-K{I + GK)~ ) < 32.6042. Using trial and error, a controller Ki{z) is found such that a{T (z)) = 32.3825. Figures 5.12a and 5.12c show the two-dimensional frequency representations t (uj,z) and t (Uj,z) o(T (z) and T {z) respectively. From Figure 5.12c, it is obvious that a lot of low frequency disturbances will not be removed by this controller based on the unstructured model uncertainty assumption. In order to incorporate the structure of model uncertainty, we first find the maximum absolute values surface of 5{v™, e ) in the two-dimensional frequency domain (for example, the dashed line in Figure 3.7f shows the uncertainty bounds surface at steady-state), then let t^iy^e ^) satisfy the RS condition (5.6). Again using trial and error tuning, another controller K {z) is found and maxsup(|f (i/\e ™) • S(vf,e )\) = 0.9912. In order to compare the two controller's performance, Table 5.3 shows the results about the RS condition from (5.6) and (5.7). The controller Ky{z) which is designed based on the unstructured uncertainties is too conservative for the structured uncertainties. Figures 5.12b and 5.12d show t (Vj,z) and t (i>j,z) in the two-dimensional frequency domain respectively for the controller K (z) designed for the structured uncertainty. l ud 1 ud ud ud yd yd l2nu> 12 2 i2 2 ud ud i2nu yd 2 K {z) has parameters: H = 8,H = l,Qi = 4J,Q = 10- J,Q = 0.007/ in (2.26). 1 3 1 K (z) 2 2 p U 2 has parameters: H = 8,H = l,Q =AI,Q p U x 3 = 10- /,Q = 0.00057 in (2.26). 4 2 3 Chapter 5. Robustness and Performance 95 Figure 5.12: (a) and (c): Illustration of £ ( t f , e ™) in (3.28) and t {v), e ™) in (3.27) in the two-dimensional frequency domain which satisfies the RS condition (5.1) for the unstructured uncertainty case; (b) and (d): Illustration of t (i>'!-, z) in (3.28) and t {v -, z) in (3.27) in the two-dimensional frequency domain which satisfies the RS condition (5.6) for the structured uncertainty case. i2 i2 ud yd v ud yd It is evident from comparing Figure 5.12c with Figure 5.12d that the controller based on the structured uncertainty K {z) has much better nominal performance than one based 2 on the unstructured uncertainty Ki(z). The main reason is that the maximum singular values of t (vj, e ™) and 5(uj, e ™) occur neither at the same dynamical frequency OJ, %2 12 ud nor at the same spatial frequency z/". We assume the M P C controller is unconstrained when we analyze the closed-loop transfer functions T (z) and T^z) yd in Chapter 3 and the above sections. Here we will test whether the closed-loop system is stable or not under constrained cases. We use the same tuning parameters above which are designed for the structured uncertainties and then test if the controller is stabilizing for all the possible models belong to n in (2.8) with 2 parameter limits 5\ = 0.25 and 8 = 0.6. Figure 5.13 shows the test results. Figure 5.13a 2 shows the realistic disturbance profile before control which is generated from a software 96 Chapter 5. Robustness and Performance Table 5.3: RS and performance test resultsforstructured and unstructured uncertainties (the single-array system) Ki(z) RS and performance Test o(T (e ™))o(A(e"™)) ud i2 max sup ( t {vi,e ™)b{vj,e ™) a ud ) i2 (for unstructured K {z)(for structured uncertainties) uncertainties) 0.9941 3.8739 (fails RS) 0.0943 (very conservative) 0.9912 2 Performance index X in (5.8) 0.0912 0.1950 simulator of Honeywell. Figures 5.13b and 5.13c show a class of the output profiles due to different possible models after control and the variances of the output profiles with time respectively, while Figures 5.13d and 5.13f show a class of the input profiles due to different possible models after control and their variances with time respectively. Note that the bending limits U in (2.30) is 20forthis case and from Figure 5.13e the actuator array apparently hits the limits. However, the closed-loop system is stable from Figures 5.13c and 5.13f. max 5.5.3 M o i s t u r e Controlled by Steam B o x and Rewet Shower A c tuators In this example, we analyze a system comprised of two actuator arrays for controlling the moisture profile of a newsprint sheet. The first array consists of H i = 64 independently operated steam showers for drying. The second array is made of n = 80 independently operated water showers to rewet the sheet and prevent over-drying. The moisture profile is measured by a scanning sensor at m = 320 locations across the sheet. The nominal model in G(z), the inputs U(z) and the outputs Y(z) in (2.20) are 2 G(z) 320x144 [Gnhn(z) 12 l2 u {z) x U(z) (5.28) G h (z))eC {z)eC \ b Ul u (z) u (z)eC , m 2 (5.29) 2 320 (5.30) yi(z) e C Y(z) The nominal parameters of Gu and G\ are 7n = —0.02358, «n = 7, £n = 13.25, /?n = 0.10, c = 2.5 + 5(i - 1) for % = 1, • • • , 64, and 7 = 0.04872, a = 7, £ = 7.5, /?i2 = 0.10, Ci2» = 2+4(^ — 1) for i = 1, • • • ,80. The parameters of the dynamical responses h (z) and 7T, (^) are a — 0.5057, q = 1 and a = 0.7575, q = 0. Although these 2 ni n 12 12 n n i2 12 l2 12 Chapter 5. Robustness and Performance Actuator numbers 97 Actuator numbers [ time s c a n s l Figure 5.13: Application of the CD-MPC controller designed for the structured uncertainties to constrained cases: (a). The output profile before control; (b) a class of output profiles after control; (c) a class of variations of the output profiles with time; (d) a class of input profiles after control; (e) the second order derivatives of the input profiles after control; (f) a class of variations of the input profiles with time (note that the y-axes are different in (a) and (b) for clear display. parameters come from simulated CD processes, they are consistent with those identified on real paper machines. The parameter uncertainty limits 6i and 6 in (2.13) for both the 2 nominal responses are assumed as 0.25 and 0.6 respectively. Based on the above parameters uncertainty, the maximum singular value for all the additive structured uncertainty a(A(e )) i2nu = 0.0868 in (5.13) (i.e., norm of A(z)). Therefore, the robust stabil- ity condition (5.26) for the additive unstructured uncertainty (in fact it is structured) is o(T (e ™)) i2 ud < 11.5229, where a controller K l » e C 1 4 4 x 3 2 0 3 was obtained through trial and error tuning . Figure 5.14a shows that all the singular values of T (e ' ' ) in (5.15) 1 !2 r J ud satisfies the RS condition (5.26). In this figure, the maximum singular value is 11.0993. Note that the discontinuity in this Figure is due to the fact that there are two actuator arrays. The sensitivity function t (v,e ™) %2 yd in the two-dimensional frequency domain is K (z) has parameters: H = 8, H = 1, Qi = il, Q = 10" /, Q = 0.006/ in (2.26). l 3 3 p u 2 3 Chapter 5. Robustness and Performance 98 illustrated in Figure 5.14c. 5.14c. Figure 5.14: (a) and (c): Illustration of the singular values of T d(z) in (5.15) in the dynamical frequency and singular vectors domain and the sensitivity function t (Vj, z) in (3.27) in the two-dimensional frequency domain for the unstructured uncertainty case; (b) and (d): Illustration of the singular values of T (z) in (5.15) in the dynamical frequency and singular vectors domain and the sensitivity function ty (u^z) in (3.27) in the twodimensional frequency domain for the structured uncertainty case (for the multiple-array system). u yd ud d From trial and error tuning, a controller K (z) GC 4 1 4 4 x 3 2 0 is found such that p(T A) 1 ud = 0.9838 in (5.25). In order to be comparable with the unstructured case, Figure 5.14b shows the singular values of T (z) ud in (5.15) as does Figure 5.14a. Note that while the maxi- mum singular value in Figure 5.14b is much larger (87.6746) than that in Figure 5.14a, it still satisfies the RS condition in (5.25) for the structured uncertainty. Table 5.4 also shows RS test results for the two controllers. The sensitivity function \t (Uj,z)\ yd at low spatial and dynamical frequency range is about 0.15 in Figure 5.14c, which means at least 15% disturbances belonging to this frequency range can not be removed by the controller K {z) has parameters: H = 8,H = \,Q = 47,Q = 10- J,<2 = 4.8 x 1CT / in (2.26). l 4 4 p U X 2 4 3 99 Chapter 5. Robustness and Performance Comparing Figures 5.14c with 5.14d, it is clear that the controller performance is much better when considering the uncertainty's structured characteristic, in other words, making use of the structure in real CD processes can remove unnecessary conservativeness. K (z). 3 Table 5.4: RS and performance test results for structured and unstructured uncertainties (the multiple-array system) RS and performance test o(T (e ™))o(A(e ™)) ud l2 l2 p(T {e ™)A(e ™)) ud i2 12 Performance index X in (5.8) K (z) (for unstructured K (z)(for structured uncertainties) uncertainties) 0.9634 7.6102 (fails RS) 0.6121 (very conservative) 0.9838 3 4 0.2585 0.3779 5.6 Summary The two-dimensional frequency bandwidth B (Vj,u), the performance index x, and the robust stability conditions for structured and unstructured uncertainties have been presented in the two-dimensional frequency domain as intuitive measures of closed-loop performance and robustness. The balance between robustness and performance of the CD-MPC can be achieved by tuning the weighting matrices Q and Q . Roughly speaking, the weight matrix Q affects the closed-loop dynamical behavior while the weight matrix Q determines the spatial performance. The example in Section 5.5.1 has shown that this two-dimensional frequency analysis method gives consistent criteria for selecting the weights of the CDMPC. Considering the structure of model uncertainty and the structure of the weight matrix Q can greatly improve the controller's closed-loop performance while maintaining the robustness. The robust stability condition for multiple-array systems has been derived in the two-dimensional frequency domain. To the author's knowledge, this is the first time this problem is explicitly addressed in the robust control framework for CD systems. Examples in Sections 5.5.2 and 5.5.3 clearly demonstrated the effectiveness of this method. The tuning parameters (weighting matrices) of the controller based on unconstrained MPC are usually satisfactory for constrained MPC. However, accurately predicting the steady-state performance is still needed for CD-MPC tuning, especially for the constrained case. This will be investigated in the next chapter. w 2 3 2 3 3 Chapter 6 Steady-State Performance Prediction As stated in Chapter 2, it is very inconvenient and time-consuming to run the large-scale CD-MPC simulation, especially for constrained cases, in order to evaluate the performance given by a set of tuning parameters. In this chapter, a static optimizer is proposed to optimally predict the steady-state performance of the large-scale CD-MPC closed-loop system [26]. First, we will introduce the state-space model for a multiple-array system and the objective function for MPC of this system. Then, we will propose a one-step static optimizer for predicting the steady-state performance. This objective can be achieved by designing the optimizing weighting matrices of the static optimizer. Finally, two examples from real working paper machines are given for demonstrating the effectiveness and efficiency of the static optimizer. 6.1 CD-MPC Implementation for Multiple-Array Systerns 6.1.1 Objective Function The CD-MPC quadratic programming (QP) problem for a multiple-array system is set up as for a single-array system in (2.26) and solved by mm AU(k) Hu-l + i ( AU k + i) Q2&U{k + T U{k + i) Q U(k (6.1) + i)} T 3 subject to equation (2.20) and the following constraints: U in {j ) m -M m a x (j) avgmax -AU m a x (j) < Uj(k) < < BmU) • j( ) < ^•avgjik) < Uj{k) u Umaxij), < k for j - 1) 100 for j = max av < for j' = AU (j), max (6.2) 1, • • • , JV„, M (j), — U gmaxi.j)i - Uj{k = 1, • 1, • • • , (6.3) N, u (6.4) , -/V, for j = u 1, • • • , NUl t (6.5) Chapter 6. Steady-State Performance Prediction 101 where Y G R y : output predictions; Y G R ' : set-points for controlled properties; H : prediction horizon; H : control horizon; For simplicity, H is chosen as 1 in the following development. Qi e R v v diagonal penalty function for controlled variables; Q G R^} i ^j=i i diagonal penalty function for the changes in manipulated variables; Q G R ^ i = i J ^ j = i j penalty function for manipulated variables; Qi, Q2, and Q are symmetric positive semidefinite. Umin(j), U ax{j) G R : lower and upper limits respectively; M (j) G bounds for bending limits; B (j) G R J J ' : band diagonal bending moment matrix; 3 j( ) , Cj = [1, • • • , 1] G R : average of the j actuator array; (k) u, •avgj U gmax(j)'- maximum value allowed for the average of the j actuator array; AUmax(j) ^ maximum change allowed between successive control actions. For each controlled property, there is a sub-penalty function Q , and for each actuator array, there is a sub-penalty function Q . In other words, Q i and Q have the following (block) diagonal structure: N Ny m m sp p u u N mxN m : =1 n x n 2 : n x n 3 : 3 nj m max n x n m e u k th nj th aV e nj: u 3j Q i = V (Qn,Q12, • • • , Q I N ) , U where V(X ,X , etc, that is, i i 2 ••• , X y N 3 Q = V(Q i, 3 } 6.1.2 2 Q , • • • , QSNJ , 3 2 (6.6) denotes a block diagonal matrix with diagonal elements Xi, X 0 0 0 Y V(Xi X , 3 • • • , XN) 0 0 X 0 0 0 0 XN 2 X N (6.7) State-Space Model It is convenient in MPC implementation to use a state-space model rather than a transfer matrix model as in (2.20). For each controlled property yi(z) = Ylf=i {Gijhij{z)uj(z)) + di(z), we can write the following state-space form X{j(k + 1) — A{jXij{k) + BijUj^k), (6.8) (6.9) 3=1 where i = 1, • • • , N , j = 1, • • • , N . y u Chapter 6. Steady-State Performance Prediction 102 Since AU, not U, is minimized in the objective function (6.1), the following augmented state-space model form is used: %aij{k + 1) — Vi{k) AaijX ij(ki) a — y ^C ijX jj(k) a (6.10) B ijAUj(k), + a (6.11) -\-dj{k), a 3= 1 where ^^X^j , CijXij Aij A •• — 6.1.3 Auj(k) 0 (6.12) = Uj(k) — Uj(k — 1), Bij Cij Bij B ij a a (6.13) 0/ aij Output Predictions The predicted future outputs Y{k +i) in (6.1) are a function of the current states, past and future changes in U of (6.1). The sub-plant state prediction can be described as [54]: Xijiji) — Pxij%aij(k) + (6.14) P {jAUj{ki), u where x j(k) is the estimated state at time k, and ai A {k) Uj Xijik) = Auj(k + 1) AUj(k) %aij Auj(k + H - 1) H -Bp) - u A„ 0 Baij P XI] (6.15) = AP H A P~ H l "•aij aij A 0 0 B, R A R .. A ~ R • "aq aij ^aij Hp 2 A H P — H •"•aij (6.16) TJ U a.ij J D Then, all sub-plants state predictions can be written by ±(k) = P x (k) x a (6.17) + P AU{k), u where X(jfc) ^n(fc), X (k), lNn X {k), Nyl X (k) NyNu T ,(6.18) Chapter 6. Steady-State Performance Prediction %a(k) = AU(fc) = x ii(k), x i (k), a ••• x i(k), a Nu I AUi(k), 103 x aNy ,(6.19) (k) aNyNu (6.20) AU (k), 2 Pun ••• 0 0 ftliV u (6.21) ft 0 PuNyl 0 ••• PuN Nu y ft — F> (P n,- • • , PxlN , x Px21, • • • , Px2N i u u • • • , PxNyl, , PxN N ) 1 1 1 y u (6.22) , where '£>' denotes the block diagonal operation defined in (6.7). Since the current states can not be directly measured, they have to be estimated from the previous states and the current measurements, L(i) x {k) = x (k) aij (6.23) + ^r(y;( aol3 f c ) - where the previous states x (k) can be inferred from aoij x (k) aoij = A x (k aij (6.24) - 1) + B jAuj{k aij ai - 1), = [0,1) is the estimator gain for each actuator array. The estimated outputs can be obtained by T N u (6.25) ]C jjX ij{k). Vi{k) — ^ a ao 3= 1 Then, the output prediction Y^k) can be obtained by Hk) Mk + 2) Cu Ci2 ••• X a X i ( k ) 2 ( k ) CiN u _ jjiik + Hp) (6.26) XiNuik) where Cij T) ( C ij, • • • , C ij ) , a a j 1)''') N u (6.27) Hn where '£>' denotes the block diagonal operation defined in (6.7). Therefore, all of the Chapter 6. Steady-State Performance Prediction 104 outputs prediction Y(k) can be obtained from Y(jfc) = = (6.28) C±(k), Y (k) Ny where X(fc) is in (6.17) and ClN C 6.1.4 0 o 0 0 0 0 u c 2 0 0 CN 1 2 0 U 0 0 0 (6.29) C NyN U J Transfer to Standard Q P Problem The objective function in (6.1) can be rewritten as J{k) = (Y(k) - Y (k)) Q,(Y(k) - Y (k)) + AV(k) Q AV(k) T + U(k) Q U(k), T sp sp (6.30) T 2 3 where Y(k) and AU(A;) are in (6.28) and (6.20) respectively, Y (k) Y i(k), sp Y (k), sp Y ••• , sp2 Y N (k) sp lT (6.31) y (6.32) Y i{k + 1), Y i(k, + 2), • • • , Y i(k + Hp) sp •* spi sp Ui(k), sp U (k), U (k) 2 Nu \ 1) + S =V U(k- (6.33) AV(k), du P -iT V(k) Uj(k), Ujik Ui(k-l), + l), Uj(k u (k-l) +H U u 2 N u 1) J , j = X (k-l) p P2 P (6.36) U (6.37) {Sdul,Sdu2, • • • , SduN ) , u ' I' > " I I I I PJ }H , U _I j = l, AU6.38) _ i I • • I > V{ Q i , Q U (6.35) V (U i, U , • • • , U N ) , U(k-1) ,N , (6.34) l s ••• , Qi ), (6.39) Chapter 6. Steady-State Performance Prediction 105 Qi = 2>( Q „ Q,, ••• , Q, ), > , (6.40) v From (6.30), after some algebra , the QP problem (6.1) can be transferred to the standard QP format: min 0.5AU(£;)HAU(A;) + </>AU(£:) (6.41) 1 r T AU(fe) subject to (6.42) QAV(k) < w(jfe) where H = P^C q CP + q +S^ Q S , T i u 2 u 3 (6.43) du 4> = -PjC Q Y (A:) + P^C^CP^k) + §? ® U U(k T 1 sp u 3 p - 1), (6.44) Q and u:(k) can be obtained from (6.2)-(6.5). 6.2 Static Optimizer As discussed in Section 2.5 of Chapter 2, the standard technique for tuning an MPC controller typically involves running simulations to analyze closed-loop performance. However the QP problem presented in (6.41)-(6.44) is typically of very large-scale and the Hessian matrix H in (6.43) can be as large as 600 x 600 in multi-array CD control applications. Even with modern computing power, the execution time of simulations is frustratingly long in industrial situations (up to 60 minutes). This section presents the framework for an optimization problem designed around the steady-state process model. The plant model (2.20) at steady-state (i.e. z = 1) is Y ss = GU SS SS +D , (6.45) SS Gu .•• ss G s s iN u G (6.46) s G sNyl • S = rZTT' i = •• G N NU SS V ,---,N , 1 y j = l,-..,N , u (6.47) where Y , U , and D are the outputs, inputs and disturbances at steady-state respectively, G is the steady-state plant model, dj and in (6.47) come from (2.21) and (2.25) respectively. ss ss ss ss 1 Refer to Appendix F for details. Chapter 6. Steady-State Performance Prediction 106 Assuming that the initial actuator profile is Uo, then the initial output profile Y is 0 Y =GU 0 SS (6.48) +D, 0 ss therefore, Y in (6.45) can be expressed by ss Y = Y + G SU ss 0 SS (6.49) SS where SU„ = U - U ss (6.50) 0 In the CD-MPC objective function (6.1), there is a cost component for penalizing the actuator move between consecutive samplings, i.e., £JJ{k) Q, AU{k). For the static optimizer, this item is removed, as it has no effect on the steady-state. In other words, there is no limit on the range of actuator move other than the minimum and maximum set-point limits. Therefore, the static optimizer's objective function is represented as T 2 J = (Y - Y ) Q (Y - Y ) + Uj Q U , T s ss sp sl ss sp s s3 (6.51) ss where Q and Q have same dimensions as Qi and Q in (6.1). Replacing Y and U from (6.49) and (6.50), the objective function (6.51) can be expressed in the standard QP format, sl s3 ss 3 min 0.56Uj U 6U ss + <&5U aa (6.52) = Gj Q G + Q , (6.53) s ss ss SU aa where H ss = GlQ (Y s sl 0 sl ss s3 - Y ) +Q U . sp s3 (6.54) Q The solution of the static optimizer can be obtained through minimizing the above objective function (6.52) subject to the constraints (6.2)-(6.4). The dynamical constraints (6.5) have no meaning in this static optimizer. The solution to the static optimization problem (6.52)-(6.54) is not in general the same as the eventual steady-state of the dynamical MPC problem in (6.41)-(6.44). The next section will discuss how the static optimization weights Q and Q in (6.53)-(6.54) may be designed in order that the static optimization problem (6.52)-(6.54) will approximate the steady-state of the dynamical MPC problem (6.41)-(6.44). sl s3 Chapter 6. Steady-State Performance Prediction 6.3 107 Two Sensitivity Functions In this section we derive the matrices that describe the steady-state of the closed-loop output disturbance rejection transfer matrices (often referred to as sensitivity functions) for the dynamical MPC problem (6.41)-(6.44) and also for the static optimization (6.52)(6.54) in the absence of constraints (6.2)-(6.5). These matrices are required by the proposed technique which will (in Section 6.4) present a design for the static optimization weights Q i and Q in (6.53)-(6.54) that is intended to match the steady-state sensitivity functions of both problems. s 6.3.1 s3 Closed-loop Sensitivity Function for Unconstrained M P C The optimal solution for the unconstrained MPC QP problem (6.41) is AU(k) = -U- (f> = K Y (k) l r sp + K x {k) + K ^U(k x a u (6.55) - 1) where AU (k) AV(k)\ =i Hu U(k - = I Aui(fe), Au (fc), ••• , 2 Au (k) Nu 1) (6.56) (6.57) ( P ^ C ^ C P u (6.58) +Q + SlQaS^^J^Qi, 2 (6.59) K u-l -(PjC^QiCft + Q + SlQ S )- Sj„Q U . L 2 3 d u 3 p (6.60) From (6.23)-(6.25), after some algebra , x (k) can be represented by 1 a .\ , x (k) =T-Y{k)+Sa AU{k) (6.61) where V and © are matrixes in (G.3) and (G.8) in Appendix G respectively. Figure 6.1 shows the block diagram of an unconstrained MPC. Since the matrix K has no contributions to the closed-loop, when calculating the closed-loop system's sensitivity function to disturbances, it can be excluded from subsequent calculations. Substituting x (k) by (6.61) into (6.55) and ignoring K , the controller K(^) from the measurements r a r 1 Refer to Appendix G for details. Chapter 6. Steady-State Performance Prediction 108 D(k) Figure 6.1: The block diagram of an unconstrained MPC diagram for a multiple-array CD system. Y(k) to U(k) can be obtained through AU(k) U(k) K(z) I — K Q — ——j-Ku - 1 j X (6.62) K TY(k), x (6.63) AU(k) = K{z)Y{k), z-l z[(z - i)i -(zi)K e - K^pXr. x (6.64) The sensitivity function S(z) (i.e., the transfer matrix from D(k) to Y(k)) can be calculated by S(z) = (6.65) (I-G(z).K(z))-\ where G(z) is the plant model in (2.21). Because the control horizon H = 1 as mentioned in Section 6.1.1, §du = V — I and Q 3 = Q 3 in (6.60), at steady-state (i.e., z=l), the controller K{z) reduces to u P Kd = — (K _i) K T 1 u = If L(j) = [0, i ] , j = T X -Q^p^c qyCP T. T x (6.66) ,N (i.e., the default observer gain in Figure 6.1), at steadyU Chapter 6. Steady-State Performance Prediction 109 state, the term CP T in (6.66) is X CP T=[ I, I, .... X / ]. (6.67) r Finally, K in (6.66) can be further simplified as d (6.68) Kd - -Q G „ Qi 3 s m where sumYN u G (6.69) s G, yl ' "' ] C"aij A j B ij G t\f l\[ SUm y u H, Gsumij ^ ai fc=i = Tij a Hp—Tdij (l-oy) k (6.70) G ij, ss S ^ ' '---> v> fe=i i=i 0 1 =1 i N 3 = h---,N . u (6.71) Refer to Appendix H for the derivation of r from (2.20)-(2.25), and (6.8)-(6.13). The closed-loop system's sensitivity function S(z) in (6.65) at steady-state is tj S = (/-G(i)-iir )- = (/-G . fir )1 d d M J (6.72) 1 d (6.73) = (/ + G -Q 7 -GL Q )- . 1 ss 1 3 m 1 6.3.2 Sensitivity Function for the Static Optimizer Next we derive the analogous sensitivity function for the static optimizer in (6.52)-(6.54). Ignoring Y as we did in the above closed-loop sensitivity function, the optimal solution for (6.52) is SP 5U SS = -FL7> S -(Gj Q G s -Uo sl + Qss)- (GlQ Y + Q U) l ss SL - (Gf Q G s sl ss 0 + S3 Q )-\GlQ D ). s3 sl ss 0 (6.74) Then, we can find the relationship between Y and D through substituting (6.74) and (6.48) into (6.49), SS Y SS — YQ + G SU SS SS SS — S D , SS SS (6.75) Chapter 6. Steady-State Performance Prediction S = /-G (Gf Q G ss ss = K s sl 110 + Q, )- (Gf Q, ) (6.76) 1 ss 3 s 1 (I-G -K )~\ ss (6.77) ss = -Q;iGlQ aa (6.78) sU where S is the sensitivity function for the static optimizer and K is the controller for the static optimizer. Note that we used the matrix inversion lemma to obtain (6.77) from (6.76) [75]. ss 6.4 aa Design for Matching Optimization Problems In this section we propose a numerical optimization technique to design the static optimization weights Q i and Q in (6.53)-(6.54) such that the static sensitivity matrix S in (6.76)-(6.77) approximates the steady-state of the 'true' closed-loop sensitivity function Sa in (6.72)-(6.73). s 6.4.1 ss s3 Approximate Solution for the Static Optimizer For each sub-plant G in (6.47) with i = 1, • • • , N , j = 1, • • • , N , there is a positive scalar in (6.71) such that G = rG in (6.70). Therefore, the relationship ssij y sumij between G sum i:j u ssij and G is ss ruG ssU riN G u sslN u (6.79) rNylG Nyi SS ••• TNyNuGssNyNu (6.80) ROG rI n ••• r I 1Nu (6.81) R r J Ny ••• r I NyNu where '©' denotes the operation of element-by-element multiplication of two matrices. Because of the block diagonal structure of Q in (6.6), the difference between the inverses of the sensitivity functions Sd in (6.73) and S in (6.77) can be written as 3 aa Sj - S- = G (R l 1 T SS ©Q i ^ G L Q i - Q2G£Q.i). (6.82) The goal is then to design symmetric positive definite Q i and Q such that this difference s s3 Chapter 6. Steady-State Performance Prediction 111 is "small". To simplify this problem, we restrict Chi Q i • E, Qa3 Qs-F, (6.83) (6.84) where E and F are diagonal matrices E = V (eil, e I, ••• , e l), 2 with scalar e< and f F = V {fj, f I, ••• , JNJ) Ny 2 (6.85) Then (6.82) may be simplified as r & ~ Sss = Gss(R QZ- T (6.86) F- ZE), l where (6.87) An upper bound on the size of the difference in (6.86) can be obtained by noting II^T - S~ \\ 1 < \\G \\ l F SS • \\R QZ- (6.88) F- ZE\\ , T F l F where '|| • ||^' denotes the Frobenius matrix norm [75], which is defined as (6.89) where pij indicates the element of P in the i row and j column, and Ok(P) denotes the k singular value. Now we consider optimizing the second term on the RHS of (6.88) via the simplification th t h th \R 0 Z - F~ ZE\\ T = \\R* QZ - l F Fj Z E \\ , l f f f F (6.90) where \Zn\ IN J V Zt F = ||^N„I||F rn (6.91) ••• II-^N„N||F V ••• r lNu (6.92) Ri fNyi • • • r N N u Chapter 6. Steady-State Performance Prediction 112 E } = diag(e e , ••• ,e ), Ny (6.93) F F = diag(f f ,--- JNJ- (6.94) u 2 u 2 The matrices on the right hand side (RHS) of (6.90) are significantly smaller than those on the left hand side (LHS). For instance, Section 6.5.1 presents an example where R in (6.81) has dimensions 264 x 153 and the corresponding R in (6.92) is only a 1 x 2 matrix. The optimization problem suggested by (6.90) has the structure f min (6.95) ||*-VW||, F diagonal V,W then, setting \I> = i?J © Z and <i> = Zf, it can be solved for Ff = V~ and E = W using (for example) the iterative technique developed in Appendix I. The prediction performance index may be evaluated by the following: l f f \S -S. ss\\F \\Sd\ (6.96) ci where S i can be calculated from (6.73), and S can be obtained by (6.76) after Q in (6.83) and Q in (6.84) are solved using (6.91)-(6.95). c ss sl s 3 6.4.2 Special Cases In general, the performance index in the optimization problem in (6.90) cannot be made exactly equal to zero for diagonal matrices Ff and Ef. In order to achieve perfect matching, the matrix Rf in (6.92) must have a certain structure, such that fr ei l / R] _ 1 2 ei fy e l 2 fe l 2 fi e l f CN Ny L 2 2 Y (6.97) I JN U then a perfect solution for the minimization problem (6.90) can be found, and \\Rj 0 Zf — Fj ZfEf\\ = \\S~ — S~}\\F = 0. In other words, the static optimizer can perfectly predict the closed-loop system's steady-state performance. However, the elements of Rf in (6.71) are a function of the identified dynamics of the process and do not follow such a predictable structure. Thus the condition (6.97) is not satisfied in general, but two special cases occur frequently in practice. If the number of measured profiles N = l'in (2.21) and (6.85), i.e., there is only one controlled property, then we can get the perfect solution for (6.90) with Ef = 1 and L l F T y Chapter 6. Steady-State Performance Prediction F 1 = diag(r , • • • ,n ), f n 113 leading to Nu Q.n = Qu, Q*j = ^-Qz,, (6.98) j = l,---,N . (6.99) u Similarly, if the number of actuator arrays N = 1 in (2.21) and (6.85), i.e., there is only one actuator array, the perfect solution can be obtained with E = diag(r ,--- > w i) and Ff = 1, then from (6.83)-(6.84) this results in u r f Q.u = rnQu, n i = 1, • • • ,N , y (6.100) y Qssi = Qz\. (6.101) The problem and proposed solution above had the goal of designing the optimization weights Q i and Q in (6.51) and (6.82) such that the static optimization problem in (6.52) would match as best as possible the steady-state of the MPC control in (6.1), in the absence of actuator constraints. In order to be useful in practice, however the constraints (6.2)-(6.5) must be reintroduced into the problem. While the dynamical rate constraints (6.5) have no meaning for the static optimization in (6.52), the actuator position constraints (6.2)-(6.3) are imposed on the static optimizer as well as the MPC controller (6.1). While the matching was performed by comparing the unconstrained problems, the following section contains simulations that compare the results achieved by the MPC controller (6.1) with those obtained from the static optimizer in (6.52) with the weighting matrices Q and Q as designed using (6.91)-(6.95) subject to (6.2)-(6.3). s s3 si 6.5 s 3 Examples In numerous simulations, the above solutions have been found to be valid for the constrained CD-MPC. We present two representative examples, one for the perfect solution and one for the approximate solution. 6.5.1 T w o A c t u a t o r A r r a y s a n d One Controlled P r o p e r t y The first example is one of the most common multi-array CD processes and consists of two actuator arrays (N = 2) controlling a single measured profile (N = 1). The first is an array of n\ = 65 individual steam shower actuators (steam box) used to dry the sheet, and the. second consists of an array of n = 88 individual water spray actuators (rewet shower) that are used to rewet and overdried sections of the sheet. The moisture u y 2 Chapter 6. Steady-State Performance Prediction 114 profile is measured by a scanning sensor at m = 264 locations across the paper sheet. The parameters describing the process dynamics in (2.25) are a = [0.7515,0.9718], T = [2,1]. (6.102) d The spatial response matrices dj in (2.21) were identified by a commercial software [35]. The lower and upper limits U (j) = 0%, U (j) = 100% in (6.2), the bound for bending limits M (j) = 7.5%U (j) in (6.3), and the rate constraints AUmaxU) = 5%i7 (j) in (6.5), all with j = 1,2. The proposed design for the static optimizer does not depend on how the MPC in (6.1) dynamically tuned, but it does need Q and Q i in (6.1) as inputs which are shown in (6.83)-(6.84). Through trial and error, the MPC controller (6.1) was tuned for acceptable performance with the following parameters: horizons H — 10 and H = 1, the weighting matrices Q i = 0.337™, Q = 5 x 10" • V(I ,I ), and Q = 10~ • V(I ,I ). Due to strict rate constraints (5% of the upper limits), the MPC controller was required to be run for 150 sampling times until when the closed-loop system reached the steady-state. Because of only one controlled property, as discussed in Section 6.4.2, the perfect solution can be achieved from (6.71), (6.98), and (6.99). The scalars r^ in (6.92), Ef and min max max ma3; max 3 p u 5 2 4 ni n2 3 ni n2 F in (6.93)-(6.94) are f Rf = [5.2839,1.1770], E = 1, } Ff = diag (5.2839,1.1770), 1 (6.103) thus, the weight matrices for the static optimizer (6.52)-(6.54) are Qsii (6.104) = Qn, Qssi = Q3i/5.2839, Q s32 = Q /1.1770. 32 (6.105) In Figure 6.2a, the dashed line is the realistic disturbance which is generated from a software simulator of Honeywell, the flat solid line is the set-point profile, and the solid curved lines are the measurements predicted from the static optimizer and the outputs from running the closed-loop simulation to steady-state (150 sampling times). Figures 6.2b and 6.2c show the steam boxes.and rewet shower actuator profiles, where the solid lines are the actuator profiles predicted by the static optimizer and the dotted lines are the actuator profiles after running the closed-loop simulation (in fact, they are indistinguishable). Figure 6.2 shows that there is no difference in the measurements or actuator profiles between using the static optimizer and using the closed-loop simulator. However, the computational time for the static optimizer is only 1.4 CPU seconds versus the closedloop simulator's 137.3 CPU seconds (on a 700 MHz CPU, 768 MB RAM computer), the Chapter 6. Steady-State Performance Prediction 115 later number is the total time for running 150 sample times. Note that two steam box actuators (the last two) in Figure 6.2b hit the lower constraints U (j) = 0% in (6.2). This means the perfect matching for the closed-loop MPC simulation in this example can be reached through the static optimizer even with active constraints. min 10 20 30 40 50 60 70 80 Actuator number Figure 6.2: (a) The moisture profiles before control (dashed) and after control (solid: static optimizer; dotted: closed-loop simulation); (b) The steam box actuator profiles (solid: static optimizer; dotted: closed-loop simulation); (c) The rewet shower profiles (solid: static optimizer; dotted: closed-loop simulation). 6.5.2 Four A c t u a t o r A r r a y s and Three C o n t r o l l e d Properties This plant model comes from a Canadian working paper machine. This paper machine has four actuator arrays: slice lip ui, steam box u , rewet shower u , and induction heating u , and the controlled properties are basis weight j/i, moisture y , and caliper (thickness) y . The numbers of the actuator arrays are n\ — 57, n = 65, n = 88, n — 84 in (2.23), and the dimensions of the controlled properties are the same as m — 264 in (2.22). In this example, note that Gi hi (z),G shis(z),G hu(z) and G h (z) of (2.20) are zero, because, for example, steam box u does not affect the basis weight j/i. Other subplants Gij in (2.21) were identified by the commercial identification software [35]. The dynamical 2 3 4 2 3 2 2 2 l u 2 A 3 2A 24 Chapter 6. Steady-State Performance Prediction 116 parameters a and T in (2.25) are d a = Td = 0.2636 0 0 0 0.4029 0.6905 0.3858 0 0.8720 0.5908 0.2636 0.9174 (6.106) 2 0 00 2 1 00 2 1 02 (6.107) The lower and upper limits U (j) = 0%, U x(j) = 100% in (6.2), the bound for bending limits Mmaxij) = 7.5%t7 (j) in (6.3), and the rate constraints AU x(j) = 2.5%U (j) in (6.5), all with j = 1,2,3,4. Through trial and error, the parameters for the MPC controller in (6.1) was set up as: the horizons H = 10 and H = 1, and the weight matrices in (6.1) are min ma max ma p max u Qi = £>(1.65/,0.33/,0.33/), (6.108) Q (6.109) m m m = 5x 10- -7J?(/ ,/ ,/„ ,/ ), 4 2 Q3 ni = 10 • 4 n2 3 n4 T>(I ,I ,I ,I ). ni n2 m n4 (6.110) Due to tight rate constraints AU in (6.5), the closed-loop simulation needs to be run 200 sample times for reaching the steady-state. From the model and the parameters of the MPC controller, the matrix Zf in (6.91) can be calculated as max 794.9 0 0 0 129.3 4.3867 3312.5 0 543.8 469.8 785.7 1590.4 Z = f (6.111) From the dynamical parameters a in (6.106) and T in (6.107), and the horizons H = 10 and H = 1, the scalar element m of Rf in (6.92) is computed from (6.71), d p u R f 7.6421 0 0 0 7.3257 6.8488 9.3719 0 3.4653 7.5690 9.6420 2.4658 - (6.112) Using the iterative approach developed in Appendix I, the optimal solution for (6.90) (after 50 iterations) is: E f = dmp(16.1762,7.5071,7.7914), (6.113) Chapter 6. Steady-State Performance Prediction 117 Ff = $00(2.1173,1.0952,0.8014,3.1598). (6.114) Thus, for the static optimizer, the optimal choices are Q = Q± • E and Q = Q • F from (6.83)-(6.85). The prediction index n = 99.12% in (5.8) indicating the static optimizer is providing good prediction. This can be further verified in Figures 6.3-6.6, which show the measurements and actuator profiles predicted from the static optimizer and those from the closed-loop simulation. All the disturbance profiles in Figures 6.3-6.5 are generated from a software simulator of Honeywell. The computational time of the static optimizer is 32.8 CPU seconds versus the closed-loop simulator's 3283.5 CPU seconds (200 sampling times were required to reach steady-state). Note that in Figure 6.6, slice lip actuators ui and induction heating actuators u hit the position constraints (6.2), which mean that the static optimizer works well even with active constraints in this example. Table 6.1 shows the standard deviations of the measurement profiles and the actuator profiles. Although there are slight differences in the steam box profile Ui and rewet shower profile u between using the closed-loop simulation and using the static optimizer, the moisture profiles are almost same. The reason is that the optimized design focussed on the output profile and different input combinations could generate the same output. s l s 3 3 4 2 50 100 150 CD locations [measurement points] 200 250 Figure 6.3: The basis weight profiles for the example in Section 6.5.2: static optimizer (solid+asterisk), closed-loop simulation (solid+circle), disturbance (dashed line), set-points (flat solid line). Chapter 6. Steady-State Performance Prediction 50 100 118 150 200 250 CD locations [measurement points] Figure 6.4: The moisture profiles for the example in Section 6.5.2: static optimizer (solid+asterisk), closed-loop simulation (solid+circle), disturbance (dashed line), setpoints (flat solid line). 6.6 Summary In this chapter, we have proposed a static optimizer to accurately predict steady-state performance of large-scale constrained CD-MPC systems. The penalty functions of the optimizer can be obtained through minimization of the difference between the inverses of the sensitivity functions related to the optimizer and an unconstrained MPC closed-loop system. The most important contribution of this chapter is that this method can be used for constrained CD-MPC. For the paper machine CD-MPC, the closed-loop system's steady-state performance can be optimally predicted. In two special cases (only one controlled property or one actuator array in CD processes), the optimal solutions happen to be perfect. The simplification of the complex minimization problem related to the large dimensional sensitivity functions has been addressed through solving for the optimal ratios. Two industrial examples have demonstrated that the computational time using the static optimizer for predicting the steady-state performance is significantly less than that using the closed-loop system simulator, while the results are almost the same for the two methods. Chapter 6. Steady-State Performance Prediction 50 100 119 150 , CD locations [measurement points] 200 250 Figure 6.5: The caliper profiles for the example in Section 6.5.2: static optimizer (solid+asterisk), closed-loop simulation (solid+circle), disturbance (dashed line), setpoints (flat solid line). 5 10 100, 15 1 10 1001 251 , 20 1 30 35 1 20 , 25 30 , 1 40 1 40 , 45 1 50 55 1 50 , i i i i i i i L_ 10 20 30 40 50 60 70 80 Actuator numbers Figure 6.6: The actuator profiles for the example in Section 6.5.2 using the static optimizer (dashed plus point line) and closed-loop simulations (solid line). Chapter 6. Steady-State Performance Prediction 120 Table 6.1: Performance and computational time for the MPC using the closed-loop simulation and the static optimizer 1.334 1.828 1.961 0 0 0 0 closed-loop 0.8705 0.86 1.247 35.23 5.247 9.413 13.41 After Control static optimizer 0.8708 0.8622 1.249 35.3 5.749 9.888 13.23 - 3253.8 32.8 Before Control Performance (two times the standard deviation of the profile) K y i Ky 2 Ky 3 Kill ^ U 2 K « 3 Computational Time [CPU second] Change 0.0345% 0.26% 0.16% 0.2% 9.57% 5.05% -1.34% -98.99% Chapter 7 Conclusions and Future Work 7.1 Summary This work has concentrated on solving the problems of industrial model predictive control (MPC) for paper machine cross-directional (CD) processes, which are known as largescale, two-dimensional, and ill-conditioned. Application of M P C to such processes has been a challenging task due to their large-scale characteristic plus the inevitable severe model uncertainty in the models. This work consists of three parts: model reduction in Chapter 4, analysis of robustness and performance in the two-dimensional frequency domain in Chapters 3 and 5, and performance prediction in Chapter 6. Wavelet packets can be used as one orthogonal basis function for paper machine C D process model reduction. This reduction is closely related to the cut-off frequency of the spatial response and the non-robustly controllable singular vectors. The indices of model reduction performance including the distortion index and the sparsity index have been proposed for model reduction evaluation. There is trade-off among these conflicting indices. Different wavelets and different wavelet filter lengths affect the reduction performance. Longer" wavelet filters can reduce the distortions but will increase the sparsity of the reduced model. Symlets with filter length between 8 ~ 16 are good choices for majority of paper machine cross-directional processes. This reduction method has been tested under numerous C D process models which simulate real working paper machines, and the results demonstrated its effectiveness. The constraints of the C D - M P C can be exactly represented in the reduced domain by using the inverse wavelet transform. A cross-directional model predictive controller based on the reduced model has been designed. This technique can reduce on-line computational time by about 50% while maintaining close to optimal control performance. It also brings some degree of robustness and will not result in actuator picketing as it does not attempt control of the non-robustly controllable directions and does not have high spatial frequencies in its actuator profile. This wavelet-based model predictive controller can be used for multiple-array systems, whereas some other reduction methods such as the singular value decomposition method can not be used. A two-dimensional frequency analysis method has been presented in Chapters 3 and 5 for cross-directional model predictive control. A rectangular circulant matrix (RCM) was defined and its properties was discussed in detail. If the column vectors of an R C M 121 Chapter 7. Conclusions and Future Work 122 are bandlimited to the spatial Nyquist frequency of inputs, its singular values can be easily obtained by pre- and post-multiplication of the complex Fourier matrices. The R C M structure can be preserved for the closed-loop transfer matrices in a linear feedback system if the open-loop plant and the controller are each an R C M and both of their column vectors are bandlimited to the spatial Nyquist frequency of inputs. The process model, the additive structured model uncertainty and the closed-loop transfer functions has been addressed in Chapters 2 and 3. They are "close" to being R C M s and can thus approximately be decoupled in the two-dimensional frequency domain. This allowed these large-scale M I M O transfer matrices to be analyzed in terms of a family of simple SISO transfer functions, one for each spatial frequency. The two-dimensional frequency bandwidth, the performance index, and the robust stability margin index were presented as intuitive measures of closed-loop performance and robustness. Analysis of the realistic C D process model uncertainty in the spatial frequency domain can give a better physical insight into the uncertainty problem. The balance between robustness and performance of the C D - M P C can be achieved by tuning the weighting matrices Q and Q . Roughly speaking, in the presented examples, 2 3 the weight Q affects the closed-loop dynamical behavior while the weight Q 2 3 determines the spatial performance. The examples in Chapter 5 have shown that this two-dimensional frequency analysis method gives consistent criteria for selecting the weights of the C D M P C using the guidelines from tradional robust control. The presented examples in Chapter 5 also demonstrated that proper selection of the weighting matrix Q ' s structure 3 and careful consideration of the structured model uncertainty will remove unnecessary conservativeness from the traditional design based on diagonal weighting matrix Q 3 and unstructured model uncertainties. This two-dimensional frequency analysis method was extended to multiple-array systems and demonstrated by the example in Section 5.5.3 of Chapter 5 for its effectiveness. A static optimizer has been presented in Chapter 6 for solving the computational time problem of predicting the steady-steady performance of the very large-scale constrained model predictive controller. The static optimizer which uses one-step optimization is much faster than the traditional closed-loop simulator. The design of the optimizing weighting parameters of this optimizer is solved by minimizing the difference between two sensitivity functions, one for the optimizer and another is for the closed-loop M P C simulator. Two examples in Chapter 6 demonstrated that the static optimizer is significantly more efficient (up to two orders of magnitude) than the conventional closed-loop simulation method while accurately predicting the steady-state performance of constrained C D - M P C . Chapter 7. Conclusions and Future Work 7.2 123 Future Work It may be necessary to investigate the parameter uncertainty limits presented in Section 2.1.1 of the process models in real working paper machines. A reasonable practical model uncertainty range is very helpful to design more effective MPC. After the model uncertainty range is known, the structures of the weighting matrix Q should be investigated for better performance and robustness, such as other high-pass spatial filters rather than the bending moment matrix. Although this work has given theoretical guidance for tuning the CD-MPC, it may be necessary to develop an automatic tuning tool of the CD-MPC for non-expert users. This tool should include automatic model reduction, optimal parameters tuning and accurate steady-state performance prediction. Such a tool will bring the maximum benefit for CD control systems. Field trials will be necessary for developing the automatic tool. It is known that weakly controllable directions exist between different actuator arrays. For example, steam boxes and rewet shower actuator may easily fight each other for controlling the moisture profile in industrial CD control. It is likely that the CDMPC objective function will,require modification for penalizing these weakly controllable directions. 3 Appendix A Closed-Loop Transfer Functions of Unconstrained M P C In this appendix, we derive the static matrix K in Figure 3.1 [54]. We have to use equations (2.16) and (2.17) for the future output predictions in (2.26). The states x (k) in (2.16) are not measurable and have to be estimated. We denote the estimated states as x (k) as generated by the observer in (3.2). From (2.16), we can predict the future states: mpc a a X(fc) = P x (k)+P AV{k), x a (A.l) u where x {k + l) Au(k) a X(fc) = x {k + 2) a , x (k + H ) a a Au{k + j p A Au(fc + 1) AU(fc) " B ' a Al > AB a P 1 a (A.2) H -1) U 0 ••• 0 B ••• 0 a (A.3) u A Hp A^~ B Aa ~ B l -rla P a 2 a ••• Aa ~ B P Hu a The current estimated states x (k) can be estimated from the previous states and the current measurement as in (3.2). The estimated outputs y(k) can be obtained by a (A.4) y(k) = C x {k). a a Equations (3.2) and (A.4) can be realized as shown in Figure 3.1. Therefore, the output predictions Y(k) can be obtained by y(k + i) Y(k) m+2) = C X(fc), 0 y(k + Hp) 124 (A.5) A. Closed-Loop Transfer Functions of Unconstrained MPC 125 where C = V(C , ••• ,C ), 'P' denotes the operation of constructing a block diagonal a a a matrix. The cost function (2.26) can be rewritten with (A.1)-(A.5) as (Y(fc) - Y (^)) Q (Y(fc) - Y,(fe)) + AV(k) q AU(k) + V(k) q V(k), T T sp T 2 1 3 (A.6) where y (k +1) y (k + 2) sp sp Y (k) (A.7) sp y (k + H ) sp p u{k) u{k + 1) U(fc) = U u{k- 1) + p u{k + H U S AV(k), du (A.8) ? Hu, (A.9) 1) I I >H, u Sdu — I I I h = I I ®i=T>(Qi,--- W i . - . Q i ) , ,Qi), 2 = 2,3. (A.10) With some algebra, the MPC problem (A.6) can be transferred to the standard quadratic programming (QP) format: min 0.5AU(fc) $AU(fc) + (f> AV(k), r (A.ll) T where $ = Pjc q c p T a 1 a u + q + sJuq3sdu, <t> = - p ^ q r Y ^ (A.12) 2 + p^qiCaP^i^ + slqsUpuik-i). (A.13) Assuming that y {k) is a constant target profile, then 4> can be rewritten as sp <t> = -Pu^QiEy (k) ap + PfCZ^CaPtXaik) + Sj q U u(k u 3 P - 1), (A.14) A. Closed-Loop Transfer Functions of Unconstrained MPC 126 where E = [I, I, The optimal solution for the unconstrained MPC QP problem (A.ll) is AU(fc) <f> = y (k) K l mpc x (k) sp K •mpc (A.15) u(k -1) a (A.16) K«-i where (PjCQ C P + Q + 1 a u +q + -(p^q.CaPu (A-17) SlQsS^P^q.E, 2 2 (A.18) sl® s )- sl® u , 3 l du 3 (A.19) p where each matrix K , Ku_i is static. But (A.15) contains control moves up to H — 1 time steps into the future. In practice MPC implements only first move, recomputing (A.ll) at each time k, thus, only first n of AU(A;) (i.e., Au (k)) are used for implementation, Au (k) can be obtained by, x u opt opt Au {k) = TAV{k), Therefore, in Figure 3.1, K ,K ,K mpc r = TKmpc, Kmpc x u r r (A.20) and K _i can be obtained by, = rK , K ,0'. T = [/, 0, opt K x — riKj, K _i (A.21) rK _i. = u u Now, Ky(z) and K in Figure 3.2 are derived as follows: substituting A ,B ,C (2.19) into (3.2), x {k) can be simplified-as a y a a in a Ax{k) x {k) a L-y(k)+6Au(k), Cx{k) (Zl xn 0 nx ~ x A)~ B l 0. n xn (A.22) x where L — [0,1] . Then, substituting x (k) in (A.22) and u{k — 1) by Au(k)/(z — 1) into (3.1), T a Au(k) u{k) I-K 9 r —Au(k) 1 -K ^ u \ {K -y + K -L- y(k)), r sp x = Ki(z) (K • y + K • y(k)) r sp y (A.23) (A.24) A. Closed-Loop Transfer Functions of Unconstrained MPC z-1 K„ = K -L. x l-KJ- •K,u-l 127 (A.25) (A.26) Appendix B Wavelet Packet Analysis Tree The next page includes the 5-stage standard wavelet packet analysis tree and gives the frequency sequence of the coefficients. ho(—n) and hi{—n) are the low-pass and high-pass filter in (4.5) and (4.6) respectively. At each scale level, the coefficients have a convolution with the corresponding filter and then are downsampled by 2 for getting the next scale level's coefficients. Note that for each scale level, the corresponding wavelet packets are arranged by natural order, but not in frequency order. Details may be found in the Wavelet Toolbox [57]. 128 Appendix C Proof of RCM's Properties Cl Proof of R C M Property 1 Assuming m = n -n and the i and (i+l) column vectors of Ni and N are a*, a ,b~i, b respectively, then the i and (i+l) column vectors of N are ~c\ = ai+h, c i = a i+b i. Because N and N are RCMs, then from the RCM definition, then a , b are equal to a,i and bi circularly down shifted n times respectively, therefore, c i is equal to c\ circularly down shifted n times. Thus N is an RCM by Definition 1. th a th th r 2 th 3 i + 2 i+1 a i+l i+ i+l i + a C.2 i+l 3 Proof of R C M Property 3 Please refer to [17, 41] for the proof. C.3 Proof of R C M Property 4 We have an RCM N G C , here assuming m = 2n for simplicity without losing generality. First, we prove that N = F NF^ is a matrix of diagonal blocks. Then, we prove that M = F*NF is an RCM. (1) "If proof. The complex Fourier matrix F and F% can be written as mxn m n m 1 1 1 1 1 Q- 1 3 J- m. V2n 1 Q 2 1 I Q Q- Q- 2 e" 3 1 3 Q- 1 Q 2 Z -{2n-l) • Q -9 2 Q Q Q -(2n-l) 3(2n-l) 2 Q 130 e _(„_l)(2„-l) 3(2n-l) 2 g-(2n-l)(n-l) 0-(2n-2)(n-l) 2 n - l 2 2 n - l 2 g _ffn-i)ffn-l) (Cl) i+ C. Proof of RCM's Properties 1 131 1 1 1 Q 1 Q Q 1 Q Q _ 1 1 1 Q 2 2 • 4 V Q • • Q Q • g3(n-l) g4(n-l) . 8 Q 6 n-1 1 Q 4 3 Q 1 9 2{n-l) 3(n-l) Q Q N . g 2 ^ (C.2) (n-l)(n-l) where g = e / , H denotes conjugate transpose. Note that g = 1, Q = g~ , F 2m n n H l M F j f are unitary matrices, i.e., F F" and = I, F F^ = L M N A n m-by-n R C M with m = 2n in Definition 1 has the following general form, TV = ty t t$ tin tl t tln-l ^1 t 3 tln-3 t^n t 2 tln-i t t t tl tl ti t% tin tln-2 3 3 5 tln-l tln-2 4 7 (C.3) Now introduce the 2n-by-2n permutation matrix 10 W 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 I>2nx2n (C.4) C. Proof of RCM's Properties 132 and form the product £ 3 £5 ti £2n-l £2n-3 fan-l ti £3 S £l N = W-N = £3 tin £5 £7 tl £4 tin £ n-2 2 ti £4 t 2 ••• £2n-2 ••• tm-4 te (C.5) R tin Note that the new matrix N £ C is formed by stacking two square circulant matrices (SCMs) S and R e C , and WW = W W = I. Similarly, we reorganize F in (C.l) as F = F -W in (C.6) in order to let F NF" — 2nxn nxn T T T m (F • W ) • (W • N)F£ = m m m 1 1 1 1 e 1 Q -(n-l)(n-l) -(n-1) V2n y-n(n— -(n+1) 1 75 m F NF», T 1 m g,-(2n-l) F F F -F . . . g-(n-l)(2n-l) Q e 2 Q e e 2n-l 3 2 n-l £ 1) .,-n(ra+l) 1 1 n Q 2 3(n-l) (n-l)(2n-l) 2 £ 2 3n n(2n-l) e3(n+l) 2 2 n+l £ 2 (n + l ) ( 2 n - l ) 2 £ 3(2n-l) 2n-l e 2 e 2 n n 2 2 (2n-l)(2n-l) 2 (C.6) n n where AF , F A = diag(l,g n n *,Q \Q ,g 2 ) . (C.7) Please note that 0 = 1 and £ 2 = — 1 have been used in the above simplification. N can be obtained by n N F NF? = m 1 V2 F F ± F NF? m F n n -L 5 n -Fn 1 F SF" n F SF" n + kF RF» n AF RF£ n (CS) C. Proof of RCM's Properties 133 Since S and R are square circulant matrices, then F SF£ and F RF" are known to be diagonal [17], and since A is diagonal in (C.7), then N in (C.8) is formed by two diagonal matrices n N n Ni N = (C.9) 2 1 = ^(F SF» i N + AF RF»), n N = -j=(F SF?-AF RF»). n 2 n (C.10) n (2) "Only if proof. From (Cl), the conjugate transpose of thera-by-m.Fourier matrix with m = 2n F% e C can be written as 2nx2n F" = fi = (Cll) 2n 1 ^ 2 \/2n Consider a matrix N £ C 2 n x n 2(i-l) 2 t-l Q • • • Q (i-l)(2n-l) 2 i = l,---,2n. (C12) of stacked diagonal form, Zi 0 ••• 0 0 l '•• 2 : : 0 '•• '•• 0 ••• 0 l ln+1, 0 • • • 0 N = l +2 0 ' '• n : 0 (C.13) n '•• '•• 0 ••• 0 l 2n _ Now, post-multiplying with the n-by-n Fourier matrix, h k JVxF„ = 1_ -= hQ- 1 ^n l>nQ -(n-1) i 0-(«-D(n-O n Vi 'n+1 (•n+l ln+2 ln+2Q~ kn 2 ln 2Q- ~ {n l) + l knQ- ~ {n v l) knQ- - ~ {n l){n Y) (C.14) C. Proof of RCM's Properties 134 Then, pre-multiplying with the m-by-ra Fourier matrix, M = F^NF . Now an expression for the matrix elements M(i,j) and M(i + 2,j + 1) can be obtained by n = h+ Q^{hQ- - ) + ••• + Q ^ {lnQ- - - ) {j Q l) i l +l + £ 2 (Jn+2£ ° 2 n »-2; + l = /i + Q VQ 2 - £ *n+l 2 5 i+2 1" £ 2 - ^ U l - P (C.15) fen 2 Q * * ^ ^ ^ ) (i+l)(2n-l) (/n+2^ ) + ••• +Q 2 « 4 ni } ( 2 n i - 2 n j + 2 j - t — 1) (i+l)(n+l) /n+l + £ 2 ) K j+ + ••• + {i+l)n t-2,, + 1 [knQ " 2 =f -V i j = Zl + £ + l) Z + •••+ Q 2 = h + Q^(hQ- ) 2 l){n n ( n t - 2 j - f t+1) M(i + 2,j + l) + £ 2 {j 0 + • • •+ Q (i-2j+l)(n-l) ZH 2 ni 11 J J ( ') (i-2j+ l)(n-l) In 2 (ni-2j+i+l) 2 n (feng 2 Zn+2 (2nt-2nj+2j-t-l) + •••+£ Z . 2 2n (C16) By inspection, M(i + 2, j + 1) = M(z, j) for i = 1, • • • ,2n - 2, j = 1, • • • , n - 1, (note that when i = 2n - l,2n, j = 1, • • • ,n - 1, M(z + 2 - 2n, j + 1) = M(i,j), and when i = 1, • • • , 2n - 2, j = n, M(« + 2, j + 1 - n) = M(i,j), and when i = 2n-l, 2n, j = n, M(i + 2 — 2ri,j .+ 1 n) = M(i,j), all of these cases can be proved in a similar fashion), indicating that column j + 1 — column j circularly down shifted n = 2 times. — a C.4 Proof of R C M Property 2 For the case n = ni, form the product 2 iV 4 = F NTN F£ ni a = {F N?F% )(F N FZ), ni 1 mi 2 (C.17) then both F^NfF^ and F A^ F^ are stacked diagonal matrices as in (3.9) and it is straightforward to verify iV is diagonal, and therefore N is an R C M by Property 4. The case for n = m-i is proved similarly. mi 2 4 4 2 C.5 Proof of RCM's Property 5 Assuming m = n -n. Since the Fourier matrix in (3.8) is a unitary matrix, then, a cr,(JV) = a^NF?) = aj(N). (C.18) C. Proof of RCM's Properties 135 By the definition of singular value decomposition [75], a {N) = y/\i(N"N), (C.19) t and since N is a matrix of diagonal blocks (by Property 4), then equation (3.10) can be easily obtained. The case for n = n • m may be proved in a similar fashion. a C.6 Proof of RCM's Property 6 We assume that m = 2n and n = 2k + 1 for simplicity. From (C.8), S. and R are SCMs, then, S = F SF" and R = F RF^ are diagonal matrices and n n S(j+ + = S(n-j + l,n-j + l) , R{j +l,j + l) = R( -j + l,n-j + l) , (C.20) (C.21) H H n where j = 1, • • • , k [17]. Therefore, N(J + 1, j + 1) and iV(2n - j + l,n - j + 1) are N(j + 1J + 1) = 5O" + U N(2n- j + l,n- (C.22) + l) + 0-*i*(j + l , j + l), j + 1) = S(n - j + l,n - j + 1) - g~^R{n = S(j + l,j + l) - j + l,n - j + 1) + (o-iR(j + l,j + l)) H H = N(j + l,j + l) . (C.23) H Note that g N(2n-3 2 = + 1,71-3 Because n (j + 1) — N(j + l,j + 1) and n {2n — j + 1) = + 1), (3.12) is confirmed. —1. d d Appendix D Proof of Theorems D.1 Proof of Theorem 1 Here we assume that m = 2n for simplicity without losing generality. First, we construct a SCM N e C from the RCM N in (C.5) by inserting the column vectors, 2nx2n s h tU t • t -lt t t • • t t- -l k t • t n-2 t -l t tl •• t t--3 2 tn 2 hn-l 3 2n 2 3 t2n t n-2 2 2n 2 t t 2 2n 2n h 2 2n 2 •^ 2 n - 3 u u to • • u u u •• t 3 2n 2n 2n (D.1) A t tl 2 2 Similarly, we insert zeros in every two elements of u, ui = u{l), 0, u(2), 0, « ( 3 ) , 0, (D.2) u(n), 0 It is apparent that N • u — N • u, and we may write s y y = (F N F%) m s • (F • ) = N • (F W W ) =N • T m Ul s m Ul s (F ), mUl (D.3) where F and W are from (C.6) and (C.4) respectively, m Ui = W • U\ u (D.4) 0„ substituting F in (C.6) into (D.3), m F -u u F -U u n n (D.5) Next, we let n € C represent the diagonal elements of N in (D.3). Because N in (D.1) is a SCM, n is equal to the Fourier transform of N 's column vectors rij in (3.21) 2n s s 3 s 136 s D. Proof of Theorems 137 [77]. Let v (i) represent the frequencies of n,-, then, r v {i) r e [0,1/j,], i = l,---,j , (D.6) b where z/ and % are from (3.23), Q (i) denote the normalized frequencies as Q (i) in (3.17). Because we assume m = 2n, then X = 2X in (3.19). Then, because 0 < v (i) <v < v^, from (3.18) and (3.19), we can get: r B y u y r B 0<a(i)<|- (D-8) Because of symmetric characteristic of n- in (3.21) [62], 3 3TT — < (D.9) n (i) < 2TT. r We can obtain i satisfying (D.8)-(D.9) using (3.17), ' . = 2 f 1."'2, • • • Jloor(n/2) + 1, \ floor(n/2) + 2 + n, floor (n/2) + 3 + n, • • • , 2n, for A; = 1 forfc= n = 2 ^ a Then, given S = [0, • • • , 0,1,0, • • • , 0] with u(j) = 1, from (D.5), r for i in (D.10) ' . otherwise ^ f y ={ ^ [ 0, , (D.ll) y Further, from (3.9) and (3.14), because of only a single frequency input on u and the diagonal characteristic in Nj in (3.9), then (3.24) is confirmed. D.2 Proof of Theorem 2 1. The matrices A,B,C in (2.14) and (2.15) can be generated through the following method: first, assuming the inputs w(k) = G -u(k), transforming y(z) = h(z)Iw(k) (h(z) in (2.2)) into state-space, that is, 0 x (k + l) = A x (k) h h y(k) h + B w(k), (D.12) h = C x (k) + d(k). h (D.13) h Because h(z) is a first-order plus dead time (p in (2.2)) system, ^2^ ChA{B or J2?=o ChA\(zI — A)~ Bh with H > p is a diagonal matrix, in other words, they are 0 l p h D. Proof of Theorems SCMs. Second, 138 C = C ,A and B = B G . = A, h h h 0 2. From Condition (1), Go is an RCM, then, because h(z) in (2.2) is a scalarfora particular z, the open-loop plant model G(z) in (2.1) is an RCM. 3. For the MPC controller, the control horizon H = 1 for simplicity, the prediction horizon H > p,p is time delay in (2.2), Q = qil , Q = <?24xn and Q = q I„ , where gi, q , and q are scalar numbers. u p x 2 mxm 2 3 3 3 From equation (3.5), the controller's transfer matrix is given by K(z) = substituting K and K „i x (i-KJV 2 - 1 } Z - l (D.14) K L, X in (A.21) into (D.14), u K(z) = 1 z-l (D.15) z-l where L and 9 are given in (A.22), and *o PuCl®iC P = a (D.16) (D.17) + Q + Q, u 2 3 Then, substituting P in (A.3), C in (A.5), Q i in (A.10), and P in (A.3) into (D.16) and (D.17), and using (2.19) and (A.22), u x a CB a CB a a CAB a a a CAB a a a a qI + 2 C Aa B p CA l a a a T CB CAB + CB r P Hp-l _1 a CB CAB EST CA B i 3 B p a + CB + qI + g J 2 E,= o g7 1 3 CA'B -j ^ ( C ^ B ^ C A ^ ) +ft/+ q I 3 qi Hp-l j ^(a4B Go) (C,Ai B,G ) + 9i r ft j=0 i=0 l 0 + fe/, (D.18) xn D. Proof of Theorems 139 T CB a ' a CA a CAB a a ' a CaAl a •qj- C Aa C A?> B p l a a a T CB CAB 1 CB + CA + 2 I CA (D.19) CA B l T -i CB CAB I CA + Eflo r r I CB Hp-1 j 3=0 i=0 ft E E ^ ^ ) ' (D.20) j=0 i=0 n #•0 T l {CA CAB + CB = ft 2 i 1 = l Y% CA {zI-A)- B E S T CA'B Hp-l CA(zI-A)~ B + CA)(zI -A)~ B r l x j Y,( ) ( CAiB T ( CAi+1 - )~ ) zI A lB j=0 i=0 = 9i E E ^ ^ ^ G o ) 7 ^ ^ 1 ^ - ^ ) - 1 ^ ^ ) . (D.21) i=0 j=0 As mentioned above, ChA B and ChA\(zI — A)~ Bh with i = 0, • • • , H is a diagonal matrix (Ah, B , and Ch in (D.12) and (D.13) ). Since G is an RCM, GQGQ is a SCM based on Property 2. Because q2l,q I are SCMs, <& in (D.18) and * • 6 in (D.21) are SCMs, and * • L in (D.20) is an RCM, then * + W + ^Q is a SCM from Property 1. The inverse of a SCM (i.e., $ + + ^ r Q ) is a SCM from the RCM's property 3. The transpose of \& • L is an RCM. According to the RCM property 2, the product of one RCM and one SCM is an RCM. Therefore, the controller K(z) G C in (D.14) is an RCM. l Y h h p h 0 0 3 3 0 0 n x m 3 D. Proof of Theorems 140 Because the column vectors 7/- of G satisfy the condition of Theorem 1 (i.e., condition 2 of Theorem 2), only the upper half diagonal elements of Ni and the lower half diagonal elements of N of G(z) and K(z) in (3.9) are non-zeros. Then, (G(z)-K(z)) is a square diagonal matrix, that is, (G(z) • K(z)) is a SCM. Finally, based on the RCM's property 3, T {z) = [I- G(z) • K(z)]~ is a SCM. Then, from the property 0 na l yd 2, T (z) = K(z) • T (z) is an RCM. ud yd Appendix E Calculation of the Eigenvalues of T^A In this Appendix, we derive an analytical expression for the eigenvalue problem in (5.18). In this derivation, we assume without loss of generality that the first actuator array's dimension n i is an even number and the second actuator array's dimension n > rt\. Denote the diagonal matrices M and M s diagonal elements by 2 a d {ai,---,a } = diag(M ), n] (E.1) a {di,-- - ,d } = diag(M ). n2 (E.2) d and the nonzero "diagonal" elements of M in (5.23) and M in (5.24) as b c {&!,••• A J = D I A G ( M ) , (E.3) {ci,-- - ,c } (E.4) 6 = DIAG(M ), ni C where D I A G is the generalized diagonal operation defined in (3.31). Then, the eigenvalues A of T A. in (5.18) can be calculated as follows: ud det(A7 - T„ A) D = det XI -M -M^ a (E.5) XI - M -Mr d = det(A/ - M ) • det ((A/ - M ) - {-M ){XI - M )- {-M )^j l a d C a b (E.6) where 'det' is the determinant and det(AJ - M ) is assumed to be nonzero. Here we use Schur's formula [75] for the determinant of a partitioned matrix. The second term in E.6 must be zero in order to calculate the eigenvalues. a det ((A/ - M ) - (-M )(XI - M )-\-M )) d C a =0 b A - (oi + di)A + (aidi - cJi) = 0/ (E.8) 2 j^ ~ (^m/2 + d )X 2 ni/2 _ + (a d ni/2 141 ni/2 (E.7) (E-9) _ -c b ) ni/2 ni/2 = 0, (E.10) E. Calculation of the Eigenvalues of T u d A 142 A = fi„/+i, (E.ll) 1 2 j (E.12) A = cZ _ / , n2 A 2 ni (E.13) 2 — (<2 /2+l + d _ /2+i)A ni n2 ni • + (fln /2+l^n -ni/2+l 1 2 ' A - (S + d )A + (a d n2 = 0, ni _ 2 ni ~ Cn^/2+\b /2+l) ni n2 (E.14) ( - c^Ai) = 0. E ' 1 5 ) (E.16) Note that each equation from (E.8) to (E.10) and from (E.14) to (E.16) has two eigenvalues (the total number is 2n\). For example, the solutions for (E.8) are Ai = (oi +di)±J(oi + dx) - Affltdy - cJi) 2 2 l2 • ( E - 1 7 ) T d A has a total of n\ + n eigenvalues. The spatial frequency analysis method allows the huge eigenvalue problem to be reduced to a family of small eigenvalue problems. If there were N actuator arrays, there u 2 u would be N u »=i eigenvalues for the matrix T„ A. d Appendix F Derivation of the Standard QP Format Let the first term of (6.30) be J from (6.28) and (6.17), (for clear notation, k is ignored in the following derivation) u = AV P C ® CP AU-2AVpTC ® £ T r T u + e ® e, T 1 u T 1 l e = CP x -Y . x a (F.l) (F.2) sp Similarly, let the third term of (6.30) be J , from (6.33), 3 J 3 = (S AU + V U(k-l)) Q (S AV + T du p 3 du = AU §lQ S „AU + 2AU§J Q ce+ r r 3 u d C = V U{k-l). 3 V U(k-l)) p f®^ (F.3) (FA) p Note that both the third terms in (F.l) and (F.3) do not relate with AU, thus they can be removed from the objective function (6.30). From (F.l), (F.3) and the second term of (6.30), the standard QP format (6.41) can be obtained. 143 Appendix G Derivation of Equation (6.61) Substituting x oij{k) a Xaiji^k) in (6.24) and y~i(k) (6.25) into (6.23), A ijX ij(k a a ' ^ ^ C jj-A. ijX jj(k 1) a a 1) ~\~ a 3= 1 N u B A (k aij Uj - 1) - ^1 Y.CaijB^Aujik - 1) + ^ydk). (G.l) 3=1 Then, from (6.19), x (k) can be formed as follows: a x (k) a r (G.2) e-AU{k), T-Y{k) + [i-z-'py'-Via^- (G.3) ,a), (G.4) P 7i ^ ( A K I , • • • ,A N ) ai u (G.5) — a • [CailA ii, • a 'L(l) L(2) L(iV )l . iV„ ' iV ' " ' ' 7V J ' r a u e Lbi (G.6) M U [z • I - /?n • [ L J I , • • • ,L ] bNy ^(•Bail, • • • , B N ) ai u — Oi • [C ilB n, a a where V is the block diagonal operation defined in (6.7). 144 (G.7) , l • j C iN B {N ) a u a u , (G.8) Appendix H Relationship Between G sum and G ss For simplicity, we consider T = 0 case (if T ^ 0, we can let w(k) = Uj(k - T ) as the new inputs and get similar results.) From (2.25), A B , and C„ in (6.8) and (6.9) can have a very simple form: dij ijy dij dij i:j A-ij CLijI, (H.l) Gij — I • Bij Then, A , B , and C j in (6.13) are aij aij -A ai Gij aij I 0 aij I I , C, 0/ (H.2) Therefore, from (6.70), Hp—Tdij Gsumij = ^ ] k=l k ^ ] Q>ij Gij 7=1 Hp —Tdij k (H.3) k=l 1=1 The ratio r in (6.71) can be computed from (H.3). tj 145 Appendix I Solution of (6.95) We propose an iterative technique for solving the optimization problem (6.95) alteratively for V and W. The individual problems are straightforward to solve. Step 1: Initialize V = I, W = I, and L = Step 2: Solve V = arg ^ min ^ |\I> - VL\\ , diagonal V Step 3: Solve W = arg Step 4: If V and W have not converged, then go to Step 2. t then form V = V • V, and L = V • L. F t t t min | | * - LW\\ , then form W = W • W diagonal W F u and L = L • W . t The optimization problem in Step 2 is straightforward to solve. If we partition \l> and L as the following: = [tf!, * , • • • , * f 2 Nu , L = [L L , • • • , L f u 2 N , (1.1) and let Vj represent the diagonal elements of V, then it is easy to see that \\*-VL\\ = J2\\*J-VjL \\ , F j (1.2) F 3=1 and note that we can minimize the optimization problem in Step 2 by solving N u inde- pendent problems, m i n ll^j _ jLj\\ , v 2 F for j = 1,2, ••• , (1.3) Vi Solving the optimization problem (1.3) is fairly straightforward, Ny J (1.4) (1.5) i=l (1.6) 146 I. Solution of (6.95) 147 where ipij and kj are elements of \& and Lj respectively. 3 We take derivative of J, ^ =- 2 ^ 3 then, setting i j l l j ) + 2 v ^ l l (1.7) i=l i=l = 0 gives . _ Y,il\{i>iilij) •_ 1 AT The optimization problem in Step 3 is solved in a similar manner. / T O N Bibliography [1] Calculation and partitioning of variance using paper machine scanning sensor measurements. Technical information paper, Technical Association of the Pulp and Paper Industry, TIP 1101-01 1996. [2] M. Alamir and G. Bornard. Stability of a truncated infinite constrained receding horizon scheme: The general discrete nonlinear case. Automatica, 31 (9): 1353-1356, 1995. [3] J. U. Backstrom, C. Gheorghe, G. E. Stewart, and R. Vyse. Constrained model predictive control for cross directional multi-array processes. Pulp & Paper Canada, 102(5):T128-T131, May 2001. [4] J. U. Backstrom, B. Henderson, and G. E. Stewart. Identification and multivariable control of supercalenders. In Proc. Control Systems 2002, pages 85-91, Stockholm, Sweden, June 2002. [5] T. A. Badgwell. Robust model predictive control of stable linear systems. International Journal of Control, 68(4):797-818, 1997. [6] B. Bamieh, F. Paganini, and M. A. Dahleh. Distributed control of spatially invariant systems. IEEE Transactions on Automatic Control, 47(7): 1091-1107, 2002. [7] R. Bartlett, L. Biegler, J. Backstrom, and V. Gopal. Quadratic programming algorithms for large-scale model predictive control. Journal of Process Control, 12:775795, 2002. [8] A. Bemporad and M. Morari. Robust Model Predictive Control: A Survey, volume 245, chapter Robustness in Identification and Control, Lecture Notes in Control and Information Sciences, pages 207-226. Springer Verlag, 1999. [9] C. J. Biermann. Handbook of Pulping and Papermaking. Academic Press Inc., 1996. [10] R. D. Braatz and J. G. VanAntwerp. Robust cross-directional control of large scale paper machine. In Proc. IEEE International Conference on Control Applications, pages 155-160, Dearborn, MI/USA, 1996. [11] C. S. Burrus, R. A. Gopinath, and H. Guo. Introduction to Wavelets and Wavelet Transforms: A Primer. Prentice Hall, 1998. [12] E. F. Camacho and C. Bordons. Mode Predictive Control. Springer, 1999. 148 Bibliography 149 [13] J. C. Campbell, J. B. Rawlings, and C. V. Rao. Efficient implementation of model predictive control for sheet and film forming processes. AICHE Journal, 44(8):17131723, 1998. [14] K. Cutshall. Paper Machine Operations, Pulp and Paper Manufacture, volume 7, chapter XVIII: Cross-direction control, pages 472-506. Atlanta and Montreal, 3rd edition, 1991. [15] P. Dave, F. J. Doyle III, and J. F. Pekny. Customization strategies for the solution of linear programming problems arising from large scale model predictive control of a paper machine. Journal of Process Control, 9(5):385-396, October 1999. [16] P. Dave, D. A. Willig, G: K. Kudva,, J. F. Pekny, and F: J. Doyle III. LP methods in MPC of large-scale systems: Application to paper machine CD control. AICHE Journal, 43(4): 1016-1031, April 1997. [17] P. J. Davis. Circulant Matrices (2nd Edition). Chelsea Publishing, 1994. [18] S. R. Duncan. The Cross-Directional Control of Web Forming Process. PhD thesis, University of London, UK, 1989. [19] S. R. Duncan. The design of robust cross-directional control system for paper making. In Proc. American Control Conference, pages 1800-1805, Seattle, WA, USA, June 1995. [20] S. R. Duncan and G. F. Bryant. The spatial bandwidth of cross-directional control systems for web processes. Automatica, 33(2):139-153, February 1997. [21] J. Fan and G. A. Dumont. A novel model reduction method for sheet forming processes using wavelet packets. In Proc. the 40th IEEE Conference on Decision and Control, pages 4820-4825, Orlando, FL, USA, December 2001. [22] J. Fan, G. E. Stewart, and G. A. Dumont. Model predictive cross-directional control using a reduced model. In Proc. Control Systems 2002, pages 65-69, Stockholm, Sweden, June 2002. [23] J. Fan, G. E. Stewart, and G. A. Dumont. Structured uncertainty in paper machine cross-directional control. Automatica, 2003. (submitted). [24] J. Fan, G. E. Stewart, and G. A. Dumont. Two-dimensional frequency analysis for model predictive control of cross-directional processes, submitted to Automatica, 2003. Bibliography 150 [25] J. Fan, G. E. Stewart, and G. A. Dumont. Two-dimensional frequency response analysis and insights for weight selection in cross-directional model predictive control. In Proc. the J^2nd IEEE Conference on Decision and Control, Maui, Hawaii, USA, December 2003. (to appear). [26] J. Fan, G. E. Stewart, G. A. Dumont, J. Backstrom, and P. He. Approximate steadystate performance prediction of large-scale constrained model predictive control systems, submitted to IEEE Transactions on Control Systems Technology, 2003. [27] A. P. Featherstone and R. D. Braatz. Control-oriented modeling of sheet and film processes. AICHE Journal, 43(8): 1989-2001, August 1997. [28] A. P. Featherstone and R. D. Braatz. Integrated robust identification and control of large-scale processes. Ind. Eng. Chem. Res., 37:97-106, 1998. [29] A. P. Featherstone, J. G. VanAntwerp, and R. D. Braatz. Identification and Control of Sheet and Film Processes. Springer, 2000. [30] C. E. Garcia, D. M. Prett, and M. Morari. Model predictive control: theory and practice - a survey. Automatica, 25(3):335-348, 1989. [31] G. Gavelin. Paper Machine Design and Operation - Descriptions and Explanations. Augus Wilde Publications Inc., 1998. [32] H. Genceli and M. Nikolaou. Design of robust constrained model-predictive controllers with volterra series. AICHE Journal, 41(9):2098-2107, 1995. [33] S. Gendron and J. Hamel. On paper streaks caused by mapping misalignment of the cross-direction control system. In Control Systems 2002, pages 117-122, Stockholm, Sweden, 2002. [34] G. C. Goodwin. When, why and how to constrain control (with application to cross directional control). In Proc. Control Systems 2002, pages 15-23, Stockholm, Sweden, 2002. [35] D. M. Gorinevsky and C. Gheorghe. Identification tool for cross-directional processes. IEEE Transactions on Control Systems Technology, ll(5):629-640, September 2003. [36] D. M. Gorinevsky and E. M. Heaven. Automated identification of actuator mapping in cross-directional control of paper machine. In Proc. the American Control Conference, pages 3400-3404, Albuquerque, New Mexico, USA, 1997. Bibliography 151 [37] D. M. Gorinevsky and E. M. Heaven. Performance-optimized applied identification of separable distributed-parameter processes. IEEE Transactions on Automatic Control, 46(10):1584-1589, October 2001. [38] D. M. Gorinevsky, E. M. Heaven, and R. N. Vyse. Performance analysis of crossdirectional process control using multivariable and spectral models. IEEE Transactions on Control Systems Technology, 8(4):589-600, July 2000. [39] D. M'. Gorinevsky and G. Stein. Structured uncertainty analysis of robust stability for spatial distributed systems. In Proc. the 39th IEEE Conference on Decision and Control, pages 3757 -3762, Sydney, Austrilia, 2000. [40] D. M. Gorinevsky and G. Stein. Structured uncertainty analysis of spatially distributed paper machine process control. In Proc. the American Control Conference, pages 2225-2230, Arlington, VA, USA, 2001. [41] R. M. Gray. Toeplitz and circulant matrices: ee.stanford.edu/~gray/toeplitz.pdf, August 2002. A review. http://www- [42] R. Hall. Coordinated cross-machines control. In TAPPI Proceedings of Process Control Conference, pages 11-19, 1989. [43] A. Halouskova, M . Karny, and I. Nagy. Adatpive cross-direction control of paper basis weight. Automatica, 29:425-429, 1993. [44] 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, January 2002. [45] W. P. Heath. Orthogonal functions for cross-directional control of web forming processes. Automatica, 32(2):183-198, 1996. [46] M. Hovd, R. D. Braatz, and S. Skogestad. SVD controllers for H - , Hoc- and uoptimal control. Automatica, 33(3):433-439, 1997. 2 [47] M. Hovd and S. Skogestad. Control of symmetrically interconnected plants. Automatica, 30(6):957-973, 1994. [48] M. V. Kothare, V. Balakrishnan, and M. Morari. Robust constrained model predictive control using linear matrix inequalities. Automatica, 32(10):1361—1379, October 1996. [49] B. Kouvaritakis, J. A. Rossiter, and J. Schuurmans. Efficient robust predictive control. IEEE Transactions on Automatic Control, 45(8):1545-1549, August 2000. Bibliography 152 [50] K. Kristinsson and G. A. Dumont. Cross-directional control on paper machines using Gram polynomials. Automatica, 32(4):533-548, 1996. [51] 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. [52] D. L. Laughlin. Control System Design for Robust Performance Despite Model Parameter Uncertainties: Application to Cross-Directional Response Control in Paper Manufacturing. PhD thesis, California Institute of Technology, 1988. [53] J. H. Lee. Recent advances in model predictive control and other related areas. In Proc. Chemical Process Control — CPC V, Tahoe City, CA, USA, 1996. [54] J. M. Maciejowski. Predictive Control with Constraints. Pearson Education, 2001. [55] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert. Constrained model predictive control: Stability and optimality. Automatica, 36(6):789-814, June 2000. [56] S. Mijanovic, G. E. Stewart, G. A. Dumont, and M. S. Davis. Robustification of a paper machine cross-directional control system. In Proc. American Control Conference, pages 2203-2209, Arlington, VA, USA, June 2001. [57] M. Misiti, Y. Misiti, G. Oppenheim, and J.-M. Poggi. Wavelet Toolbox. The MathWorks, Inc., 2.2 edition, 2002. [58] M. Morari and N. L. Ricker. Model Predictive Control Toolbox. The MathWorks, Inc., 1.0.7 edition, 1998. [59] K. R. Muske. Steady-state target optimization in linear model predictive control. In Proc. of the American Control Conference, Albuquerque, NM, USA, 1997. [60] K. R. Muske and J. B. Rawlings. Model predictive control with linear models. AIChE Journal, 39(2):262-287, 1993. [61] Z. Nesic. Paper machine data analysis using wavelets. Master's thesis, Dept. of Electrical and computer Engineering, The University of British Columbia, Vancouver, BC, Canada, 1996. [62] A. V. Oppenheim, R. W. Schafer, and J. R. Buck. Discrete-Time Signal Processing (2nd Edition). Prentice Hall, 1999. [63] S. J. Qin and T. A. Badgwell. An overview of industrial model predictive control technology. In Proc. Chemical Process Control—CPC V, Fifth International Conference on Chemical Process Control, Tahoe City, CA, USA, 1996. Bibliography 153 [64] C. V. Rao and J. B. Rawlings. Steady states and constraints in model predictive control. AIChE Journal, 45(6):1266-1278, 1999. [65] J. B. Rawlings. Tutorial overview of model predictive control. IEEE Control Systems Magazine, 20(3):38-52, 6 2000. [66] J. B. Rawlings and K. R. Muske. The stability of constrained receding horizon control. IEEE Transaction on Automatic Control, 38:1512-1516, 1993. [67] A. Rigopoulos. Application of Principal Component Analysis in the Identification and Control of Sheet-Forming Processes. PhD thesis, Georgie Institute of Technology, USA, 1999. [68] O. J. Rojas, G. C. Goodwin, and G. V. Johnston. Spatial frequency antiwindup strategy for cross-directional control problems. Control Theory and Applications, IEE Proceedings, 149(5):414-422, 2002. [69] E. L. Russell and R. D. Braatz. The average-case identifiability and controllability of large scale systems. Journal of Process Control, 12:823-829, 2002. [70] H. Sarimveis, H. Genceli, and M. Nikolaou. Design of robust nonsquare constrained model-predictive control. AICHE Journal, 42(9):2582-2593, 1996. [71] J. Shakespeare. Identification and Control of Cross-Machine Profiles in Paper Machines: A Functional Penalty Approach. PhD thesis, Tampere University of Technology, Tampere, Finland, November 2001. [72] J. Shakespeare, J. Pajunen, V. Nieminen, and T. Metsala. Robust optimal control of profiles using multiple CD actuator systems. In Proc. Control Systems 2000, pages 305-310, Victoria, BC, Canada, May 2000. [73] J. C. Skelton, P. E. Wellstead, and S. R. Duncan. Distortion of web profiles by scanned measurements. In Proc. Control Systems 2000, pages 311-314, Victoria, BC, Canada, May 2000. [74] 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. [75] S. Skogestad and I. Postlethwaite. Multivariable Feedback Control: Analysis and Design. John Wiley & Sons, Inc., 1996. [76] R. Soeterboek. Predictive Control: A Unified Approach. Prentice Hall, 1992. Bibliography 154 [77] G. E. Stewart. Two Dimensional Loop Shaping Controller Design For Paper Machine Cross-Directional Processes. PhD thesis, The University of British Columbia, Vancouver, Canada, 2000. [78] G. E. Stewart, J. U. Backstrom, C. Gheorghe, and R. N. Vyse. Controllability in cross-directional processes: Practical rules for analysis and design. In 87th Annual Meeting, PAPTAC, pages C35-C42, Montreal, PQ, Canada, Feburary 2001. [79] G. E. Stewart, D. M. Gorinevsky, and G. A. Dumont. Robust GMV cross directional control of paper machines. In Proc. American Control conference, pages 3002-3007, Philadelphia, PA, USA, 1998. [80] G. E. Stewart, D. M. Gorinevsky, and G. A. Dumont. H loop shaping controller design for spatially distributed systems. In Proc. IEEE Conference on Decision and Control, pages 203-208, Phoenix, AZ, USA, December 1999. 2 [81] G. E. Stewart, D. M. Gorinevsky, and G. A. Dumont. Spatial loop shaping: A case study on cross-directional profile control. In Proc. American Control Conference, pages 3098-3103, San Diego, CA, USA, June 1999. [82] 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, ll(5):612-628, September 2003. [83] G. E. Stewart, D. M. Gorinevsky, and G. A. Dumont. Two-dimensional loop shaping. Automatica, 39(5):779-792, May 2003. [84] G. E. Stewart, D. M. Gorinevsky, G. A. Dumont, C. Gheorghe, and J. U. Backstrom. The role of model uncertainty in cross-directional control systems. In Proc. Control Systems 2000, pages 337-346, Victoria, BC, Canada, May 2000. [85] G. Strang and T. Nguyen. Wavelets and Filter Banks. Wellesley Cambridge Press, 1996. [86] P. P. Vaidyanathan. Multirate Digital Signal Processing. Prentice-Hall, 1992. [87] J. G. VanAntwerp and R. D. Braatz. Fast model predictive control of sheet and film processes. IEEE Transactions on Control Systems Technology, 8(3):408-417, May 2000. [88] J. G. VanAntwerp, A. P. Featherstone, and R. D. Braatz. Robust cross-directional control of large scale sheet and film processes. Journal of Process Control, 11:149-177, 2001. Bibliography 155 [89] P. Vuthandam, H. Genceli, and M. Nikolaou. Performance bounds for robust quadratic dynamic control with end condition. AICHE Journal, 41:2083-2097, 1995. [90] R. N. Vyse, J. King, and M. Heaven. Consistency profiling — a new technique for CD basis weight control. In Proc. 81th Annual Meeting, technical section CPPA, pages A267-A273, 1995. [91] X. G. Wang, G. A. Dumont, and M. S. Davis. Estimation in paper machine control. IEEE Control Magazine, 13(8):34-43, 1993. [92] A. G. Wills and W. P. Heath. Analysis of steady-state performance for crossdirectional control. IEE Proceedings: Control Theory Application, 149(5):433-440, 2002. [93] K. Zhou, J. C. Doyle, and K. Glover. Robust and optimal control. Prentice Hall, Inc., 1996.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Model predictive control for multiple cross-directional...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Model predictive control for multiple cross-directional processes : analysis, tuning, and implementation Fan, Junqiang 2003
pdf
Page Metadata
Item Metadata
Title | Model predictive control for multiple cross-directional processes : analysis, tuning, and implementation |
Creator |
Fan, Junqiang |
Date Issued | 2003 |
Description | In this thesis, practical techniques for analyzing, tuning and implementing an industrial model predictive controller (MPC) for paper machine cross-directional (CD) processes are developed. These techniques include model reduction, performance and robust stability analysis of a linear closed-loop system, and optimal prediction of steady-state performance for constrained MPC. Paper machine cross-directional processes are large-scale two-dimensional systems. In order to make the online computational load feasible in real time, it is necessary to reduce the dimensions of the model, especially for multiple-array systems. Due to the very large dimension and ill-conditioning of the process, the plant model can be effectively reduced by wavelet matrices based on a modified wavelet packet method. The reduced model captures the controllable spatial components of CD processes. After model reduction, MPC based on the reduced model is implemented in the wavelet domain. Meanwhile, the constraints can be exactly represented in the wavelet domain. The main benefit of implementing MPC in the reduced domain is that the on-line computational time for solving the quadratic programming (QP) problem is greatly reduced compared to the implementation of MPC based on the original model while good control performance is preserved. This method is applicable for multiple-array systems. For large-scale spatially-distributed systems such as CD processes, the process model, the additive structured uncertainty, and the linear portion of the MPC controller are approximated as linear, spatially-invariant, and time-invariant. The transfer function analysis will hold as long as the disturbances are small enough magnitude that the MPC does not hit constraints. To analyze the relevant closed-loop transfer matrices, the novel concept of rectangular circulant matrices (RCMs) is proposed. RCMs can be diagonalized by complex Fourier matrices, allowing analysis in terms of a family of single-input single-output (SISO) transfer functions across the spatial frequencies. Familiar concepts from control engineering such as bandwidth and stability margin are extended into the two-dimensional frequency domain, providing intuitive measures of closed-loop performance and robustness. Consistent criteria are given for the analysis of the closed-loop effect of the industrial CD MPC tuning weights based on standard robust control theory. Properly tuning the magnitude of weights, choosing the structure of weights, and considering the structure of model uncertainty in the MPC design stage can greatly improve the performance. This method can be used to design robust CD MPC for multiple-array systems. In order to assess the steady-state performance of constrained CD MPC, the state of the art requires to run closed-loop simulations. Due to the large-scale characteristic of CD processes, especially for multiple-array systems, it is very time-consuming and inconvenient for use as a practical tuning tool. However, fast and correct prediction of the steady-state performance is necessary for tuning the CD MPC, especially for the cases with active constraints. A one-step static optimizer is proposed to predict the steady-state closed-loop performance of CD MPC. Two examples are given for demonstrating that the static optimizer is significantly more efficient (up to two orders of magnitude) than the conventional closed-loop simulation method while reliably and accurately predicting the steady-state closed-loop performance. |
Extent | 13527762 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0091448 |
URI | http://hdl.handle.net/2429/15135 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2003-860167.pdf [ 12.9MB ]
- Metadata
- JSON: 831-1.0091448.json
- JSON-LD: 831-1.0091448-ld.json
- RDF/XML (Pretty): 831-1.0091448-rdf.xml
- RDF/JSON: 831-1.0091448-rdf.json
- Turtle: 831-1.0091448-turtle.txt
- N-Triples: 831-1.0091448-rdf-ntriples.txt
- Original Record: 831-1.0091448-source.json
- Full Text
- 831-1.0091448-fulltext.txt
- Citation
- 831-1.0091448.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0091448/manifest