UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Multivariable predictive control of a TMP plant Du, Huaijing 1998

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

Item Metadata

Download

Media
831-ubc_1998-345203.pdf [ 9.63MB ]
Metadata
JSON: 831-1.0065334.json
JSON-LD: 831-1.0065334-ld.json
RDF/XML (Pretty): 831-1.0065334-rdf.xml
RDF/JSON: 831-1.0065334-rdf.json
Turtle: 831-1.0065334-turtle.txt
N-Triples: 831-1.0065334-rdf-ntriples.txt
Original Record: 831-1.0065334-source.json
Full Text
831-1.0065334-fulltext.txt
Citation
831-1.0065334.ris

Full Text

M U L T I V A R I A B L E P R E D I C T I V E C O N T R O L OF A T M P P L A N T By Huaijing Du B . A . Sc. Beijing University of Chemical Technology (1983) M . A . Sc. Beijing University of Chemical Technology (1988) A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E STUDIES D E P A R T M E N T O F E L E C T R I C A L A N D C O M P U T E R E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A October 1998 © Huaijing Du, 1998 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 gieefric*/ U^-J (Z&^rSfc,^ /_-~^,r/^y The University of British Columbia Vancouver, Canada Date / 9 9 % . / O. _2._? DE-6 (2/88) Abstract This thesis describes the development of a novel control strategy for a two-stage thermo-mechanical pulping (TMP) plant. Desired pulp quality is achieved by selecting the set-points of specific energy and refining intensity at each stage. The targets of specific energy and refining intensity are obtained through the control of motor load, production rate and refining consistency by manipulating closing pressure, chip flow rate and dilution flow rate at the inlet of each stage. A constrained predictive controller is developed based on the generalized predictive control (GPC) algorithm because of its simplicity, ease of use and ability to handle problems in one algorithm. Future control actions are determined by minimizing the predicted errors without violating input and output constraints. A multi-input multi-output (MIMO) C A R I M A model identified through identification ex-periments is used to predict the future process outputs. Model parameters are estimated on-line to handle a time-varying nature of the process. A n analytical solution of a constrained MIMO G P C subject to input and output con-straints is derived by solving a quadratic programming (QP) problem. The computation required by the analytical solution is substantially lower than that required by an algo-rithmic solution. For general cases of constrained MIMO G P C , an optimal solution is derived by solving a mixed-weights least-squares (MWLS) problem. In the use of M W L S , a control performance index can be easily augmented and the weights can be modified in a manner that encompasses both the requirements for the future control movements to lie inside the feasible region and to minimize the control performance index. If the constrained optimization problem is unfeasible, the M W L S will converge to the point that minimizes the maximum constraint violation. The proposed control schemes were i i tested on the simulation model developed using the mechanistic and empirical methods to describe the behavior of a real process. Simulations demonstrated the proposed control schemes' efficiency and capability of handling problems in one algorithm. A linear model-based control strategy may lead to system oscillation or even insta-bility if a refiner in the process is operated near maximum load point because at this point the nonlinearity between refiner motor load and plate gap becomes severe. To overcome the problem, a nonlinear Laguerre model - a type of orthonormal functions - is used to represent the nonlinear relationship. A MIMO Laguerre model-based predictive controller is then derived as an alternative for the control of a wood chip refiner. The Laguerre model can be arranged in linear form in model parameters so that the standard recursive least squares (RLS) can be used for on-line parameter estimation according to which controller parameters are adjusted. Simulations demonstrated that in the use of the Laguerre model-based control scheme, the nonlinearity of the process can be repre-sented appropriately and the plate clashes resulting from the severe nonlinearity can be avoided. In addition, in the use of the Laguerre function representation the dynamics of an actual process can be described appropriately without the need for assumptions about the plant order and the time delay, i.e., accurate assumptions about their true values are not necessary. i i i Table of Contents Abstract ii List of Tables ix List of Figures x List of Symbols xiv Acknowledgment xvi 1 Introduction 1 1.1 General Introduction 1 1.2 Background and Literature Review 2 1.2.1 Mechanical Pulping Processes 2 1.2.2 Refining Theory 6 1.2.3 Control of A Refiner Mechanical Pulping Process 8 1.2.4 Control Problems 14 1.3 Objective of this Study 16 1.4 Contributions of this Thesis 17 1.5 Outline of the Thesis 20 2 Process Variables, Modelling and Simulations 22 2.1 Introduction 22 2.2 Process Variables 23 iv 2.2.1 Key Manipulated Variables 23 2.2.2 Key Operating Variables 25 2.2.3 Key Pulp Quality Variables 27 2.2.4 Disturbance Variables 28 2.2.5 Variable Interactions 29 2.3 Modelling a T M P Process 30 2.3.1 Modelling Throughput 30 2.3.2 Modelling Consistency 31 2.3.3 Modelling Energy Input 35 2.3.4 Modelling Pulp Properties 40 2.4 Simulation Results and Discussions 40 2.4.1 Simulating the Primary Refiner 41 2.4.2 Simulating the Secondary Refiner 49 2.5 Summary and Conclusions 54 3 Industrial TMP Process Identification 55 3.1 Introduction 55 3.2 Process Model 57 3.3 Analysis Tools 58 3.4 Identification Procedures 58 3.4.1 Model Structure Identification 59 3.4.2 Parameter Estimation 61 3.4.3 Model Checking 61 3.5 Industrial T M P Process Identification 62 3.5.1 Process Description and Experiment Design 62 3.5.2 Identification Results 63 v 3.5.3 Result Discussions 76 3.6 Summary and Conclusions 77 4 Constrained Multivariable Model-Based Predictive Control 78 4.1 Introduction 78 4.2 Unconstrained Generalized Predictive Control 81 4.2.1 Basic SISO Generalized Predictive Control 81 4.2.2 Basic M I M O Generalized Predictive Control 87 4.3 Recursive Parameter Estimation in Adaptive Control 90 4.4 Control under Constraints 91 4.5 Constrained Control via Quadratic Programming 95 4.5.1 Quadratic Programming via Kuhn Tucker Multipliers Computation 96 4.5.2 Analytical Solution for Constrained SISO G P C 98 4.5.3 Analytical Solution for Constrained M I M O G P C 99 4.6 Constrained Control via Mixed-Weights Least-Squares Algorithm . . . . 102 4.6.1 Preliminary 104 4.6.2 Constrained Control via Mixed-Weights Least-Squares 107 4.6.3 Control Law 110 4.7 Simulations 112 4.8 Summary and Conclusions 122 5 Laguerre Function-Based Adaptive-Predictive Control 123 5.1 Introduction 123 5.2 The Laguerre Functions 126 5.2.1 The Linear Laguerre Functions 126 5.2.2 The Nonlinear Laguerre Functions 129 5.3 SISO Laguerre Function-Based Predictive Control 131 vi 5.3.1 Linear Model-Based Control 131 5.3.2 Nonlinear Model-Based Control 135 5.4 M I M O Laguerre Function-Based Predictive Control 138 5.4.1 Mixed Linear-Nonlinear Laguerre Model 138 5.4.2 M I M O Laguerre Function-Based Predictive Control Law 140 5.5 Summary and Conclusions 143 6 Constrained Multivariate Control Of A TMP Plant 144 6.1 Introduction 144 6.2 Choice of Variables and Overall Control Strategy 146 6.2.1 Choice of Variables 146 6.2.2 Constraints 148 6.2.3 Control Strategy Over Two Stages 149 6.3 Control of T M P Refiners 151 6.3.1 Process Model 151 6.3.2 Controller Development 156 6.3.3 Simulation Results and Discussions 166 6.4 Laguerre Function-Based Control of T M P Refiners 177 6.4.1 Nonlinearity Description 177 6.4.2 Past Control Solutions 178 6.4.3 2 x 2 Mixed Linear-Nonlinear Laguerre Function-Based Control Of T M P Refiners 180 6.4.4 Simulation Results and Discussion 181 6.5 Summary and Conclusions . 184 7 Conclusions and Future Work 190 7.1 Conclusions 190 vii 7.2 Future Work . 192 Bibliography 194 Appendices 204 A Time Series Analysis Tools 204 A . l Auto-Correlation Functions 204 A.2 Partial Auto-Correlation Functions 205 A . 3 Cross-Correlation Functions 206 B Optimal Solution via Analytical Quadratic Programming 208 B. l Constraints 208 B.2 Optimal Solution via Analytical QP 210 C Constraint Mapping 218 viii List of Tables 2.1 Primary Refining Conditions 42 2.2 First Stage Refining at Fd = A.7kg/sec 43 2.3 First Stage Refining at E = ZAGJ/t 45 2.4 Effects of First Stage Refining on Pulp Properties, E=3.4 G J / t 48 2.5 Secondary Refining Conditions 49 2.6 Effects of Secondary Stage Refining at E=2.5 GJ / t 53 2.7 Effects of Secondary Stage Refining on Pulp Properties, E=2.5 G J / t . . . 53 3.8 M L / C P Parameter Estimates, Estimated Standard Deviations and F P E . 64 3.9 M L / T S Parameter Estimates, Estimated Standard Deviations and F P E . 68 3.10 M L / D W Parameter Estimates, Estimated Standard Deviations and F P E 72 3.11 Identification Results 76 4.12 Kuhn-Tucker Multipliers for Different Cases 98 ix List of Figures 1.1 A simplified version of two-stage wood chip refining process 5 1.2 A double-disc refiner 6 1.3 Plate surface of a typical primary refiner 7 2.4 Interactions between manipulated variables, operating variables, pulp qual-ity variables and disturbances for a two-stage refining 29 2.5 Mass and energy flows in the lst-stage refiner 32 2.6 Model structure of the lst-stage wood-chip refining 41 2.7 The 1st stage refining at F d=4.7 kg/sec: effects of closing pressure Pc and chip screw speed Sc ( - - :SC = 25 rpm (Fc=271 t/d), -.-.: Sc = 27.5 rpm (F c=298 t/d) and —:5 C = 30 (Fc=325 t/d)) 44 2.8 The 1st stage refining at E = 3.4 GJ / t : effects of dilution flow rate Fd and chip screw speed Sc (- -: 5C=25 rpm (Fc=271 t/d), 5C=27.5 rpm (F c=298 t/d) and — : Sc=30 (Fc=325 t/d)) 45 2.9 The 1st stage refining at E=3.4 GJ / t : effects of rotational speed w and inlet consistency c{ ( — : a = 20%, d = 23.5% and - -: a - 27%) . . 47 2.10 The 2nd stage refining at E2 = 2.5 GJ / t : effects of the rotational speed w and the inlet consistency Cj ( — :c ; = 26%, -.-.: Q = 30.5% and - -: c» = 35%) 51 2.11 Two-stage refining at E2 = 2.5 GJ / t : effects of the 2nd-stage rotational speed w and the inlet consistency Cj ( — : c, = 26%, -.-.:c; = 30.5%, and --: C i = 35%) 52 3.12 Input to and output from a dynamic process 56 3.13 Motor Load (ML) response to Closing Pressure (CP) 65 3.14 M L / C P impulse response estimate 66 3.15 M L / C P residual correlations (d=0) 66 3.16 M L / C P residual correlations (d=l) 67 3.17 C P / M L input, output and model prediction 67 3.18 Motor Load (ML) response to Transfer Screw Speed (TS) 69 3.19 M L / T S impulse response estimate 70 3.20 M L / T S residual correlations 70 3.21 M L / T S residual correlations (d=l) 71 3.22 M L / T S input, output and model prediction 71 3.23 Changes in Valve Positions (VP), Dilution Water (DW) and Motor Load (ML) 73 3.24 Motor Load (ML) response to total Dilution Water (DW) 74 3.25 M L / D W impulse response estimates 74 3.26 M L / D W residual auto- and cross-correlations 75 3.27 M L / D W residual auto- and cross-correlations (d=l) 75 3.28 Dynamic and static relationships between inputs and outputs 77 4.29 Basic structure of constrained model-based predictive control 95 4.30 Constrained G P C via analytical QP (t > 300 sampling intervals: -1.5 < u t < 1.5, -2.0< Aut <2.0) 113 4.31 Constrained G P C via M W L S ( t > 300 sampling intervals: -1.5< ut < 1.5, -2.0 < Aut < 2.0) 114 4.32 Criterion for stopping the iteration 116 4.33 Constrained G P C via analytical QP (t > 300 sampling intervals : -4.0 < «i < 4.0, -4.0< u2 < 4.0, | A u i | < 4.0, | A u2\ < 4.0) 118 xi 4.34 Constrained G P C via M W L S (t > 300 sampling intervals: -4.0 < ux < 4.0, -4.0< u 2 < 4.0, | A M I | < 4.0, | A u2\ < 4.0) 119 4.35 Output constrained G P C via analytical QP (t > 300 sampling intervals: -1-0 < yi < 1.0) 120 4.36 Output constrained G P C via M W L S (t > 300 sampling intervals: -1.0< y i < 1.0, -1.0 <y2< 1.0) 121 5.37 Laguerre ladder networks 127 5.38 The discrete-time nonlinear Laguerre network 131 5.39 Linear model-based control will lead to oscillation when the nonlinearity is represented by the dashed curve [1] 135 6.40 Specific energy is the area under curve of N plotted against / 145 6.41 Block diagram for the control of a two-stage T M P plant 149 6.42 Constrained multivariable adaptive-predictive control of the primary refinerl58 6.43 Input-output dynamic and static relationship for the primary refiner . . . 166 6.44 Constrained control via M W L S , t> 200 sec: | A ux\ <50, | A u2\ <0.2, | A u31 <0.5 ( ?i = q2 = 1, p2 = 1.0, pi = p 3 = 0.0001) 168 6.45 Constrained control via analytical QP, t>200 sec: 820 < ux < 1150, 3.4 < u 3 < 5.2 (qi = q2 = l , p 2 = 1.0, P l = p 3 = 0.001) 169 6.46 Unconstrained controller (qx = q2 = l , p 2 = 1-0, pi = P3 = 0.0001) . . . . 171 6.47 Unconstrained controller (qx = q2 = l , p 2 = 10, pi = p 3 = 0.0001) 172 6.48 Unconstrained controller (p2 = l , p i = p 3 = 0.0001) 175 6.49 Unconstrained controller (p2 = l , p i = p 3 = 0.0001) 176 6.50 The gain between motor load and plate gap is subject to slow drift and a sudden change in the sign from negative to positive 177 6.51 Structure of a 2 x 2 Laguerre model for a T M P refiner 180 xii 6.52 Motor load setpoint unreachable 185 6.53 A plate clash was avoided by the proposed control scheme 186 6.54 Response to the changes in consistency setpoint 187 xiii List of Symbols The variables appearing in Section 2.3 are given in the follows: Q : the volumetric flow rate of inlet wood chips (dry wood + moisture) ( v : the chip velocity in chip transfer screw direction(m/s) Sc : the chip screw rotational speed (revolutions/min or rpm) pc : the chip bulk density of incoming wood chips(A;p/m3) fc : the wood chip flow rate at inlet of a refiner(fcg/s) s : the solid content of incoming chips(%) ec : the energy in inlet chips(/sJ/s) ed : the energy in dilution water(A;J/s) esw : the energy in the seal water(/cj/s) P : the net power drawn by a motor(k J/s) 77 : the power efficiency coefficient(%) e3i : the energy in back flow steam(k J/s) ea2 : the energy in forward steam(A;J/s) ep : the energy in outflow pulp(fc J/s) Fc : the dry chip flow rate to a refiner(kg/s) ff : the dry fiber flow rate in outlet \oulp(kg/s) fwi : the water flow rate in inlet chips(fc^/s) fw2 • the water flow rate in outlet pulp (kg/s) Fd : the dilution water flow rate(kg/s) Fsw : the seal water flow rate(fc#/s) xiv Tc : the inlet chip temperature (°C) Tp : the outlet pulp temperature (°C) Td : the dilution water temperature(°C) Tsw : the seal water temperature(°C) c c : the specific heat of dry chips(A;J/A;5'.0C) Cf : the specific heat of dry fiber(A:J/A;(7.0C) cw : the average specific heat of water(20° — 150°C)(kJ/kg.°C) L : the enthalpy of steam(kJ/kg) E : the specific eneTgy(kJ/kg) T : the residence time(s) w : the refiner disc rotational speed(radian/s) Ci : the inlet consistency(%) v\ : the inlet radii of refining zone(m) r2 : the outlet radii of refining zone(m) Ls : the latent heat of steaxa(kJ/kg) N : the number of refiner bar impacts received by unit pulp J : the average specific energy per refiner bar impact or refining intensity(A; 3j kg) e : the specific refining power(k J/kg/s) xv Acknowledgment M y greatest thanks go to my senior supervisor, Prof. Guy A . Dumont for his invaluable guidance and advice in my progress throughout the research, for his steady assistance, supervision and support in carrying out this study, for his continuous encouragement and help in my Ph.D. study, and for his valuable time and attention he devoted to me. I would also like to thank Dr. Patrick Tessier, Dr. Bruce Allison and Dr. Ye Fu for their stimulating discussion and helpful suggestions in my thesis work. I wish to express my gratitude to the members of my supervisory committee, Prof. Guy A . Dumont, Prof. Michael Davies and Prof. Ezra Kwok for their valuable discussions and assistance which proved to be invaluable for me to finish my thesis. I thank all the members of Control Group from the U B C Pulp and Paper Centre for their valuable inputs in my thesis work during my stay at U B C . Finally the financial support from Dr. Guy A . Dumont and NSERC Wood-Pulps Network Centres of Excellence is gratefully acknowledged. xvi To my husband, parents and brothers Huaijing, with love Chapter 1 Introduction 1.1 General Introduction In the pulp and paper industry, wood pulping refers to the process by which wood chips are reduced to individual fibres. This task can be accomplished chemically (chemical pulping) or mechanically (mechanical pulping) or by combining these two methods. Each of them produces substantially different fibre characteristics. In chemical pulping, lignin which binds the wood fibres together is dissolved chemi-cally. The resultant fibres are longer, more flexible and considerably stronger than the mechanical fibres. During the manufacturing of chemical pulp, most of the lignin is suc-cessfully removed. Components other than lignins, such as hemicellulose and cellulose, are removed simultaneously, hence the yield of chemical pulping is lower (40-45 %) than the yield of mechanical pulping (up to 95 %) [2]. In mechanical pulping, the lignin is fractured out and not removed. Therefore, more wood constituents, such as fibre bundles, damaged fibres and fibre fragments with some whole fibres, are retained. This results in a higher yield (up to 95 %) in converting the wood into pulp [2]. Because of the mixed nature of the particles in the pulp, mechanical pulp has high light scattering coefficients as well as high bulk and relative stiffness [3]. However, mechanical pulp has small average particle length. The resulting paper conse-quently has low strength and high bulk, and discolors easily. In order to produce pulp 1 Chapter 1. Introduction 2 with the physical properties desired for higher quality printing, mechanical pulping com-bined with some chemical treatment can provide an alternative to mechanical pulping. The major drawback of mechanical pulping is that a greater amount of power is required than in chemical pulping. Hence equipment utilization, pulp quality level and product uniformity must be increased. This can only be accomplished by modifying the process, and by improved equipment and process control [3]. This thesis partially addresses the problems. In particular, the thesis attempts to develop a suitable control strategy for a two-stage thermo-mechanical pulping (TMP) plant. The purpose of this chapter is to introduce the problems and state the thesis ob-jectives. Section 1.2 includes the background material relevant to mechanical pulping processes and their control. Section 1.2.1 gives a brief introduction to mechanical pulp-ing processes including unit operations; Section 1.2.2 describes refining theories; Section 1.2.3 briefly reviews the past and existing control schemes; and Section 1.2.4 states some existing control problems. Thesis objectives, contributions and outline are given in Sec-tions 1.3, 1.4 and 1.5. 1.2 Background and Literature Review 1.2.1 Mechanical Pulping Processes Depending on the machinery and action utilized, mechanical pulps can be produced by two different processes: grinding and refining. Grinding or stone groundwood is the oldest method, in which wood logs are forced against a rapidly revolving roughened grindstone and converted to individual fibres. In refining, wood chips are fed between two metal discs (at least one of which rotates) of a refiner and converted to individual fibres. These two processes result in significantly different pulp characteristics. Groundwood pulp has a higher content of fine material due to the abrasive action, whereas refiner Chapter 1. Introduction 3 pulp has a smaller content of fine material [4] but a higher content of long fibres. As a result, refiner pulping produces much stronger fibres than stone groundwood. Various types of refiner mechanical pulp can be obtained by modifying a refiner pulping process. Thermo-mechanical pulping (TMP) is a modification of a refiner mechanical pulping (RMP) process. It involves steaming the raw material for a short period of time prior to refining to soften the chips. The resulting pulp has a greater percentage of long fibres and less shives than those of an R M P . Chemi-thermo-mechanical pulping (CTMP) is a modification of a T M P with some chemical treatment in order to obtain pulp with desired physical and chemical properties for higher grade applications. Due to the thermal pretreatment to the chips and chemical action in separation of fibres, C T M P produces pulp with a higher content of longer fibres and a higher brightness. Unit Operation A thermo-mechanical pulping process generally involves three main operation areas: 1. wood chip pretreatment 2. wood chip refining 3. pulp processing Wood chip pretreatment consists of chip screening to remove under or oversize ma-terial; chip washing to remove rocks, metal and sand; chip steaming to soften lignins binding the fibres so that produced pulps have a greater percentage of long fibres and less shives. In chemi-thermo-mechanical pulping, the wood chips are impregnated with chemicals to improve pulp brightness and strength. Wood chip refining aims at breaking chips into individual fibres (see following descrip-tion). Chapter 1. Introduction 4 Pulp processing is aimed at enhancing and controlling pulp quality. It consists of: latency removal to straighten the fibres; pulp screening and reject refining to remove unrefined fibre bundles; pulp cleaning to remove heavy contaminants; pulp washing to remove wood resins and metallic ions; and pulp bleaching to increase brightness. Wood Chip Refining Generally, wood chip refining is carried out in two stages as shown in Figure 1.1. Wood chips (dry wood and water moisture) are introduced into the open eye of the primary refiner by the chip transfer screw. As the wood chips pass through the refining zone between the rotating plates (or discs), they are progressively broken down into smaller particles and finally into individual fibres. In the secondary refiner, the fibres are further refined according to the pulp quality requirements. Simultaneously, dilution water is added to both primary and secondary refiners to control the consistencies in the refining zones. Plate gap for each refiner is controlled by either an electro-mechanical or hydraulic loading system mounted in the refiner. Following each of the refiners, cyclones are in-stalled to separate steam from the produced pulp and to cool the hot pulp to prevent brightness loss. A wood chip refiner, a key element in the pulping process, is either a single-disc (one rotating disc and one stationary disc), a double-disc (two discs rotating in oppo-site directions, Figure 1.2) or a twin (a rotating double-sided disc between two station-ary discs) configuration. A refiner disc, driven by a large motor, usually rotates at 1200 — 1800 r/min. Refiners are currently available with up to 70 inch (1.8m) diameter and 18, OOOif P(13,400kW) of power supplied to each plate [2]. A large-diameter high-powered refiner is generally developed for large-capacity refining. Plate gap is accurately controlled to about 0.1mm and is always in the range of 0.05 to 0.2mm [5]. Refiner plates are generally divided into three sections: the breaker bar section, the intermediate Chapter 1. Introduction 5 chip silo pre-steaming heating cyclone cyclone v/ / chip impres- transfer primary 2nd stage latency washer safiner screw refiner refiner tank Figure 1.1: A simplified version of two-stage wood chip refining process bar section and the fine bar section (Figure 1.3). The breaker bar section (closest to the eye of the refiner) has wide bars with deep grooves, serving to break up the chips as they enter the refiner. The intermediate and fine bar sections consist of progressively narrower bars and shallower grooves, which serve to refine continuously the partly refined pulp. The three sections may vary in size depending on the refining stage. Plate life is generally in the range 500 - 1000 hours and varies with different material and refining stage. The plate life of a primary refiner is longer than that of a secondary refiner because refining action in the secondary refiner is greater than that in the primary refiner. In the refining zone (or the fine bar section), most refining is done at consistencies in the range of 16% to 50% [2]. Most of the energy introduced by a refiner is converted into the heat which produces hot water and steam. The steam produced in the refining zone develops pressure substantially above those at either the inlet or the discharge. This amount of steam can be indicated by the temperature measurement in the refining zone. The peak temperature may range from 115 °C to 145 °C, equivalent to a peak pressure in the order of 170 kPa to 446 kPa [6]. Additional information on refiner mechanical pulping can be found in [2, 4, 7, 8]. Chapter 1. Introduction 6 chips steam water *• pulp primary refiner Figure 1.2: A double-disc refiner 1.2.2 Refining Theory Intensive theoretical studies have been made by many researchers towards developing a unified theory of refining. Much work done in the past was based on the product quality-energy consumption relationship, i.e., pulp quality such as freeness was determined by the energy applied. However, recent work showed that specific energy alone does not fully determine pulp properties and refining action [9, 10, 11, 12, 13]. Many researchers agree that pulp quality is determined by two basic variables: the energy per unit mass of pulp or specific energy (E) and the specific energy per bar impact or refining intensity (/), rather than by the specific energy alone [5, 14, 3, 15, 9, 10, 11, 12, 13]. For a given refiner design and a constant disc rotational speed, the refining intensity is largely dependent on refining consistency while it is insensitive to the specific energy [16]. The influence of the refining intensity on pulp quality in low consistency refining has been known for a long time [10], but it has only recently been pointed out by Lunan et Chapter 1. Introduction 7 Figure 1.3: Plate surface of a typical primary refiner al. [17] that it may also be important in determining the pulp quality in high consistency refining. Lunan's idea has been further investigated by Miles [10], which showed that for a given specific energy under a wide variety of operating conditions, the refining intensity was strongly correlated to the pulp quality. The work by Miles et a/.[10,11] and Rodarmel [18] also showed that for a given specific energy, low intensity refining produces pulps of higher strength but lower opacity, while high intensity refining results in pulps of lower strength but better printing properties. It was suggested by Miles et al. recently that at a given total specific energy, the combination of high intensity-low energy refining in the first stage followed by low intensity-high energy refining in the second stage could produce pulp with improved printing properties with no loss in strength [11]. Despite significant progress by many researchers last few decades, the refining action in the refining zone is still not well understood. To optimize the refining processes, more theoretical and experimental work on the mechanics and dynamics of the fibres between two plates of a refiner are needed. Chapter 1. Introduction 8 1.2.3 Control of A Refiner Mechanical Pulping Process The control of mechanical pulping processes has been in a state of ongoing develop-ment since the mid-1970s [19]. Significant progress has been made towards refining optimization, cost reduction and product quality improvement. The following briefly re-views some control strategies and applications available in the literature. Additional information on the control of refiner mechanical pulping processes can be found in [20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Past Control Strategies Specific Energy Control Specific energy is a main variable affecting pulp properties. A traditional pulp quality control strategy is by manipulating specific energy based on the quality-energy relation-ship. Specific energy is defined as the ratio of net power or motor load applied to dry fibre mass throughput or production rate. Therefore, the target of specific energy is often achieved through the control of production rate or motor load. Many mills use chip screw speed to set the desired production rate. However, the production rate may still vary due to changes in the quality of raw material such as chip bulk density and moisture content. Motor Load Control The motor load is most sensitive to plate gap, and therefore one way of setting and maintaining the motor load is to adjust plate gap. Various control approaches [27, 30, 31, 32, 33] have been proposed to solve motor load control problems. With these control approaches, the motor load was controlled through manipulating the plate gap while the throughput was kept constant and the refining consistency was not under control. Although load controls via gap adjustment will reduce load variations, the motor load can still vary due to the variations in raw material quality such as chip size, Chapter 1. Introduction 9 moisture and density. To deal with the fluctuations in raw material quality, an alternative has been proposed [34, 35, 36, 37], where the motor load was controlled by automatically adjusting chip feed indicated by the chip screw speed. Using this control approach, the lower frequency load fluctuations induced by wood quality variations can be reduced, but not the higher frequency fluctuations in the motor load [34]. Therefore, further work needs to be done to reduce the higher frequency fluctuations in the motor load. Consistency Control Refining zone consistency is known to have a major effect on pulp properties at given specific energy. A n improved consistency control will significantly enhance the refiner performance by eliminating the unavoidable feeding disturbances [38], while a poor consistency control can lead to unstable refiner operation [28]. Various control strategies have been proposed with the help of on-line consistency sensors [39, 38, 40, 41, 35, 36], where the consistency at the blow-line was controlled by manipulating the dilution flow at the inlet. Improved consistency control can provide a significant reduction in the fluctuation of the motor load and obtain large improvements in stabilizing T M P / C T M P plants, thus resulting in energy savings, control efficiency and improved pulp quality [41, 40]. Pulp Quality Control The pulp quality can be defined by many property variables, but in many cases freeness is still seen as an important variable to characterize the strength of mechanical pulp [42]. Various control strategies have been proposed [26, 27, 43, 44, 42, 45] in the area of freeness control. The objective of these works was to control the freeness to the quality target by adjusting the specific energy while minimizing freeness deviation from the target as well. Freeness was controlled by adjusting the specific energy setpoint through manipulating the motor load or the production rate. Simulation studies and industrial applications showed that the proposed controls brought about freeness variation reduction Chapter 1. Introduction 10 and some improvement in pulp quality. However, measuring and controlling pulp quality by freeness alone is not satisfactory because pulp quality is determined by at least two property variables, rather than by freeness alone. A n alternative to the single-variable control has been proposed [19], where both freeness and long fibre content were controlled by selecting the optimal setpoints of a T M P operation through a model-based analysis system. It showed the potential for multivariable quality control of the refining operation, but one of the weaknesses of this approach was maintaining accurate information flow to the model [19]. Advanced Control Applications Proportional-Integral-Derivative (PID) control is the most prevalent form of industrial control. However, standard PID techniques are not suitable for the control of the pro-cesses with time delay, nonlinearity, multivariable and time-varying gain/dynamics. A l -ternatively, advanced control methods will be able to handle the problems more efficiently and help to ensure optimum control performance. Following is a brief review of some advanced controls available in the literature. Additional information in advanced control applications can be found in [21, 46, 23, 47]. Adaptive Control Motor load is most sensitive to plate gap adjustments. However, control of the motor load by adjusting the plate gap is difficult due to the fact that the gain between the load and the gap is subject to a slow drift as well as a sudden change in sign. To deal with the problem, various control schemes have been proposed [24, 27, 30, 48]. Dumont first applied an adaptive control scheme to the problem [30], where a Dahlin regulator was tuned using the gain estimate through a recursive least-squares estimator with a variable forgetting factor. Industrial trials showed that the proposed regulator was able to track the slow gain drift and a sudden change in its sign [30, 43]. However, one problem with Chapter 1. Introduction 11 this technique is the tedious tuning of the forgetting factor algorithm parameters, so it may be unreliable in a continuous usage environment [32]. Dual Adaptive Control To improve the reliability of the adaptive controller proposed in [30], Dumont and Astrom [32] proposed a dual control strategy, where additional terms reflecting the non-linear nature of the process were introduced into a performance index, rather than using an index function reflecting the output error only. In the case of an unreachable motor load setpoint, the controller will tend to keep the gain negative instead of eliminating the output error. The proposed dual control strategy was further modified with an improved criterion by Allison et al. [48]. Preliminary results of an industrial refiner control imple-mentation showed the general success. However, one problem with this technique was the computational difficulty attached to solving the dual control problem on-line [31]. Adaptive Inferential Control A n adaptive inferential control strategy was proposed by Kooi [49] for the closed loop control of freeness through manipulating the gap. The control strategy was for the cases when the measurement of freeness is not readily available either due to lack of a reliable on-line sensor or the long sampling time to obtain the test result. Wi th this strategy, rapidly measured specific energy was used to infer the freeness based on a dynamic-noise model which linked the plate gap and the specific energy to the freeness. Due to the time-varying and stochastic natures of the process, the controller parameters were tuned based on on-line estimates of the model parameters. However, this single-variable control approach was limited by its definition of pulp quality. Minimax Robust LQ Control A linear quadratic (LQ) control based on multimodelling techniques was proposed by Toivonen and Tamminen [45] for the control of freeness in a T M P plant. The process was assumed to be described by a set of linear models presenting different operating Chapter 1. Introduction 12 situations. The control problem was then reduced to a single constant control law so that satisfaction of control performance in terms of the output error was obtained for all set of the models. The idea was achieved by formulating a minimax control problem, in which the maximum output error of a given sets of models was minimized subject to a constraint on the maximum input variance. The proposed control strategy was tested in a T M P plant for the control of freeness through manipulating the chip feed. The results showed that the performance of the proposed controller was better than that of a manual control or the control by a Pi-regulator [45]. However, this approach is only suitable when the range of the parameter changes can be estimated and when the changes are not too large. Neural Networks In general, the objective of a neural network based control approach is to develop an algorithm to adjust the parameters of the network based on a given set of inputs and outputs. The error between the desired output and the computed output at the output layer of the network is back propagated to the hidden and input layer for weight updating. The use of neural networks with back-propagation learning was proposed by Kooi and Khorasani [50] for the control of the specific energy in a T M P plant. The proposed controller consisted of one processing element in each input and output layer, and two hidden layers with 10 processing elements in each hidden layer. It was reported in [50] that this neural network-based controller was robust in the sense that it needed no prior knowledge of process order, time delay, dynamics and noise. However, the use of a neural network model usually requires a large number of weights and the weights are not in a linearized form, thus the estimation converges slowly [51]. The performance index via model parameters also becomes complicated. These characteristics are the main difficulties when applying a neural network model in adaptive control. Laguerre Function-Based Control Chapter 1. Introduction 13 A predictive control strategy based on a nonlinear Laguerre model was introduced by Fu and Dumont [51] to handle the nonlinear relationship between refiner motor load and plate gap. The controller parameters were tuned based on the estimates of the La-guerre model parameters. Simulation results showed that in the case of an unreachable motor load setpoint, the proposed controller was always trying to keep the gain negative instead of eliminating the output error [51]. The above single-variable control idea was later extended into a 2 x 2 system by Du et al. [52] to overcome the limitation with a single-variable control philosophy. With this control approach, the motor load and the outlet consistency were controlled by adjusting the gap through manipulating the closing pressure, and the dilution flow rate. The process model was based on a mixed linear-nonlinear Laguerre function representation, where the nonlinear Laguerre functions were used to model the gap-load nonlinear relationship whereas the linear Laguerre functions were used to model the dilution-consistency relationship. The predictive controller pa-rameters were tuned based on on-line estimates of the model parameters. Simulation results showed that the use of this approach was not only able to handle the nonlinearity between the load and the gap, but also able to control the refining consistency to its setpoint according to the requirement of pulp quality [52]. Also, the use of Laguerre function-based approaches can eliminate the requirement for the accurate knowledge of the model order, time delay and process dynamics. Constrained Predictive Control For the purpose of control applications, it is always necessary to take process con-straints into account in a controller design, due to safety, equipment physical limitations and product quality requirement. The controller should not be allowed to drive input-output variables outside their specified bounds. Du and Dumont [53] first proposed a constrained control approach for the control of a wood chip refiner. In the use of this approach, a wood chip refiner was modelled as a 3 x 2 system. A constrained predictive Chapter 1. Introduction 14 control law based on G P C was derived by taking constraints into account. Optimum solution was calculated through quadratic programming. A n adaptive control scheme was used to handle the time-varying and nonlinear natures of the process. The simula-tion results showed [53] that the proposed controller was able to maintain the outputs at their setpoints without violating the constraints. It was also shown that in the use of the predictive control strategy, input-output constraints can be handled explicitly and efficiently. 1.2.4 Control Problems Extensive research effort towards better control of a T M P plant has been made over the last two decades. However, the control of today's refiner mechanical pulping plant is still relatively primitive [33]. This may be due to three major problems: (i) a lack of on-line sensors for measuring key process variables, (ii) unknown process mechanism and (iii) limitations of the past control methods. A lack of sensors On-line sensors are the key to successful control of industrial processes, but they are not available sometimes. Sensors required for the control of refiner mechanical pulping can generally be divided into three groups: sensors for raw material, sensors for refining conditions and sensors for pulp qualities. The main difficulty in the control of specific energy is our inability to measure fibre feed rate or throughput accurately. For a given chip transfer screw speed, throughput variations are mainly resulted from changes in inlet chip bulk density and moisture con-tent. Even when throughput is constant, pulp quality may still vary because of changes in wood species. Different wood species will require different energy input due to differences in fibre length, coarseness, etc. The variations in raw material should be compensated in order to produce pulp with the desired quality. However, due to a lack of on-line sensors Chapter 1. Introduction 15 neither these variations nor their effects on the pulp quality can be measured accurate and fast enough at present [39, 54]. Sensors for refining conditions are also required to help stabilize refiner operation and compensate for immeasurable changes. Consistency in the refining zone is perhaps the most important variable in determining the quality of the pulp [55], but it is difficult to measure the consistency in high consistency refining. Presently, several blow-line consis-tency sensors are available[56, 41, 40] with varying degree of success in their applications. There has been no publication assessing the performance of these sensors. This may be partially due to the fact that these sensors are all relatively new and their reliability needs to be investigated. Refining intensity is also an important variable in determining the quality of the pulp for a given specific energy. If there was a technique available for on-line measuring the refining intensity, there would be one more degree of freedom in the control of refining and pulp quality. On-line sensors for measuring pulp properties are absolutely necessary in manufactur-ing high quality pulp, but there is still a lack of fast, continuous, reliable on-line sensors for pulp quality [33]. A relatively sophisticated analyzer: the Pulp Quality Monitor or P Q M [57, 25, 58] is currently available. It can be used for measuring pulp drainage, fibre length distribution and shive content. The P Q M has the distinct advantage of being ca-pable of performing several tests per hour in comparison to the relatively low frequency manual testing [59]. However, the P Q M has not been widely accepted in North Amer-ica because of its high cost and maintenance requirements [33]. Relatively simple and inexpensive quality analyzers are needed by many mills. Unknown Process Mechanism The main impediment to improving control of refiner mechanical pulping process is a lack of understanding of refining mechanisms. Over the last few decades, tremendous efforts have been made by many researchers and engineers towards a better understanding Chapter 1. Introduction 16 of the wood chip refining process. Despite intensive activity in the field and improved knowledge of the refining process, there is still a lack of understanding of the refining mechanism in a disc refiner. More theoretical and experimental work are needed in the fields of the mechanism of energy transfer and mass flows in the refining zone during wood chip refining. Limitations of Control Methods A traditional method for the control of pulp quality is based a quality-energy rela-tionship. Associated controllers are single-input single-output linear model based with fixed controller parameters. In many cases, control actions are made and implemented on the assumption that the resulting control actions are feasible and can be implemented, However, an actual process is mainly a function of two variables: the input energy per unit of mass pulp or specific energy, and the specific energy per refiner bar impact or refining intensity. The process is characterized by nonlinearity, variable interactions, pro-cess constraints, time varying nature and external disturbances. Therefore, conventional control methods are not satisfactory and limited by their control performance. 1.3 Objective of this Study In the past, various control strategies have been proposed for the control of wood-chip refining, but the associated control performances are limited for certain reasons. First, past control methods were more focused on controlling the process as a single-variable, time-invariant and linear process. They were not satisfactory because the real process is a stochastic one with multivariable, nonlinear and time-varying characteristics. Secondly, in many cases, past control approaches were developed based on a quality-energy rela-tionship, i.e., pulp quality is determined by specific energy alone. However, pulp quality is a function of two variables: the specific energy and the rate of the transferred specific Chapter 1. Introduction 17 energy or the refining intensity. For a given specific energy, pulp quality will vary as the refining intensity changes. Thirdly, the control calculations in past control applications were made and implemented with the assumption that the resulted control actions were feasible and could be implemented. This would limit the control performance because the inability to implement the control signal exists due to the reasons of safety, equip-ment physical limitations and product quality requirements. The limitations of the past control approaches motivated the current study. The objective of this study was to develop a control strategy suitable for the control of an industrial T M P plant. To achieve this goal, three tasks were completed. The first of these was to model the process using mechanistic and empirical methods. The second task was to develop a constrained multivariable control strategy for the control of the T M P plant. The third was to investigate the proposed control strategy on a simulation model developed using the mechanistic and empirical models to approximate the characteristics of an industrial situation. 1.4 Contributions of this Thesis This thesis describes the development of a novel control strategy for the control of a two-stage T M P plant. Desired pulp quality is maintained by selecting specific energy and refining intensity setpoints which are achieved through the control of production rate, motor load and refining consistency at each stage. The manipulating variables are hydraulic closing pressure on refiner plates and flowrates for both chip feed (at the lst-stage) and dilution water. A constrained multivariable controller is developed based on the G P C algorithm because of its simplicity, flexibility and capability of handling prob-lems in one algorithm. A n optimal solution of the constrained M I M O G P C is obtained by solving a quadratic programming or mixed-weights least-squares problem. A M I M O Chapter 1. Introduction 18 C A R I M A model is used for the controller development. However, in the case of process severe nonlinearity, i.e., the incremental gain between refiner motor load and plate gap is subject to a sudden change in the sign, a mixed linear-nonlinear Laguerre model is used for the controller design. The major contributions of this thesis are as follows: • A steady state model is developed to link process input variables to output variables. The outlet consistency at each stage is modelled using mass and energy balances, showing that outlet consistency is a function of production rate, motor load and inlet consistency. If the measurements of refiner inlet and outlet consistency are available, the consistency model may be used to back calculate production rate. • Process dynamics and disturbances are modelled through identification experiments performed on an industrial process. The process is identified as a multivariable dynamic-stochastic process with lower-order dynamics plus time delay. In the use of process identification (or the empirical method) rather than the mechanistic method, process dynamic modelling can be realized more easily without the need for the development of complex differential equations. The proposed control scheme will be based on the model structure obtained from process identification stage. • A multivariable control strategy is developed for the control of the T M P plant. Desired pulp quality is achieved through selecting specific energy and refining in-tensity setpoints which are achieved through the control of production rate, motor load and outlet consistency at each stage. The manipulated variables are closing pressure on refiner plates and the flowrates for both chip feed and dilution water. The proposed control strategy will have more degrees of freedom in the control of the process, and the capabilities of optimizing process operation and improving product quality. Chapter 1. Introduction 19 • A constrained multivariable controller for T M P is developed based on the G P C algorithm for the control of the process subject to input and/or output constraints. A n optimal solution of the constrained controller is obtained through solving a quadratic programming (QP) problem. Simulations demonstrate that the control law given by input and output constraints is realizable. In the case of an unfeasible constrained control problem, an optimal solution is obtained by solving a mixed-weights least-squares (MWLS) problem and the M W L S will converge to a point that minimizes the maximum constraint violation. In addition, in the use of the M W L S the control performance index for any case of constrained control problems can be augmented easily and the weights can be modified in a manner that encompasses both the requirement for the future control movements to lie inside the feasible region and to minimize the control performance index. Computer simulations show that the performance of the constrained controller is better than the unconstrained one. • For the case of the severe nonlinearity of the process, i.e., the incremental gain between refiner motor load and plate gap is subject to a sudden change in the sign, Laguerre functions - a type of orthonormal functions - are considered as an alter-native to represent plant dynamics and nonlinearity. The multivariable predictive controller is then derived based on a mixed linear-nonlinear Laguerre model. The nonlinear Laguerre functions are used to represent the severe nonlinearity between refiner motor load and plate gap, while the linear Laguerre functions are used to approximate the relationships between other inputs and outputs. The major ad-vantages of the Laguerre approximation are: (i) easy to model and reconstruct, (ii) flexible model structure and simple representation, (iii) understandable structure similar to transient signals, and (iv) the capability of representing a plant properly Chapter 1. Introduction 20 without the need for assumptions about plant order and time delay. Computer sim-ulations demonstrate the robustness and efficiency of the proposed control scheme, along with the ability to avoid plate clashes due to sudden changes in the sign of the incremental gain between refiner motor load and plate gap. • The efficiency and applicability of the proposed control schemes are tested on a simulation model developed using a combination of the empirical and mechanistic methods to approximate an industrial situation. The simulation results demon-strate various advantages of the proposed control schemes with the potential to offer an automated control of mechanical pulp manufacturing. 1.5 Outline of the Thesis After a brief review of the fundamentals related to refiner mechanical pulping processes, Chapter 2 describes static modelling of the process using the mechanistic method. A qual-itative relationship between raw material quality, operating conditions and pulp property is presented. Simulations show the variable interactions and nonlinearity of the process. Chapter 3 presents dynamic modelling of the process using the empirical method. To avoid the complexity of deriving differential equations and the inability of modelling process disturbances through the mechanistic modelling, process identification is used to obtain the characteristics of process dynamics and disturbances. Chapter 4 describes the development of a constrained multivariable predictive controller for the control of a T M P plant. A n optimal solution of the constrained M I M O controller for simple cases is derived by solving a quadratic programming problem, while for general cases an optimal solution is derived by solving a mixed-weights least-squares problem. Chapter 5 presents the development of M I M O Laguerre model-based controller to avoid system oscillation resulting from a linear-model-based controller in the case of severe nonlinearity of the Chapter 1. Introduction 21 process. Using the Laguerre model-based control approach, the dynamics of the process can be represented appropriately without the need of assumption of the model order and time delay. The proposed control strategy for the control of the two-stage T M P plant is detailed in Chapter 6. Desired pulp quality is achieved through selecting specific energy and refining intensity setpoints which are obtained by changing production rate, motor load and outlet consistency. Simulation results demonstrate the efficiency and applica-bility of the proposed control schemes. The thesis conclusion and future work are given in Chapter 7. Chapter 2 Process Variables, Modelling and Simulations 2.1 Introduction A mathematical model, generally speaking, can be a completely mechanistic one, an empirical one or a combination of these two. A mechanistic model based on physical principles (or the first principles) is more desirable since it is generally applicable over a wide range of operating conditions and gives more insight into a process, but it may be time consuming to develop. By comparison, an empirical model derived by fitting a mathematical equation to a set of experimental data is easier to develop and possesses more simplicity, but because it approximates the process over a limited range of opera-tion, it will present the real process with some limitation. Which approach is chosen for modelling depends on the purpose of the model. In the past, much work has been done on modelling refining process using different approaches. Most of the work focused on modelling the process in order to obtain more insight into the mechanism of the process, rather than for control purposes. The mathematical model developed in this study is for control purposes and will give a trade-off between the model accuracy and simplicity. The model is developed by using a combination of mechanistic and empirical methods. Such a mechanistic-empirical model not only gives some insight into the process mechanism, interactions, and nonlinearity, but also represents the characteristics of process dynam-ics and disturbances. Mechanistic modelling is given in the following, while empirical modelling will be presented in Chapter 3. 22 Chapter 2. Process Variables, Modelling and Simulations 23 This chapter is organized as follows. Section 2.2 lists a number of main process variables which will be considered in process modelling, such as manipulated variables, operating variables, pulp quality variables and disturbance variables. Section 2.3 presents modelling of a two-stage refiner mechanical pulping process. The developed model will show process interactions and nonlinearity. Section 2.4 includes simulations of the model and result discussions. Summary and conclusions are given in Section 2.5. 2.2 Process Variables Refining is a series of interrelated events designed to produce pulp. The interrelationships of these events can be expressed in terms of process variables. The variables, under this study, are divided into four groups: manipulated variables, operating variables, pulp quality variables and disturbance variables. Following are the key variables which have strong effects on wood chip refining. Other process variables such as steam flow, refiner plate pattern design, chip feed temperature, dilution water temperature, etc. can also affect chip refining, but they are not included ether because the effects of the variables are small or because the variables are relatively constant during refining. 2.2.1 Key Manipulated Variables Manipulated variables are the input variables which can be adjusted during the process operation to control the refining conditions. Changes in manipulated variables cause changes in the operating conditions and thereby changes in the end product quality. Chip Transfer Screw Speed The chip transfer screw, located prior to the primary refiner, is used to adjust the volumetric flow rate of wood chips to the process. Currently, many mills use transfer screw speed to set desired production rate or throughput. Any changes in the screw Chapter 2. Process Variables, Modelling and Simulations 24 speed will change dry fibre flow to the process. As a result, the energy input per unit dry fibre or specific energy will be varied. Hydraulic Closing Pressure The hydraulic closing pressure applied by a hydraulic loading system mounted in a refiner is adjusted to position the plate gap of the refiner. A n increase (or decrease) in the hydraulic pressure results in a decrease (or increase) in the gap and hence an increase (or decrease) in motor load because motor load is more sensitive to the gap. Dilution Water Flow Rate The consistency in refiner refining zone is important to pulp quality. For a given spe-cific energy, any change in the consistency may alter the rate of the specific energy input or refining intensity and as a result alter pulp quality. Also, excessive consistency varia-tions can lead to unstable refiner operation (plate gap, power, etc.) [28]. Manipulating the dilution flow at the inlet of each refiner is the most convenient way of maintaining the consistency in the refining zone. Disc Rotating Speed Disc rotating speed is defined as the number of disc revolutions per minute. Any change in rotating speed will directly vary the centrifugal force on the pulp in the refining zone and thus change the residence time for a fibre to pass through the refining zone. Any change in residence time will alter the number of refiner bar impacts per unit pulp and for a given specific energy, this will in turn lead to changes in specific energy per bar impact or refining intensity. As a result, pulp quality will be changed. Generally, disc rotation speed is fixed once it is determined according to production requirement and refiner physical reliability. However, strong effect of disc rotating speed on the refining intensity suggests that manipulating disc rotating speed, at a given specific energy, will be a more sensitive way of adjusting pulp quality. Chapter 2. Process Variables, Modelling and Simulations 25 2.2.2 Key Operating Variables The quality of the pulp produced depends very much on the operating conditions in the refining zone between plates. The operating conditions in the refining zone will change with changes in the manipulated variables. Motor Load A motor integrated with a refiner is a power actuator which delivers electrical energy to the refiner. Motor load or net power is defined as the electrical energy applied to the fibre flow in the refiner. At a given throughput, any change in motor load may vary pulp quality because motor load is proportional to specific energy. Motor load may vary with changes in throughput, plate gap, refining consistency, raw material quality, etc. Refining Consistency Refining consistency is important in a refiner mechanical pulping process, because it changes the energy-quality relationships [60]. At a given specific energy, refining at different consistencies results in the pulp with different quality [61]. Refining consistency may change with changes in throughput, dilution flow, motor load, etc. but the dilution flow is the most convenient way to be used to adjust the consistency. Plate Gap Plate gap is the separation of two discs of a refiner and is controlled by either a hydraulic or electro-mechanical loading system integrated into the refiner. The gap mea-surement can be obtained by a gap sensor in the refiner or indicated by the relative shaft position of the loading system [49]. In practice, the gap measurement is used for indi-cating the plate movement and/or interlocking purpose to prevent plate clashing, rather than for closed-loop refiner control. Any changes in plate gap will change the mechanical force exerted by the wood material in the refining zone on the plates, and as a result alert motor load and specific energy for a given throughput. Chapter 2. Process Variables, Modelling and Simulations 26 Production Rate Production rate or throughput is defined as oven-dry 1 wood flow fed to a refiner. Many mills use chip screw speed to set desired production rate. However, this is not accurate due to variations in chip density and moisture content. For a given motor load, any variations in throughput may result in changes in the energy input per unit fibre or specific energy and thus pulp quality. Specific Energy It is known that pulp quality is strongly affected by specific energy. Since specific energy is proportional to motor load and inversely proportional to throughput, the tra-ditional method of adjusting pulp quality is through manipulating the motor load or the throughput. Refining Intensity At a given specific energy, refining intensity determines the rate of the energy inputs. For a given specific energy, different refining intensity will produce the pulp with quite different quality [17, 62, 63, 10, 11]. It is only in the last few years [11] that the concept of refining intensity has been identified and its importance has been realized. It has been suggested [64] that refining intensity must now be recognized as an important variable in all types of refining. Refining intensity is a function of specific energy, residence time, rotating speed and refiner design. However, for a given specific energy residence time is only available way to alert the refining intensity. Residence Time The residence time for a fibre to pass the refining zone is the most important parame-ter in the determination of the final pulp quality for a given specific energy [10, 15]. This is because changes in the residence time will bring about changes in the number of bar 1 oven-dry: condition of cellulosic material that has been dried to constant mass (or weight) at a temperature of about 105 °C. Chapter 2. Process Variables, Modelling and Simulations 27 impacts received by a fibre. As a result, average specific energy per impact or refining intensity will change if specific energy is constant. Residence time is a function of specific energy, inlet consistency, rotational speed and refiner design. For a given specific energy, manipulating inlet consistency is the most convenient way to alert the residence time. 2.2.3 Key Pulp Quality Variables The pulp properties in terms of papermaking like freeness, fiber length, fiber specific surface, shive content, coarseness, flexibility, strength, etc. are important. The choice of different property variables depends on the different end users for the different pulp quality and also relies on the techniques available for obtaining the measurements of these variables. In the following, the Canadian Standard Freeness, long fibre content and shive content are considered as property variables. Such choice is not only because these variables can be used to predict handsheet strength and pulp drainage properties, but also because they can be measured using the techniques available on the market, such as Pulp Quality Monitor (PQM). P Q M is an automated device for on-line measuring freeness, fibre length distribution and shive content. Pulp quality is dependent on the conditions of refining and can be consequently adjusted by changing the operating conditions through manipulating the manipulated variables. Canadian Standard Freeness The degree of refining done on pulp is measured by the drainage of pulp, called free-ness. However, pulp drainability is usually obtained by means of the Canadian Standard Freeness (CSF) test, which is defined as the number of ml water collected from the side orifice of the standard tester when a dilute stock drains through a perforated plate under carefully controlled conditions [65]. Since freeness characterizes the strength of mechanical pulp, it is used to describe the strength of the produced pulp [42]. Long Fibre Content Chapter 2. Process Variables, Modelling and Simulations 28 Long fibre content is another variable used to describe the strength of the produced pulp. The longer the individual fibre is, the stronger the pulp is. The long fiber content is defined as the percentage by weight of fibres retained on the 48-mesh screen of the Bauer-McNett fiber length classifier [66]. Shive Content Shive content is defined as the percentage of oven-dry pulp retained on a standard slotted fractionating plate, usually with slit width of 0.15mm [65]. It is possible for pulps with identical long fibre content to have different strength and drainage properties due to the presence of the shives. So shive content is also used to indicate pulp strength and drainage properties. 2.2.4 Disturbance Variables Disturbance variables are the variables which cannot normally be controlled and often cannot even be measured. The variations in wood quality and refiner plate conditions are the normal disturbances [67] influencing refining conditions and pulp quality. Raw material quality Very often, a major disturbance to a refiner mechanical pulping process is the varia-tions in raw material quality, such as wood species, chip size, chip bulk density and chip moisture content. However, chip bulk density variation is a main cause of the disturbance [54]. Any change in the chip bulk density will change the flow rate of the dry wood mass or the throughput. For a given motor load, this will lead to variations in the specific energy and hence in pulp quality. Also, variations in the chip density reflect changes in the nature of the wood which may directly affect the pulp quality [61]. Refiner plate wear Refiner operation may not have a stationary behavior partially due to refiner plate wearing. Plate wearing normally occurs gradually over a period of several hundred hours Chapter 2. Process Variables, Modelling and Simulations 29 due to erosion and corrosion, but it can also occur suddenly in the case of plate clash-ing [30]. For a given specific energy, new plates tend to produce the pulp with higher freeness than worn plates. When the plates become extremely worn, short-term load variations increase and it becomes increasingly difficult to fully load a refiner. Plate wearing cannot be avoided and often not even measured, but it can be compensated for by selecting some advanced control techniques. DISTURBANCES: Wood species Chip size Chip moisture . Chip density Chip temperature Plate age MANIPULATED VARIABLES: Dilution flow Hydraulic pressure Transferscrew speed Disc rotation speed OPERATION VARIABLES: Motor load Production rate Outlet consistency Specific Engergy Refining intensity Pulp OPERATION VARIABLES: Motor load • Plate gap Outlet consistency Specific engergy Refining intensity QUALITY VARIABLES: Freeness Long fibre content Shive content Ist-stage refiner 2nd-stage refiner MANIPULATED VARIABLES: Dilution flow Hydraulic pressure Disc rotation speed Figure 2.4: Interactions between manipulated variables, operating variables, pulp quality variables and disturbances for a two-stage refining 2.2.5 Variable Interactions Interactions between the variables are given in Figure 2.4. As can be seen, the final pulp quality is determined by the refining operating conditions and the raw material quality. Since variations in raw material quality cannot be controlled and often cannot even be Chapter 2. Process Variables, Modelling and Simulations 30 measured, desired pulp quality can only be achieved by setting correct operating condi-tions. Changes in wood quality and operating conditions can be compensated through adjusting manipulated variables. 2.3 Modelling a TMP Process The purpose of this section was to develop a steady state relationship between the vari-ables described previously. The relationship developed not only gives some insight into process interactions and nonlinearity but also provides a valuable tool for process quanti-tative analysis. The variables appearing in the following derivation can be found in List of Symbols in the beginning of the thesis. 2.3.1 Modelling Throughput Throughput or production rate is denned as the oven-dry wood flow fed to a refiner. Wood chips delivered by the chip transfer screw are usually metered on a volumetric basis. If assuming that the chip transfer screw is fully loaded with chips, the volumetric flow rate Q of the wood chips (dry wood + moisture) to the process may be expressed as a function of the chip velocity v in the screw direction: where A is the screw cross sectional area available to the chip flow. Generally, the chip velocity v in the screw direction is proportional to the chip screw rotational speed Sc, i.e. where A: is a constant factor, determined by the size of the screw systems. If assuming that the chip bulk density pc of the incoming wood is known, then the wood chip flow Q = A-v (2.1) v — k • Sc (2.2) Chapter 2. Process Variables, Modelling and Simulations 31 rate fc can be expressed by: fc = Pc-Q (2.3) If assuming that the solid content s of the incoming chips is known, the oven-dry (moisture free) chip mass flow rate or throughput Fc is: Fc = 8-fe (2.4) = s- pc-Q Substituting Equations 2.1 and 2.2 into Equation 2.4 gives: Fc = s • pc • A • k • Sc (2.5) The above equation shows that for given chip screw systems, throughput or production rate Fc is a function of chip bulk density p c , chip solid content s and chip screw speed Sc. Generally, chip bulk density and solid content are not constant. Chip bulk density varies with changing in wood density and chip size, whereas solid content changes with changing in wood species and age of the chips. In practice, the chip transfer screw speed is the only available way to alter throughput or production rate. 2.3.2 Modelling Consistency The consistency at the refiner outlet is derived based on the conservation principle of two fundamental quantities: the total mass and the total energy. The lst-stage refining consistency The mass and energy flows in the primary refiner are illustrated in Figure 2.5. The energy balance in the primary refiner is: Chapter 2. Process Variables, Modelling and Simulations 32 motor load forward steam where backward steam -wet chip flow dilution flow wet pulp flow seal water Figure 2.5: Mass and energy flows in the lst-stage refiner energy in = energy out + energy losses energy in = ec + ed + eaw + r/P energy out = esl + es2 + ep energy losses = 0 (assuming negligible energy losses) — Fc • cc ' Tc -f- fwi ' cw ' Tc (2-6) (2.7) — Fd • cw • Td esi =val-L eS2 = vs2-L £p — ff ' Cf ' Tp -\- f w 2 ' Cw ' Tp where u s l is the backward steam flow rate, vs2 is the forward steam flow rate, and L is the enthalpy of the steam. Note that the temperatures given in Equation 2.7 are the relative temperatures from a reference temperature where the enthalpy of the liquids is assumed to be zero [68]. Chapter 2. Process Variables, Modelling and Simulations 33 The mass balances in the primary refiner are: water in = water out + water losses water in = fwl + Fd water out = fw2 + v3l + vs2 water losses = —F s t u(seal water and leakage) fibre in = fibre out + fibre losses fibre in = Fc fibre out = ff fibre losses = 0 (assuming 100% yield pulp) and the water flow rate fw\ carried by the inlet chips is calculated by: (2.8) /»! = — - Fc (2.9) s (2-10) It can generally be assumed that: c c = Cf Substituting Equations 2.7, 2.8, 2.9 and 2.10 into Equation 2.6 gives: nP - [Fc cc(Tp - Tc) + (Fc/s - Fc)cw(Tp - Tc) + (Fsw + Fd)cw(Tp - Td)} V s l + V a 2 = ( L _ c T p ) (2.H) Consistency is defined as the ratio of dry fibre flow to the dry fibre flow plus water flow. So the outlet consistency C of the primary refiner is expressed as: C = -^j- (2.12) // + Jw2 Substituting Equations 2.8, 2.9 and 2.10 into 2.12 gives: F (j = i i (2 13) ^ + Fd + Fsw - {val + vs2) Chapter 2. Process Variables, Modelling and Simulations 34 The above equation shows that the outlet consistency is a function of the dry chip flow Fc, chip solid content s, dilution flow rate Fd, seal water rate Fsw and total steam flow vsi + vS2- Substituting Equation 2.11 into Equation 2.13 yields: In general (2.15) L n P » Fc cc(Tp - Tc) + (Fe/s - Fc)cw(Tp - Tc) + (Fsw + Fd)cw(Tc - Td) Equation 2.14 can be then simplified as: F L ^ ~ ^fL + FdL + FSWL — nP ( 2 " 1 6 ) Similarly, the inlet consistency c, is defined as the ratio of the dry wood flow rate to the total mass flow at inlet of a refiner: Ci = — (2.17) *f + Fd + Fsw K ' The 2nd-stage refining consistency The consistencies in the 2nd stage refiner are derived in a similar fashion as were the consistencies in the 1st stage. The outlet consistency in the 2nd-stage refiner is given as: F T Q C (2 18") Yr^L> + Pd^i + r s w L , c w v p 1 L where Tj (°C) is the temperature of the inlet pulp and m (%) is the moisture content of the inlet pulp. If considering the conditions given in Equation 2.15, the above equation can be further simplified as: F T ^ . L + FdL + FswL-nP V " ; Chapter 2. Process Variables, Modelling and Simulations 35 The inlet consistency in the 2nd stage is calculated as: Ci = -p § — (2.20) -ts. L TP, -L. Jp V ' 1 - m ' r d * r a w Note: • The consistency equations derived above may be used to predict inlet or outlet consistency for each stage. If a consistency sensor is available and if significant errors between the measurements and model estimations exist, the errors should be taken into account in sensor calibration or consistency model improvement. • The relationship between the inlet and outlet consistencies, e.g. in the lst-stage, can be easily derived by substituting Equation 2.17 into Equation 2.16 as: C - T j ^ m <2:21> If consistency measurements are available, the above equation can be used to back calculate the production rate Fc, i.e. < 2 - 2 2 > Static modelling of the outlet consistency can also be found in [69, 70], but their models were derived with more complexity or unknown parameters. The consistency model given in [16] was derived based on the radial velocity of the pulp in the refining zone from a calculation of the forces that govern the flow of the pulp. However, simulations of the consistency model (Equation 2.16) based on mass-energy balances predict the same results as that given in [16]. 2.3.3 Modelling Energy Input The degree of refining, a critical factor to the final product, is strongly influenced by the energy input to the pulp flow in wood chip refining. Chapter 2. Process Variables, Modelling and Simulations 36 Specific energy (E) is defined as the ratio of the net power or the motor load P to the throughput or production rate Fc: E=§- (2.23) The above equation shows that the specific energy E is proportional to the motor load P and inversely proportional to the throughput Fc. Therefore, the motor load and/or the throughput can be adjusted according to the desired specific energy. M o t o r load (P) is defined as the electrical energy applied to a refiner. In a pressur-ized refiner, the total axial thrust force Ftotai applied on refiner plates is balanced by the steam reaction force Fs acted on the plates and the wood reaction force Fw exerted on the plates [71, 3, 72], i.e. Ftotai = FW + FS (2.24) The experimental study by Miles and May shows that the motor load P for each refiner linearly depends on the wood reaction force Fw which is proportional to the closing pressure Pc of the hydraulic loading system. So we have: P = a-Fw (2.25) and Fw = (3-Pc (2.26) where a and /3 are the coefficients, a may be determined by the refiner size, disc rotational speed, refining zone consistency etc., whereas /? is mainly determined by the size of the hydraulic loading system [71, 3, 72]. Combining the above two equations gives: P = a-f3-Pc (2.27) Many experimental studies show that the motor load is also proportional to the through-put, i.e. P = k-Fc (2.28) Chapter 2. Process Variables, Modelling and Simulations 37 This is because changes in the throughput will vary the wood reaction force on the plates and thereby alter the motor load. Due to the complex relationship, little work has been done towards deriving the static quantitative relationship between the throughput, closing pressure and motor load. To start with a simple model, this relationship may be considered as: For control design purposes, the above model can also be linearized about some operating point and modified as: where P is the motor load deviation from its operating point, Pc and Fc are the deviations from their corresponding steady-state values, and kx and k2 are the coefficients which will be determined by the operating conditions, plate conditions, refiner design and raw material quality. Residence t ime (r) is defined as the time a pulp takes to pass through the refining zone. A theoretical model of the residence time has been derived by Miles and May [16] through calculation of the radial velocity of pulp in a high-consistency chip refiner. Through this model, the residence time is linked to a number of variables in terms of major operating variables and refiner design parameters as follows: P = kPcFc (2.29) P = hPc + k2Fc (2.30) L8-L s CiE (2.31) Chapter 2. Process Variables, Modelling and Simulations 38 where w : disc rotational speed Cj : inlet consistency r i , r 2 : the inlet and outlet radii of the refining zone fir,Ht '• radial and tangential friction coefficients between pulp and discs La : latent heat of steam a : 4 for a single-disc refiner, 2 for a double-disc refiner As can be seen from the above equation, the residence time r is a function of refining consistency, disc rotational speed, specific energy and refiner design. Although the effect of rotational speed is far more significant, refining consistency, for a given specific energy, is the most readily available means of altering residence time in industrial operations. Ref ining intensity (/) is defined as the average specific energy per refiner bar impact [16]: (2.32) where N is the number of refiner bar impacts received by unit pulp and is calculated by [55]: N = nh W^Y^- T (2.33) where n = the number of bars per unit length of arc h = 1 for a single-disc refiner, 2 for a double-disc refiner. It can be seen from Equation 2.32 that for a given specific energy E, the number of bar impacts N governs the average amount of energy transferred at each impact, i.e. the refining intensity / . The concept of refining intensity implies that refining will depend not only on the specific energy applied but also on the intensity of its application. Equations 2.32 and 2.33 show that for a given specific energy appropriate manipulation of the Chapter 2. Process Variables, Modelling and Simulations 39 residence time r by adjusting inlet consistency will provide a desired combination of N and / . Specific refining power (e) is denned as the rate of energy delivered per unit mass of pulp [6]: E e = — (2.34) r Equations 2.32 and 2.34 can be combined as a ratio: e N (2.35) or , {n + r 2 ) . . e = nhw1—-—-I (2.36) As can be seen from the above equation, the specific refining power is another measure-ment of the severity or intensity of refining for a given specific energy. Energy input over two stages. The quality of final product is dependent not only on the energy input at each stage but also on total energy split and refining intensity distribution over two stages. The total motor load, specific energy, residence time and number of bar impacts can be expressed as: Ptotal = P\ + ?2 Etotal — Ei + E2 (2.37) Ttotal = Ti+T~2 Ntotal = Nt + N2 where the subscripts 1 and 2 refer to the first and the second stage respectively. The total refining intensity Itotai and total specific refining power etotai are calculated by: Itotai = § ^ (2.38) ^total etotai = — (2.39) ^total Chapter 2. Process Variables, Modelling and Simulations 40 2.3.4 Modelling Pulp Properties Recent experimental work by Qian and Tessier [70] has shown that pulp properties such as freeness CSF(ml), long fibre content LF(%) and shive content SC(%) are mainly dependent on the specific energy E, refining intensity / and specific refining power e. The empirical relationships between these variables are expressed as [70]: CSF = (CSF0 - h(E - Eb))10 f c a ( / - 7 o ) LF = LF0 + k3(e - e0) + k4(E - E0)) (2.40) SC = SC010^k5(-E^Eo^+k6^~^ where CSF0, LFQ and SCQ are the initial values of pulp properties. E0,10 and e0 are the initial conditions of operating variables. A l l the coefficients k{ (i = 1 • • • 6) are dependent on the type of refiner and can be determined or updated if on-line measurements of pulp properties and operating variables are available. As can be seen from Equation 2.36, the specific refining power e can be determined from the refining intensity / . Therefore, the pulp quality in Equation 2.40 is mainly a function of the specific energy E and the refining intensity / . 2.4 Simulation Results and Discussions The objective of simulations was trying to demonstrate process interactions, nonlinearity and input-output quantitative relationship. Simulations were performed on the math-ematical relationship illustrated in Figure 2.6, which was proposed in Section 2.3. In order to approximate features of an industrial plant, simulations were performed under the operating conditions of an industrial C T M P production fine. Associated operating conditions and refiner design parameters are given in Tables2 2.1 and 2.5. 2 see references [73] and [70] Chapter 2. Process Variables, Modelling and Simulations 41 MPc<Fc) f2(Sc,Pc,s) F d Pc s f3(Fd,Fc,s) f4(Fd>Fc,P,s) T s f5(E.cr,a>] CO f7(E-N> M M e empirical model CSF > LF SC Figure 2.6: Model structure of the lst-stage wood-chip refining 2.4.1 Simulating the Primary Refiner Effects of hydraulic closing pressure and chip transfer screw speed Figure 2.7 shows for a given dilution flow rate (Fd = A.lkg/sec) the effects of the first stage closing pressure Pc and chip transfer screw speed Sc. The dashed line indicates Sc = 25 rpm (Fc = 271 t/d), the dashdot line indicates 5 C = 27.5 rpm (Fc = 298 t/d) and the solid line indicates Sc — 30 rpm (Fc — 325 t/d). In practice, an increase in the closing pressure for a given throughput will increase the motor load and hence the specific energy. This has been indicated in Figure 2.7. The reason is that increased closing pressure increases the total axial thrust force applied to the refiner plates and, as a result, increases the motor load which is proportional to the total thrust. The outlet consistency (Figure 2.7) increases with an increase in the closing pressure because a higher motor load resulting from increased closing pressure brings about a higher temperature in the refining zone and, therefore, more steam and less wet mass in the refining zone. Chapter 2. Process Variables, Modelling and Simulations 42 Table 2.1: Primary Refining Conditions Variables Typical Value Throughput, Fc 325 a.d.t/d Specific Energy, E 3.4 G J / t Motor Load, P 12.5 M W Dilution Water, Fd 280 1/min Seal Water, Fsw 75 1/min Disk Rotation Speed, w 1800 r/min Chip Screw Speed, Sc 30 r/min Chip Moisture Content 45.6±3.2 % Chip Solid Content, s 54.4 % Chip Bulk Density, pc 98.9±5.8 kg /m 3 Closing Pressure, Pc 1100 psi Refining Pressure 290 psi Inlet Radii of Refining Zone, r\ 51.4 cm Outlet Radii of Refining Zone, r2 76.2 cm Bars per Unit Length of Arc, n 100 Single Disc Refiner a=4, h = l As can be seen from Figure 2.7, an increase in the screw speed (or the throughput) increases the motor load and the outlet consistency. This is because increased throughput brings about a higher wood reaction force acted on the plates and hence a higher motor load. The higher the motor load is, the higher the temperature in the refining zone is, and as a result the less the wet mass and the higher the consistency are. As the specific energy is proportional to the motor load and inversely proportional to the throughput, the direction of changes in the specific energy is dependent on both the motor load and throughput. At a given closing pressure, the specific energy increases is due to a larger increase in the motor load, while the specific energy decreases is due to a larger increase in the throughput (see Figure 2.7(b)). Table 2.2 summarizes the simulation results given in Figure 2.7. For a given throughput 325 t/d, 40% increase in the closing pressure (from 1000 to 1400 psi) results in about 63% increase in the load (Figure 2.7(a)), 64% increase in the specific energy (Figure 2.7(b)) Chapter 2. Process Variables, Modelling and Simulations 43 and 36% increase in the outlet consistency. Table 2.2 also shows that for a given closing pressure=1200 psi, 20% increase in the throughput (from 271 to 325 t/d) brings about 25% increase in the motor load, 5% increase in the specific energy, and 17% increase in the consistency. The simulation results show that the motor load seems more sensitive to the closing pressure changes than to the throughput. Table 2.2: First Stage Refining at Fd = 4.7kg/sec Closing Pres.=1000-1400psi Throughput=271-325 t /d A Closing Pres.=40% AThroughput=20% (throughput=325t/d) (closing pressure=1200psi) Motor Specific Outlet Motor Specific Outlet Load Energy Cons. Load Energy Cons. (MW) (MJ/kg) (%) (MW) (MJ/kg) (%) 9.5 to 15.5 2.5 to 4.1 42 to 57 10 to 12.5 3.2 to 3.35 41 to 48 +63% +64% +36% +25% +5% +17% Effects of dilution water and chip transfer screw speed Figure 2.8 shows for a given specific energy (E=3.4 GJ/ t ) the effects of the dilution flow rate Fd and the chip screw speed Sc. As can be seen, the inlet consistency decreases with increases in the dilution flow rate. For a given specific energy, the outlet consistency decreases with increases in the dilution flow rate for a given throughput. The explanation is that refining under constant specific energy and a given throughput determines a constant load applied. In this case, an increase in the dilution flow rate determines more wet mass in the refining zone and as a result decreases the outlet consistency. Table 2.3 summarizes the simulation results given in Figure 2.8. For a given specific energy E=3.4 GJ / t , increasing dilution flow by 19% (from 4.2 kg/sec to 5.0 kg/sec) results in a 5.0% drop in the inlet consistency and a 7.8% drop in the outlet consistency. For a given dilution flow rate Fd = 4.6 kg/sec, increasing the chip screw speed by 20% (from 25 r /min to 30 r/min) brings about a 9.3% increase in the inlet consistency and a 16.5% increase in the outlet consistency. The simulation results indicate that the consistencies Chapter 2. Process Variables, Modelling and Simulations 44 Figure 2.7: The 1st stage refining at F^=4.7 kg/sec: effects of closing pressure Pc and chip screw speed Sc ( - - :SC = 25 rpm (Fc=271 t/d), Sc = 27.5 rpm (F c=298 t/d) and — :SC = 30 (F c=325 t/d)) Chapter 2. Process Variables, Modelling and Simulations are more sensitive to the chip feed or the throughput. 45 4 . 4 4 . S fc> - O d i l u t i o n w a t o r ( k g / s e c . ) d i l u t i o n w a t e r ( k g / s e c . ) Figure 2.8: The 1st stage refining at E = 3.4 GJ / t : effects of dilution flow rate Fd and chip screw speed Sc (- -: 5C=25 rpm (Fc=271 t/d), -.-.: Sc=27.5 rpm (F c=298 t/d) and — : 5C=30 (F c=325 t/d)) Table 2.3: First Stage Refining at E = 3AGJ/t Fd = 4.2 - 5.0 kg/sec Sc = 25 - 30 r /min AFd = 19 % A 5 C = 20% (Se = 30 r/min) (F d = 4.6 kg/sec) Ci C C (%) (%) (%) (%) 30.1 to 28.6 51 to 47 27 to 29.5 42.5 to 49.5 - 5 % -7.8% +9.3 +16.5% Effects of rotational speed and inlet consistency Figure 2.9 shows for a given specific energy E = 3.4 G J / t the influences of the rotational speed w and inlet consistency Cj on operating variables and pulp properties. The solid line indicates Cj = 20 %, dashdot line c{ = 23.5 % and dashed line Q = 27 %. Chapter 2. Process Variables, Modelling and Simulations 46 As can be seen from Figure 2.9, an increase in the rotational speed results in decreases in the residence time and the number of bar impacts, and increases in the refining inten-sity and the specific refining power. The explanation is that increased rotational speed increases the centrifugal force on the pulp flow in the refining zone and thus increases the radius velocity of the pulp. As a result, the residence time decreases and thus the number of refiner bar impacts received by unit pulp decreases. For a given specific energy, the specific energy per impact, i.e. the refining intensity increases due to a decrease in the number of bar impacts. The specific refining power increases as a result from a decrease in the residence time. Increased refining intensity will bring about more fibre cutting. As a result, there are less long fibres and shives but more fines, resulting in a decrease in the pulp drainage or freeness. For a given specific energy (E = SAGJ/t), an increase in the inlet consistency in-creases the residence time and thus the number of bar impacts. As a result, the refining intensity and the specific refining power decrease. This is because increased inlet consis-tency reduces the wet mass in the refining zone, resulting in a decrease in the centrifugal force on the pulp flow in the refining zone and hence increases in the residence time and the number of bar impacts. For a given specific energy, increased number of bar impacts brings about a decrease in the refining intensity, whereas an increase in the residence time decreases the specific refining power. Decreased refining intensity results in less fibre cut and thus more long fibres but less fines. Therefore, the freeness, long fibre content and shive content are increased. Table 2.4 summarizes the simulation results given in Figure 2.9. At a given inlet consistency c, = 27%, a 29% increase in the rotational speed (from 1400 rpm to 1800 rpm) brought about a 58 % drop in the residence time r (from 0.6 sec to 0.25 sec). As a result, the number of impacts N was reduced by about 44 % (from 5500 to 3100) and the refining intensity J increased by about 67 % (from 0.6 x l 0 ~ 3 M J / k g to 1.0 x l O - 3 Chapter 2. Process Variables, Modelling and Simulations 47 residence time (s) no. of impacts ! : r 1 6000 | : > 1400 1500 1600 1700 1800 1400 1500 1600 1700 1800 1400 1500 1600 1700 1800 1400 1500 1600 1700 1800 C S F (ml) LF (%) 1400 1500 1600 1700 1800 1400 1500 1600 1700 1800 rotational speed (r/m) S C (%) 1400 1500 1600 1700 1800 rotational speed (r/m) Figure 2.9: The 1st stage refining at E=3.4 GJ / t : effects of rotational speed w and inlet consistency Cj ( —: Cj = 20%, -.-.: ct = 23.5% and - -: Cj = 27%) Chapter 2. Process Variables, Modelling and Simulations 48 MJ/kg) . The changes in the specific refining power were even more pronounced, which was a 127 % increase from 5.5 MJ/kg/sec to 12.5 MJ/kg/sec. As a result of a 28.6% increase in the rotational speed, freeness reduced by 38%, long fibre content reduced by 23% and shive content decreased by 50%. Table 2.4 shows that at a given rotational speed, e.g. w = 1800, increasing the inlet consistency by about 35% (from 20 % to 27 %) results in a 25% increase in the residence time r , 44% increase in the number of impacts AT and 33% decrease in the refining intensity. A 35% increase in the inlet consistency also brings about a 67% increase in the freeness CSF, 25 % increase in the long fibre content L F and 100 % increase in the shive content SC. The simulation results indicate that for a given specific energy refining conditions and pulp properties are strongly related to the rotational speed and refining consistency. The effect of rotational speed seems to be more apparent. If there were a technique available for on-line adjusting the rotational speed, there would be a more sensitive way for the control of refining process. However, the consistency at present is the only available way to alter pulp properties for a given specific energy. Table 2.4: Effects of First Stage Refining on Pulp Properties, E=3.4 G J / t w = 1400 - 1800 r/min Aw = 29% (ct = 27%) r N I CSF L F SC (sec) ( x l O - 3 MJ/kg) (ml) (%) (%) 0.6 to 0.25 5500 to 3100 0.6 to 1.0 800 to 500 65 to 50 4 to 2 -58% -44% +67 % -38 % -23% -50% Ci = 20 - 27 % Aci = 35% (w = 1800 r/min) T N I CSF L F SC (sec) ( x l O " 3 MJ/kg) (ml) (%) (%) 0.2 to 0.25 2250 to 3250 1.5 to 1.0 300 to 500 40 to 50 1.0 to 2.0 +25% +44% -33% +67% +25% +100% Chapter 2. Process Variables, Modelling and Simulations 49 2.4.2 Simulating the Secondary Refiner Effect of rotational speed and inlet consistency Figures 2.10 and 2.11 show for a given secondary specific energy E2 — 2.5GJ/t, the effects of the 2nd-stage rotational speed w and the inlet consistency Q . The solid line indicates C* = 26 %, dashdot line Cj = 30.5 % and dashed line C* = 35 %. In all simulations, the specific energy and inlet consistency in the 1st stage were set to E\ — ZAGJ/t and cx = 27% respectively which produced the residence time T\ = 0.2753 sec and the number of impacts iVi = 3.3 x 103. Table 2.5: Secondary Refining Conditions Variables Typical Value Specific Energy, E 2.5 G J / t Motor Load, P 9.3 M W Dilution Rate, Fd 90 1/min Seal Water, Fsw 20.4 1/min Disk Rotation Speed, w 1800 r/min Closing Pressure, Pc 957 psi Inlet Radii of Refining Zone, rx 51.4cm Outlet Radii of Refining Zone, r2 76.2cm Bars per Unit Length of Arc, n 100 Single Disc Refiner a=4, h = l Figure 2.10 shows that an increase in the rotational speed (or inlet consistency) results in decreases (or increases) in the residence time and the number of bar impacts, but increases (or decreases) in the refining intensity and the specific refining power. The explanation is the same as that given for Figure 2.9. Decreases in the residence time r 2 and number of bar impacts -/V2 in the 2nd-stage (Figure 2.10) resulted in decreases in the total residence time rt — rx + r 2 and total number of bar impacts Nt — Nx + N2. As a result, the refining intensity It and specific refining power et over two stages increased. Increased total refining intensity (or specific refining power) brought up more fibre cut Chapter 2. Process Variables, Modelling and Simulations 50 and hence less long fibres and shives but more fines. As a result, the freeness, long fibre content and shive content decreased. Table 2.6 summarizes the simulation results given in Figure 2.10. For a given inlet consistency c, = 30.5%, a 28.6% increase in the rotational speed w2 (from 1400 rpm to 1800 rpm) decreased the residence time r 2 by about 53% (from 0.475 sec to 0.225 sec) and the number of impacts by about 36% (from 4250 to 2700). As a result, the refining intensity I2 increased by about 69.6% (from 5.6 x l 0 ~ 4 M J / k g to 9.5 x l O - 4 M J / k g . The increase in the specific refining power e2 was even more pronounced, which is about 105% from 5.5 MJ/kg /s to 11.3 MJ/kg/s . The effect of the inlet consistency is smaller compared to that of the rotational speed. Table 2.6 shows that at a given rotational speed, e.g. w = 1800 r/min, a 35% increase in the inlet consistency (from 26% to 35%) resulted in a 39% increase in the residence time, a 41% increase in the number of impacts, a 29% decrease in the refining intensity and a 28.7% decrease in the specific refining power, respectively. Table 2.7 shows that a 28.6% increase in the secondary rotational speed resulted in a 33% decrease in the total residence time rt, a 24% decrease in the total number of impacts Nt and a 28% increase in the total refining intensity It. In terms of pulp properties, a 28.6% increase in the secondary rotational speed decreased the 2nd stage freeness by 50%, the long fibre content by 18% and the shive content by 44%. At a given rotational speed w2 = 1800 r/min, Table 2.7 shows that a 35% increase in the inlet consistency (from 26 to 35%) on the 2nd stage brought about a 22% increase in rt, an 18% increase in Nt and 17.4% decrease in It. As a result, the freeness, the long fibre content and the shive content on the secondary stage increased by 21%, 11%, and 31% respectively. Chapter 2. Process Variables, Modelling and Simulations 51 4 1 ; i i 1 4 1 1 1 1 1 1400 1500 1600 1700 1800 1400 1500 1600 1700 1800 Figure 2.10: The 2nd stage refining at E2 = 2.5 GJ / t : effects of the rotational speed w and the inlet consistency Cj ( — : C J = 26%, -.-.: Cj = 30.5% and - -: q = 35%) Chapter 2. Process Variables, Modelling and Simulations 52 total residence time (s) 0.4 1400 1500 1600 1700 1800 10000 8000 6000 total no. of impacts 4000 1400 1500 1600 1700 1800 x 1 0-4total ref. intensity (MJ/kg) total specific ref. power (MJ/kg/s) 1400 1500 1600 1700 1800 1400 1500 1600 1700 1800 C S F (ml) 250 100 1400 1500 1600 1700 1800 LF (%) 1400 1500 1600 1700 1800 rotational speed (r/m) S C (%) 1400 1500 1600 1700 1800 rotational speed (r/m) Figure 2.11: Two-stage refining at E2 = 2.5 GJ / t : effects of the 2nd-stage rotational speed w and the inlet consistency Cj ( —: = 26%, - . - . :C j = 30.5%, and - - : C; = 35%) Chapter 2. Process Variables, Modelling and Simulations 53 Table 2.6: Effects of Secondary Stage Refining at E=2.5 G J / t w2 = 1400 - 1800 r/min Aw2 =29% = 30.5%) N2 I2 e 2 (sec) ( x l 0 ~ 4 MJ/kg) (MJ/kg/sec) 0.475 to 0.225 4250 to 2700 5.6 to 9.5 5.5 to 11.3 -53% -36% +69.6% +105% Ci = 26 - 35% Ad = 35% (w2 = 1800 r/min) N2 I2 e2 (sec) ( x l O - 4 MJ/kg) (MJ/kg/sec) 0.18 to 0.25 2200 to 3100 11.3 to 8.0 13.6 to 9.7 +39% +41% -29% -28.7% Table 2.7: Effects of Secondary Stage Refining on Pulp Properties, E=2.5 G J / t w2 = 1400 - 1800 r/min Aw2 =29% (Ci = 30.5%) Nt It CSF L F SC (sec) ( x l O " 4 MJ/kg) (ml) (%) (%) 0.75 to 0.5 7900 to 6000 7.8 to 10 200 to 100 45 to 37 0.27 to 0.15 -33% -24% +28% -50% -18% -44% Ci = 26 - 35 % Aci = 35% (w2 = 1800 r/min) n Nt It CSF L F SC (sec) ( x l O " 4 MJ/kg) (ml) (%) (%) 0.45 to 0.55 5500 to 6500 10.9 to 9.0 140 to 170 35 to 39 0.13 to 0.17 +22% +18% -17.4% +21% +11% +31% Chapter 2. Process Variables, Modelling and Simulations 54 2.5 Summary and Conclusions The purpose of this chapter was to obtain the quantitative relationships between pro-cess inputs and outputs. A set of equations was derived using the mechanistic method, through which major manipulating variables are linked to major output variables. As can be seen from process modelling, for given raw material and refiner design the pulp quality is mainly determined by specific energy and the rate of transferred energy, i.e., refining intensity (or specific refining power). Specific energy is a function of refiner motor load and throughput, while refining intensity is a function of specific energy, disc rotational speed and refining consistency. Although simulations showed that the effect of the rotational speed is the most apparent one, for a given specific energy refining con-sistency appeared to be the only readily available means of adjusting refining intensity in an industrial operation. Simulations also demonstrated the quantitative relationships between the variables and process nonlinearity. Chapter 3 Industrial TMP Process Identification 3.1 Introduction The model developed in the previous chapter provides a valuable tool for qualitative analysis and a better understanding of the process. For control purposes, process dy-namics and disturbances should be considered in the model development. However, it is not easy to model the dynamics of the process through mechanistic modelling and it is even impossible to characterize the disturbances by an analytical (or mechanistic) function [74]. A dynamic model which possesses maximum simplicity and a minimum number of parameters consonant with representational adequacy is usually achieved from process identification, i.e., a model is derived experimentally by using input-output data collected from a real process (Figure 3.12). One of the advantages in using process identification or the empirical method is that only quantitative knowledge of the process mechanism is required to find rough estimates of process input-output relationships. Secondly, using process identification, dynamic modelling can be realized more easily in less time than using the mechanistic analysis where the development of complex differential equations may be required. Lastly, due to the impossibility of modelling disturbances by the mechanistic method, process identification becomes a successful tool in characterizing disturbances. Conventional identification methods are based on the choice of special inputs to a process, such as step, pulse and sinusoidal inputs. These methods have been useful 55 Chapter 3. Industrial TMP Process Identification 56 v(t) u(t) Dynamic y(\ Process Figure 3.12: Input to and output from a dynamic process when a process is affected by a small disturbance. However, they have not always been successful [75] and usually provide poor results on processes characterized by long time constants or excessive measurement noises. This is because large process input manipu-lations are needed to generate the required information at process outputs, thus risking sustained deviations from product quality targets and unsafe operating conditions. The approach adopted in this study is time series analysis which refers to the analysis of data from random observations with time (Figure 3.12). The reasons for using time series techniques are capabilities of (i) identifying accurate models even when the input vari-ations are small, (ii) modelling the effects of all unmeasured external disturbances and measurement noises, and (iii) providing a general methodology for identification of the multivariate processes. The objective of industrial T M P process identification is to find a dynamic-stochastic model suitable for the control of a T M P plant. To achieve this goal, a number of mill trials have been performed at a B .C. pulp and paper mill. A large number of on-line operational data have been collected. With the help of identification techniques and time series analysis, intensive data analysis has been performed to interpret the data. Section 3.2 describes the choice of a model structure for process identification. Section 3.3 details the time series analysis techniques used for process identification. Section 3.4 presents the procedures for identifying, fitting, and checking models when simultaneous Chapter 3. Industrial TMP Process Identification 57 pairs of observations («i, yi), (u2, y2), • • •, (u>t, Ut) of input and output are available. Sec-tion 3.5 presents identification results and discussions of the results. The summary and conclusions are given in Section 3.6. 3.2 Process Model Before performing process identification, it is always necessary to assume the type of model according to the prior knowledge of the process and the purpose of the model development. For control purposes, we assume that a stochastic process (simply referred to a process) can be described by a discrete-time dynamic-disturbance model as follows for a SISO case: Vt = Gp(z-l)uM + Nt (3.41) where yt is the output deviation from its setpoint at time t and ut is the input deviation from its corresponding steady-state value. The parameter k(> 0) is the number of sam-pling intervals of pure process delay. One additional period of delay is introduced by the sample and hold operation of the computer. Gp{z~x) presents the input-output dynamics of the process, such as the process gain and time constant. Nt characterizes all unmea-sured external disturbances and measurement noise as observed at the output. Box and Jenkins [75] suggest that the transfer function-stochastic model above be parameterized as: Biz-1) Ciz-1) ,0 As can be seen, the disturbance term Nt in Equation 3.41 is modelled by passing an independent, normally distributed (white) noise sequence et with variance cr2 through a linear dynamic transfer function D C ^ - i ) ^ d (A = 1 - z~x). The advantages of using the Box-Jenkins model are that first, the allowance of d roots equal to unity (d = 0 or 1 usu-ally) enables one to characterize the type of nonstationary drifting behavior that process Chapter 3. Industrial TMP Process Identification 58 variables tend to exhibit when uncontrolled; secondly, because an integrator function is modelled in the disturbance model, it naturally leads to integral action in a model-based controller design. Finally, by using the Box-Jenkins model, the effects of both manipu-lated input and disturbance can be modelled separately. Terms A(z~1), B{z~l), C(z~1) and D(z~1) are the polynomials in the backward shift operator z~x (i.e., z"xyt = y t-i) and they are described as follows: A(z~x) = 1 + a i Z ' 1 + a2z~2 + ••• + anaz-na Biz'1) = b0 + hz-1 + b2z-2 + ••• + bnbz-nb (3.43) C{z~x) = 1 + cxz~x + c2z~2 + ••• + cncz~nc D(z~l) = 1 + dyz~x + d2z~2 + ••• + dndz —nd 3.3 Analysis Tools Once the data has been collected, a number of statistical analyses are performed, con-sisting of (i) the auto-correlations, (ii) the partial auto-correlations, and (iii) the cross-correlation functions. These are aimed primarily at identifying the model structure (Equation 3.42) from input-output data, and then performing model parameter estima-tion and model checking. The calculation of the auto-, partial auto- and cross-correlations are given in Appendix A , and associated properties can be found in the book [75] by Box and Jenkins. 3.4 Identification Procedures Once the data has been collected, the steps involved in model identification are (i) iden-tification of the model structure (na,nb,nc,nd,k), (ii) model parameter estimation, and (iii) identified model checking. The following presents the details of each step in process model identification. Chapter 3. Industrial TMP Process Identification 59 3.4.1 Model Structure Identification Differencing Data to Obtain Stationarity The time series analysis tools given in Section 3.3 are based on the assumption of the stationarity of a process, which implies that the process has constant mean fj, and con-stant variance a2. In practice, however, this is not always true. The nonstationarity of a process exists when time series auto- and cross-correlations do not die away, but rather are composed of many significant correlations even at high lags. Alternatively, the nonstationary behavior is suspected when a time series does not vary about a fixed mean value, but rather slowly meanders from one operation point to another. The stationarity of a process can be achieved by differencing the process data (Ut, Yt) pair d times (d is usually 1 or 2) until the stationary data (u t, yt) pair is obtained, where ut = AdUt and yt = A < % Identifying the Structure of a Transfer Function The structure (na, nb, k) (Equation 3.43) of a transfer function is obtained by estimating the impulse responses of the process and then comparing them with theoretical ones. Theoretical examples of impulse responses for some common processes can be found in [75]. In the case of open-loop process identification, if the input series ut is white (uncorre-cted), then the process impulse response is proportional to the cross-correlation ruy(k) and can be obtained by multiplying the cross-correlation coefficients ruy(k) by ay/au. However, if the input signal ut is not white, prewhitening the input has to be done before the estimation of impulse response. Prewhitening can be obtained by fitting the input ut to a time series model, i.e., u't = 0 ( ; z _ 1 )0 - 1 ( ; z - 1 )M T and transforming the correlated input series ut to the uncorrelated white series u't. The same transformation is applied to the output series yt to obtain y't. After prewhitening, the cross-correlation function between Chapter 3. Industrial TMP Process Identification 60 the input u't and correspondingly transformed output y't is proportional to the impulse response. Alternatively, the impulse response can be obtained without prewhitening by relating the input-output cross-covariance to the input auto-covariance. Suppose that a process can be adequately approximated by the following transfer function: yt = v(0)ut + t ; ( l )« t _i + v(2)ut-2 H = [v(0) + v(l)z-1 + v(2)z-2 + ---]ut (3.44) where the weights v(0),v(l),v{2) • • • in the above equation are called the impulse re-sponse. Multiplying the above equation by ut-k for k > 0 and taking expectations gives: cuy(k) = v(0)cuu(k) + v(l)cuu(k - 1) + k = 0,1,2,- •• (3.45) Supposing that the weights v(i) are effectively zero beyond k = K, then the above equation with k = 1, 2,3 • • • K can be written in a more compact form as: C U y ^ U U (3.46) or where cuy(0) c„ u(0) Cuu(1) Cuu(K) «(0) C u y Cuu(1) CUu(0) • CuU(K - 1) V = v(l) CUy(K) Cuu(K) cuu(K - 1) • Cuu(0) v(K) (3.47) Identifying the Structure of a Disturbance Model Chapter 3. Industrial TMP Process Identification 61 The idea of identifying the structure of a noise model ( Equation 3.42) is the same as the idea of identifying the structure of an input-output transfer function. Once an acceptable input-output transfer function is obtained (see Section 3.4.3), the structure (nc,nd) of the noise model can be obtained through calculating the residual auto- and partial auto-correlations and comparing them to the theoretical ones of known processes. Theoretical examples of the auto- and partial auto-correlations of some known processes can be found in [75]. 3.4.2 Parameter Estimation Once an acceptable structure (na, nb, nc, nd, k) of a dynamic-noise model is obtained, es-timation of the parameters of A(z~x), -B(z _ 1 ) , C(z~x) or D(z~x) polynomials in Equation 3.42 is carried out by minimizing the sum of prediction error as follows: m i n ^ e 2 (3.48) t=i where where e t is calculated at given (yt,ut-k) data pair and a set of initial parameters of A(z-x),B(z-x),C(z-x) and D(z~x) polynomials. 3.4.3 Model Checking The tools used for checking the accuracy of a model are the residual auto-, partial auto-and cross-correlations. A n adequate input-output model is obtained when the residual cross-correlations die away within a defined confidence interval. Similarly, a disturbance model is adequate when the residual auto- and partial auto-correlations die away within the confidence interval. Chapter 3. Industrial TMP Process Identification 62 3.5 Industrial TMP Process Identification The T M P process identification aimed at identifying the dynamics between inputs and outputs, including model order, time-delay, time constant and process gain. A multi-variable relationship between inputs and outputs was obtained through identifying the relationship for each input-output pair by doing experiments separately. The following presents the details of the mill trials, on-line data and model identification results. 3.5.1 Process Description and Experiment Design A number of mill trials were conducted on the primary refiner of a two-stage C T M P production line at a B . C . pulp and paper mill. Heated chips were fed from a vertical steaming vessel to a stream splitter conveyor by a pressurized plug screw discharger. From the stream splitter conveyor, the chips were metered by two load sensing conveyors into the drive and tail ends of the pressurized primary refiner (a twin-disc refiner). Dilution water was sent into the drive and tail ends of the primary refiner to obtain the desired consistency. From the primary refiner, pulp was blown into the pressurized primary cyclone and was then fed to the pressurized secondary refiner. Because of pulp quality specifications, the process was operated under the conditions of minimizing the deviations of the process variables from their targets. This made the identification more difficult because of the lack of excitation for the process identification. Therefore, Pseudo Random Binary Sequences (PRBS) - a series of uncorrelated input changes, were added as input perturbation signals during identification. The advantage of using P R B S signals rather than a classical bump is that PRBS signals of much lower amplitude can be used, thereby minimizing the effects of the test on the process. Fur-thermore, because a PRBS signal is continually switching states, the process is excited during the entire time history rather than just at the start of the record as is the case Chapter 3. Industrial TMP Process Identification 63 during traditional bump tests. The amplitude of a PRBS perturbation signal is an important design variables and may be designed by placing hard limits on the input variations, i.e., < ut < uu in order to minimize process deviations from the desired targets while generating as much information on the process dynamics as possible. The upper and lower limits should be chosen based on how much variation is tolerable. The other design variable is the switching period T of a PRBS signal. In order to have an input with a proper frequency content, the switching period T should, in general, be some integer multiple of the sampling period T3, i.e., T = n x Ts [76]. For the duration of our mill trials, the amplitudes of PRBS signals were limited to ±0.8% of each input and the switching period T was chosen as two times of the sampling rate, i.e., T = 2 Ta. Key operating data was collected from Bailey Distributed Control Systems (DCS), where the sampling rate was one second, i.e., Ts = 1 sec. 3.5.2 Identification Results The initial plan of process identification was to identify the dynamic relationship between inputs (closing pressure, chip screw speed and dilution water) and outputs (motor load and consistency). The simplest way for doing so is to identify the relationship between each input and output separately by doing the experiments separately, i.e., perturbing one input at a time while keeping other inputs constant. Closing Pressure (CP) to Motor Load (ML) Figure 3.13 shows the refiner motor load (ML) response to the changes in the hydraulic closing pressure (CP). Figure 3.14 shows the unit impulse response and its 99% confidence level, estimated using the motor load/closing pressure data from Figure 3.13. In practice, increasing C P increases M L . This has been indicated by the first few weights of the impulse responses. The second, third and forth impulse weights are significant, indicating Chapter 3. Industrial TMP Process Identification 64 Table 3.8: M L / C P Parameter Estimates, Estimated Standard Deviations and F P E Model d C(z-X) D{z-X) k F P E (Eqn.) std std std std 3.50 0 1 -0.8678 0.0061 1 0.3941 1 -0.5096 1 0.132 0 0.1159 0.0026 0 0.0867 0 0.0831 3.51 1 1 -0.8852 0.0060 1 0.1473 1 1 0.1719 0 0.3141 0.0030 0 0.0711 0 that the process may have one or more complete sampling intervals of pure time delay (k> 1). The remaining impulse weights suggest one of several alternatives for the orders of the transfer function polynomials and B(z~x). The two most likely alternatives are first order or second order. Note that the impulse response (Figure 3.14) also shows a periodic variation of 5 sec in the motor load. This variation might result from the periodic variation in the chip feed. A n iterative model-building procedure consisting of structure identification, fitting, and model checking stages identified the following dynamic-stochastic model: 0.0061 1 + 0.3941Z"1 ,„ . * = l - O ^ f e - i " * - 3 + l - 0 . 5 0 9 6 ^ e t ( 3 - 5 0 ) Table 3.8 shows the parameter estimates, their estimated standard deviations and Akaike's F P E (Final Prediction Error) [77]. The above equation was further tested on a fresh data set that was not used in the model estimation. Associated residual tests in Figure 3.15 in-dicate that the above dynamic-stochastic model is adequate because the cross-correlations all fall in the 99% confidence intervals and the residuals are nearly white. The input-output transfer function part, however, may be of limited use if it was used to design a minimum variance controller or a generalized predictive controller. The controller would not contain integral action because of the lack of nonstationarity in the disturbance model. If this were the intended use of the model then it would be necessary Chapter 3. Industrial TMP Process Identification 65 to force an integrator into the disturbance model even though this may not be entirely supported by the data. Fitting the noise model with an integrator gave the following parameter estimates: 0.006 1 + 0 .1473Z- 1 Vt = --, n o c r o r - l " t - 2 + £ e< (3.51) 1 - 0.8852z-Table 3.8 shows the parameter estimates, estimated standard deviations and Akaike's F P E . The residual correlations calculated using a fresh data set are given in Figure 3.16, indicating that the above model is adequate. Figure 3.17 shows the fresh data set and the model prediction using the same data set. 200 250 150 200 250 Time (sec.) 400 Figure 3.13: Motor Load (ML) response to Closing Pressure (CP) Chapter 3. Industrial TMP Process Identification 66 Impulse response estimate c ' 0 ) CD ) ? o ? ! . CP I I ? 4 I4 H 1 i A i i i i i 11 1 1 . . 1 0 5 10 15 20 25 lags Figure 3.14: M L / C P impulse response estimate Chapter 3. Industrial TMP Process Identification 67 Correlation function of residuals. Output # 1 5 10 15 20 lag C r o s s corn function between input 1 a n d residuals from output 1 - 2 5 - 2 0 - 1 5 - 1 0 25 Figure 3.16: M L / C P residual correlations (d=l) 1 5 r 14 -s 13 * 12 -11 -~ i V i r -_ J L_ 20 40 60 80 100 120 140 160 180 200 Figure 3,17: C P / M L input, output and model prediction Chapter 3. Industrial TMP Process Identification 68 Table 3.9: M L / T S Parameter Estimates, Estimated Standard Deviations and F P E Model d A(z-X) Biz'1) Ciz-1) D(z-1) k F P E (Eqn.) std std std std 3.52 0 1 -0.8186 0.0616 1 0.5900 1 -0.6159 8 0.1148 0 0.2527 0.0563 0 0.0643 0 0.0665 3.53 1 1 -0.8997 0.0548 1 0.5768 1 8 0.1425 0 0.4630 0.0672 0 0.0561 0 Transfer Screw Speed (TS) to Motor Load (ML) Figure 3.18 shows the motor load (ML) response to the changes in the transfer screw speed (TS). To have an initial guess of the model structure, the unit impulse responses are estimated and given in Figure 3.19, using the raw data in Figure 3.18. In practice, increasing the chip feed speed increases the motor load. The first eight impulse weights in Figure 3.19 are insignificant, indicating that the process may have eight sampling period intervals of pure time delay (k = 8). This delay is attributed to the finite time required to deliver wood chips to the inlet of the primary refiner. The remaining impulse weights suggest one of several alternatives for the orders of the transfer function polynomials. The most likely alternatives are first order with fractional delay, or second order. A n iterative modelling procedure consisting of structure identification, model fitting, and model checking stages identified the following dynamic-stochastic model structure: 0.0616 1 + 0 .59Z- 1 , 0 V t = 1 - 0 . 8 1 8 6 , - ^ - 9 + 1-0.6159*-!* ( 3 ' 5 2 ) Table 3.9 shows the parameter estimates, their estimated standard deviations and Akaike's F P E . A check on the above model was made by calculating the residual auto- and cross-correlations with the use of a fresh data set that was not used for the model estima-tion. The residual correlations in Figure 3.20 indicate that the identified the dynamic-stochastic model (Equation 3.52) is adequate. Chapter 3. Industrial TMP Process Identification 69 To design a controller containing integral action, we forced an integrator into the disturbance model and obtained: Vt = 0.0548 1 + 0.57682- 1 - w t _ 9 H et 1 - 0.89972" A (3.53) Table 3.9 shows the parameter estimates, estimated standard deviations and Akaike's F P E . A check on the above model was made by calculating the residual auto- and cross-correlations with the use of a fresh data set. The residual correlations in Figure 3.21 indicate that the identified dynamic-stochastic model (Equation 3.53) is adequate. Figure 3.22 shows the fresh data set and the model prediction. 150 200 250 Time(sec) 400 Figure 3.18: Motor Load (ML) response to Transfer Screw Speed (TS) Chapter 3. Industrial TMP Process Identification 70 Impulse response estimate c > < > c p ( j. ' > <j ) a a ? T T 1 1.21 ' ' ' ' 1 O 5 10 15 20 25 lags Figure 3.19: M L / T S impulse response estimate Correlation function of residuals. Output # 1 ) 21 i i 1 .—i 1 1 1 1 1 1 - 2 5 - 2 0 - 1 5 - 1 0 - 5 0 5 10 15 20 25 lag Figure 3.20: M L / T S residual correlations Chapter 3. Industrial TMP Process Identification 71 Chapter 3. Industrial TMP Process Identification 72 Dilution Water (DW) to Motor Load (ML) Because the primary refiner under this trial is a twin-disc refiner, two sets of PRBS signals were added to perturb the valve positions on both the drive and tail ends in order to vary the total dilution flow (Figure 3.23). To find the relationship between the motor load (ML) and total (total = drive end + tail end) dilution water (DW), M L response to the changes in the total D W is plotted out in Figure 3.24. Figure 3.25 shows the unit impulse weights estimated using the motor load/ dilution flowrate data from Figure 3.24. A n increase in D W will decrease M L . This has been indicated by the first few negative impulse weights. The first and second impulse weights in Figure 3.25 are insignificant, indicating the process may have two sampling periods of pure time delay (A; = 2). The remaining impulse weights in Figure 3.25 suggest one of several alternatives for the orders of the transfer function polynomials. A n iterative model-building procedure consisting of structure identification, fitting, and model checking stages identified the following dynamic-stochastic model: -0.0089 1 + 0.4878Z-1 ' , n r A . V t = l - 0 . 2 4 5 , - ^ - 3 + 1 - 0.545,-1 * ( 3 " 5 4 ) Table 3.10 shows the parameter estimates, estimated standard deviations and Akaike's F P E . The residual correlations in Figure 3.20 indicate that the above identified model is Table 3.10: M L / D W Parameter Estimates, Estimated Standard Deviations and F P E Model d Biz'1) Ciz~l) Diz-1) k F P E (Eqn.) std std std std 3.54 0 1 -0.245 -0.0089 1 0.4878 1 -0.545 2 0.1155 0 0.201 0.002 0 0.0614 0 0.0607 3.55 1 1 -0.2653 -0.0087 1 0.3765 1 2 0.1468 0 0.2178 0.0020 0 0.0535 0 adequate. Chapter 3. Industrial TMP Process Identification 73 Forcing an integrator into the disturbance model gave: Vt -0.0087 1 + 0.3765Z-1 -«t_3 + et 0.2653z- A (3.55) Table 3.10 shows the parameter estimates, estimated standard deviations and Akaike's F P E . The residual correlations in Figure 3.20 indicate that the above identified model is also adequate. 150 200 T i m e (sec.) 350 350 350 350 350 Figure 3.23: Changes in Valve Positions (VP), Dilution Water (DW) and Motor Load (ML) Chapter 3. Industrial TMP Process Identification 74 15 2 5 01 , , , , , j 0 50 100 150 200 250 300 Time (sec.) Figure 3.24: Motor Load (ML) response to total Dilution Water (DW) Impulse response estimate ! j , , CP T T ? ! ? I * A l i > 0 I , , , 1 1 0 5 10 15 20 25 lags Figure 3.25: M L / D W impulse response estimates Chapter 3. Industrial TMP Process Identification 75 Correlation function of residuals. Output # 1 1t , , • i — lag Figure 3.26: M L / D W residual auto- and cross-correlations Chapter 3. Industrial TMP Process Identification 76 Table 3.11: Identification Results Model Input Output Gain Time Constant Time Delay (eqn.) (sec.) (sec.) 3.50 C P M L 0.0461 (MW/psi) 7.052 1 3.52 TS M L 0.3396(MW/rpm) 4.995 8 3.54 D W M L -0.0118(MW/L/min) 0.711 2 3.5.3 Result Discussions As mentioned earlier, the simplest way to find rough estimates of a multivariable relation-ship between inputs and outputs is to identify each input-output relationship separately. Identification results showed that wood chip refining is a dynamic-stochastic process with lower order dynamics plus time delay. To have some physical knowledge of the process, the identified discrete dynamics of the process were transformed into the continuous ones. Associated process gains, time delays and time constants are given in Table 3.11. A higher static gain between the transfer screw speed and the motor load indicates that the motor load is strongly affected by the chip feed. By comparison, the dilution water seems to have less impact on the process. Shorter time constant between the dilution flow and the motor load was identified, while longer time constants between the closing pressure, screw speed and the motor load were observed. A longer time delay between the screw speed and the motor load was identified. This was attributed to transferring wood chips to the inlet of the primary refiner. Because of the lack of an on-line consistency sensor during the mill trials, no identification was performed for the relationship between the consistency and other variables. Chapter 3. Industrial TMP Process IdentiEcation 77 d i s t u r b a n c e s + C P T S 1 st o r d e r •+- de lay N . L . M e c h . m o d e l M L D W C o n s C P : hydraul ic c l o s i n g p r e s s u r e T S : c h i p transfer s c r e w s p e e d DW:di lu t ion water flowrate M L : motor load C o n s : c o n s i s t e n c y Figure 3.28: Dynamic and static relationships between inputs and outputs 3.6 Summary and Conclusions For control purposes, this chapter has tried to find dynamic relationships between input and output variables of a T M P plant. To avoid the complexity of developing differential equations in mechanistic modelling, process identification (an experimental method) was used to find rough estimates of the dynamics and disturbances of the process. Identi-fication results indicated that wood chip refining is a dynamic-stochastic process with excessive stochastic disturbances. Due to product quality concerned, the input pertur-bation signal was limited to a certain rage therefore signal to noise ratio was limited. For this particular case, the model was limited to lower order such as second order in order to obtain rough estimates of process dynamics with representational adequacy. It was found that first-order dynamics plus time delay gave more adequate approximation in terms of minimizing Akaike's F P E and residual cross-correlations. As mentioned earlier, the simplest way to find estimates of a multivariable relationship between inputs and outputs is to identify the relationship for each input-output pair separately. Figure 3.28 shows identified multivariable dynamics of the process. As can be seen from mechanistic modelling, the process has nonlinear and time-varying characteristics. Therefore, the process may be described by its dynamics with unit gain plus the nonlinear relationship obtained from mechanistic modelling stage. Chapter 4 Constrained Multivariable Model-Based Predictive Control 4.1 Introduction Model-Based Predictive Control (MBPC) refers to a class of algorithm that computes a sequence of manipulated variable adjustments in order to optimize the future behav-ior of a process. Originally developed to meet the specialized control needs of power plants and petroleum refineries, M B P C technology is probably the most important ap-proach to the advanced control of complex interacting industrial processes [78] and is gaining widespread acceptance in various industries. M B P C ' s effectiveness in practice was discovered in the early 80's by the development and application of model predictive heuristic approaches: I D C O M (Identification/Command) by Richalet et al. [79] and D M C (Dynamic Matrix Control) by Cutler and Ramaker [80]. In both algorithms an explicit dynamic model of a plant is used to predict the effect of future actions of the manipu-lated variables on the output (thus the name "Model Predictive Control"). The future control actions are determined by minimizing the predicted error subject to operating constraints. The calculation of future control actions is repeated at each sampling time based on updated measurements from the plant. A recent paper written by Froisy [81] provides an interesting industrial perspective of current M B P C technology and summa-rizes likely future developments. Clarke and Zhang [82] compared the different M B P C approaches and found them all suited to plants that exhibit variable dead-time and non-minimum phase characteristics. A comprehensive survey of the history and theory of 78 Chapter 4. Constrained Multivariable Model-Based Predictive Control 79 model-based controls can be found in various papers [83, 84, 85, 86, 87]. Generalized Predictive Control (GPC) first introduced by Clarke et. al. [88] is one of a class of M B P C algorithms. Few advanced control methods have had as much influence, widespread acceptance and success in industrial applications as the G P C approach. The success of this technique is due to its capabilities of controlling a process with: • variable time delay and model order; • over-parameterization (plant/model mis-match); • unstable zeros (non-minimum phase); • unstable poles; • load-disturbances. In addition, G P C can successfully overcome difficulties encountered by other control methods. Minimum-variance self-tuning has been used for a long time, but it is highly sensitive to the assumption made about the dead-time. Pole-placement and Linear Quadratic Gaussian (LQG) self-tuners perform badly if the order of a process is overes-timated because of pole/zero cancellations in the identified model. Although self-tuning and adaptive controls have made much progress over the previous decade in terms of theoretical understanding and practical applications, neither one is suitable as a general purpose algorithm for the stable control of the majority of real processes [88]. The G P C algorithm possesses simplicity and flexibility, appearing to overcome the problems in one algorithm. As far as process constraints are concerned, G P C has the capability of incor-porating and handling the constraints explicitly. However, relatively little information is available on constraint handling using other methods. The refiner mechanical pulping process under this study is a rather complex process. Good control will improve the production, efficiency and pulp quality, while poor control Chapter 4. Constrained Multivariable Model-Based Predictive Control 80 can lead to operating difficulties and poor physical and papermaking properties. Devel-oping a suitable control strategy for the process is a very challenging task, because of the multivariable nature, nonlinearity, process interactions, time varying, time-delay, vari-able constraints and stochastic disturbances of the process. A controller for the process must be able to handle the problems together with the robustness to modelling errors. Constrained multivariable adaptive-predictive control approach presented in this chapter will be used for the control of the process. The controller is developed based on the G P C approach, where future control actions are determined by minimizing the predicted error subject to input-output constraints. A process model is used to predict future process output under the influence of future control actions. The model parameters are estimated at each sampling time based on update measurements from the process. Quadratic programming and mixed-weights least-squares algorithms are introduced to solve the constrained control problem. The proposed controller will possess a number of important features as follows: • The controller is able to handle process/model mismatch (inaccurate time-delay and model order) and disturbances. • The controller adopts an integrator as a natural consequence of its assumption about the basic process model to obtain offset-free close-loop behavior of the con-trol. • A n embedded adaptive scheme can help to handle process nonlinear and time-varying characteristics under the help of a proper parameter estimation method. • A n explicit model is used in the controller development so that fewer parameters are needed to be estimated on-line than equivalent implicit schemes and thus it could be very beneficial to the parameter estimation of an M I M O model. Chapter 4. Constrained Multivariable Model-Based Predictive Control 81 • The proposed predictive control strategy can help to handle variable constraints explicitly and naturally by taking in account the constraints into the controller development. This chapter is structured as follows. Section 4.2 gives a brief introduction to uncon-strained Generalized Predictive Control (UCGPC) for both single-input single-output (SISO) and multi-input multi-output (MIMO) cases. Section 4.3 describes recursive parameter estimation methods that may incorporate with an adaptive control scheme to compensate for the time-varying and nonlinear features of the process. Section 4.4 presents a predictive control problem subject to input and output constraints. In Sec-tion 4.5, an analytical solution of the constrained G P C is derived by solving analytical quadratic programming (QP) problem. Section 4.6 presents the derivation of the con-strained G P C for the general cases by solving a mixed-weights least-squares (MWLS) problem. The associated weights in M W L S can be modified in a manner that encom-passes both the requirement for the future control moves to lie inside the feasible region and minimization of desired control objective function. Simulations of the proposed con-trol strategy are presented in Section 4.6. Summary and conclusions are given in Section 4.7. 4.2 Unconstrained Generalized Predictive Control 4.2.1 Basic SISO Generalized Predictive Control Model and Output Prediction A model is the center for any kind of model-based control designs. The model used in G P C design is the Controlled Autoregressive and Integrated Moving Average (CARIMA) Chapter 4. Constrained Multivariable Model-Based Predictive Control 82 model as: Aiz-^yit) = B(z-l)u(t -k) + (4.56) where y(t) is the process output, u(t) is the input (y and u are deviation variables), and is uncorrelated, normally distributed zero-mean random variable (or white noise) with the variance <r|. A: is the time delay. The difference operator A = 1 — z _ 1 in the denominator of the noise term is widely assumed, as it forces an integrator into the controller in order to eliminate offset between the measured output and its setpoint. A(z_1), B(z~l) and C(z~l) are the polynomials in the backward-shift operator z"1 with the orders of na, nb and nc respectively: A(z-X) = 1 + axz'x + a2z~2 + ••• + anaz-na Biz'1) = b0 + blz-1 + b2z-2 + --- + bnbz-nb (4.57) Ciz'1) = 1 + cxz" 1 + c2z-2 + • • • + cncz-nc The above C A R I M A model has a more general application. It has also been used by Tuff and Clarke to derive the generalized minimum variance control (GMV) and pole-placement self-tuners with inherent integral action [89]. To derive a j-step ahead predictor y(t+j\i) from Equation 4.56, consider the following Diophantine identity: Ciz-1) = Ejiz'^Aiz'1) A -rz-iFjiz-1) (4.58) Ejiz"1) and Fjiz'1) (order j — 1 and max(na, nc — j), respectively) are determined by the given polynomials A ( z - 1 ) , C ( z _ 1 ) and the prediction interval j. Note that the disturbance in Equation 4.56 is modelled as a nonstationary process by A ( z - 1 ) . If the model is obtained using identification techniques, C ( z _ 1 ) is the polynomial affected by a bigger error. Therefore, in G P C , the Ciz'1) polynomial is considered as a design polynomial chosen by the user to obtain good behavior of the controller [88]. Chapter 4. Constrained Multivariable Model-Based Predictive Control 83 Multiplying Equation 4.56 by Ej(z *) A z* and combining with Equation 4.58 gives: y(t + j) = Gjiz'1) A u'(t + j - l ) + Fjiz-Wit) + Ejiz'1)^ + j) (4.59) where Gj(z~x) — Ej[z~1)B(z~x') in the order of j + nb— 1. Note, A; — 1 leading zeros have been added to B(z~l) in Equation 4.56 to simplify the notation. Hence the polynomial B(z~x) is in the order of nb — nb+ k — 1. The superscript ^ in the above equation denotes filtering by 1 /C(z _ 1 ) . As the polynomial Ej(z~x) is of the degree j — l(j > 1), the noise terms Ej(z~x)^(t + j) in the above equation are all future unknown values and thereby the output prediction in Equation 4.59 can be modified as follows: y(t + j\t) = G^z-1) A U*(t + j - l ) + V ( * ) (4.60) To separate past known filtered control actions from current control actions yet to be determined, consider the following identity: Gjiz'1) = G'jiz-^Ciz'1) + z-JTjiz-1) (4.61) where G'^z'1) and Yj(z~x) (order j — 1 and max(nb — l ,nc — 1), respectively) are defined by Gj^z'1), C ( z _ 1 ) and j. Substituting the above equation into Equation 4.60 gives: y(t + j\t) = G ' ^ A u i t + j-V + Tjiz-^Aufit-V + Fjiz-Wit) (4.62) = G'jiz-1) Au(t + j - l ) + f(t + j) where f(t + j) = Tjiz-1) A - 1) + V W (4-63) As can be seen from the above equation, the output prediction y(t + j\t) consists of two types of process responses: • G'^z'1) A u(t + j — 1) is called "forced response" of the process. It depends on the future control actions yet to be determined. Chapter 4. Constrained Multivariable Model-Based Predictive Control 84 • fit + j) is called "free response" of the process. It depends on the past filtered known controls together with filtered measured outputs under the assumption that no future control actions are performed. Predictive Control Law The control law of G P C for a SISO case is derived from a quadratic optimal cost function or the function based on a designed trade-off between the control performance and soft constraints on the control actions as follows: { N2 NU ~\ E [Vi* + J) ~ w(t + J)f + E PU) A u\t + 3-1)) (4.64) where y(t + j) : the future output w(t + j) : the future setpoint u(t + j-l) : the future control action JVx : the minimum predictive horizon, (JVi > 1) N2 : the maximum prediction horizon, (iV2 > iV l ) NU : the control horizon, (N2 > NU > 1) Pti) : the control-weighting factor The objective of the predictive control is to drive future plant output y(t+j) close to the future setpoint w(t+j) in some sense by selecting a control movement so that the desired performance index of Equation 4.64 is minimized. The future plant output y(t + j) is predicted based on the plant model with knowing the current and past control actions and past plant output. To simplify the following derivation, the control weighting factor and the minimum predictive horizon are set to p(j) = p and Ni — 1, respectively. The objective function Chapter 4. Constrained Multivariable Model-Based Predictive Control 85 (Equation 4.64) and the output prediction (Equation 4.62) can then be stated in more compact forms as: JGPC = (y - w ) T ( y - w) + p u T u (4.65) and where y = G ' u + f y = [y(t + l),y(t + 2),---,y(t + N2)}T u = [Au(t),Au(t + l),---,Au(t + NU-l)}T f = [f(t + l),f(t + 2),...,f(t + N2)}T w = [w(t + l),w(t + 2),---,w(t + N2)]T and the matrix G ' is lower-triangular of dimension N2 x NU: G ' = 9o 0 9i 9o 9NU-I 0 0 0i 9o 9N2-I 9N2-2 ' ' ' 9N2-NU (4.66) (4.67) Substituting Equation 4.66 into 4.65 and minimizing with respect to ii yields: u = ( G ' T G ' + p I ) - 1 G ' T ( w - f) (4.68) As G P C is a receding-horizon control strategy, only the first control increment Au(i) in the control vector u is implemented into the plant under the assumption that beyond the control horizon NU, further increments in control actions are zero, i.e., Au(t + NU) = Au(t+NU+l) =,•••,= Au{t+N2) = 0 (or u(t+NU-l) = u(t+NU) = u(t+NU+l) = Chapter 4. Constrained Multivariable Model-Based Predictive Control 86 ,•••,= u(t + N2)). At the time t + 1, the computation is repeated with the horizon moved by one time interval. The G P C control algorithm described above represents an extremely flexible approach for the control of different processes due to its various design and tuning parameters incorporated into the algorithm. The role of setting these design and tuning parameters is to obtain a satisfactory closed-loop performance for different processes. Predictive horizon (Ni,N2): predicting the output response over the predictive hori-zon enables compensation for any variable or unknown time delays. The minimum pre-dictive horizon Ni can often be taken as 1 if the plant dead-time is unknown. A small value of the maximum predictive horizon N2 will result in more active control effort and the performance is quite similar to the performance of a minimum-variance control. While with a larger value of N2 -» oo, the control performance will be very sluggish. Control horizon (NU): a decrease in the control horizon NU will significantly reduce the computational effort and simplify control derivation, while an increase in NU will result in excessive movements in the control action and plant output. Control weighting (p): it is designed to reduce the control movements. Zero control weighting (i.e. p = 0) gives a dead-beat response, while with a very large value of p, the control action is dampened significantly. To reduce computation complexity, it is always advisable to take prior knowledge of plant delay and order into account so that Ni and N2 can be selected. It is also possible to reduce the computational burden by imposing a constant control input after some control horizon NU. It is found that a large class of plant models can be stabilized by G P C with Ni = 1, N2 = 10 [88] and that it will give an acceptable control by setting NU = 1 for simple processes. The details of guidelines for the principal design parameters and stabilization theorems are available in [90, 91, 92]. Chapter 4. Constrained Multivariable Model-Based Predictive Control 87 4.2.2 Basic MIMO Generalized Predictive Control In the following, the development of an M I M O G P C is a direct extension of the SISO G P C described in the previous section. The derivation of a multivariable G P C can also be found in the papers by Shah et al.[93] and Mohtadi et al.[94]. The following is a brief description of the derivation of the MIMO G P C . Model and Output Prediction Consider an M I M O C A R I M A model for a m-input n-output process as follows: A{z-1) A y(t) = Biz-1) A u(t - 1) + C{z~x)e{t) (4.69) where • u(£),y(£) and e(t) are the vectors of the process inputs, outputs and disturbances of the dimension m, n and n respectively. • A and C are diagonal polynomial matrices of the dimension n x n, while B is a polynomial matrix of the dimension n x m. The time delay of the process has been included in B matrix. If any of the input-output pairs of the process has a non-zero dead-time, the leading elements of the corresponding polynomial in B are zeros. To derive a j-step ahead predictor f(t + j\t), the following matrix Diophantine equa-tion is considered: C ( z _ 1 ) = Ejiz-^Aiz-1) A +z-jTj(z~1) (4.70) where £j(z~x) and ^ ( z _ 1 ) are diagonal polynomial matrices, defined by A(z-X) and C ( z - 1 ) . £ j ( 2 - 1 ) is has a degree of j - 1. Multiplying Equation 4.69 by £j(z-1)zj and Chapter 4. Constrained Multivariable Model-Based Predictive Control 88 combining with Equation 4.70 yields: y(t + 3) = GAz-1) A uf(t + j - 1) + ^ ( * - V ( * ) + £i (* _ 1 )e(* + 3) (4-71) where Qj = £jB. The superscript ^ denotes filtering by C(z~x)~x. As £j(z~x) is a degree of jf — l(j > 1), the noise terms £j(z~x)e(t + j) are all future unknown values. Therefore, the above output predictor can be modified as: y(* + il*) = Sii*-1) A u'(t + j - 1) + ^ ( z - V W (4.72) To separate the past known filtered control actions from the current control actions yet to be determined, the following polynomial matrix identity is considered: Gjiz-1) = gr(z-x)C(z-x) + z-'CAz-1) (4.73) Substituting Equation 4.73 into 4.72 gives: y(t + j\t) = ^ ( O A u ^ i - l l + ^ l A a ' f t - l l + ^ r V W = 0j(*-i)Au(t + j - l ) + /(* + i ) where /(* + j) = C^z-1) A u'(t - 1) + ^ " V C ) (4.75) M I M O Predic t ive Con t ro l L a w The optimal movements of control actions are selected such that the following desired performance index is minimized: minu J = E f i ^ [y (* + j) ~ w(t + j)]TR[y(t + j) - w(t + j)} + Z™AuT{t + j-l)QAu{t + j - l ) where Chapter 4. Constrained Multivariable Model-Based Predictive Control 89 • R is the output weighting matrix with positive elements on its diagonal, i.e., R diag{rx,r 2 ,-- - , r n } . • Q is the control weighting matrix with positive elements on its diagonal, i.e., Q diag{qi,q2,--- ,qm}. • w(t + j) is the output setpoint vector of dimension n. Substituting Equation 4.74 into 4.76 gives: min J = [Qu + f - w]TR[<7u + f - w] + uTQii (4.77) where u = [Au(t) A u(t + 1) • • • A u(t + NU - 1)]T f = [fl(t+l)---fn(t + l)---f1(t + N2)---fn(t + N2)]T w = [Wl(t + 1) • • • wn(t + 1) • • • W!(t + N2) • • • wn(t + N2)]T and Q is the matrix of dimension n(N2) x m(NU) in the case of Nx — 1: Q = Q'i 0 Si NU-l 0 0 "G'x & Gm-l G'N2-2 ' ' ' G'N2-NU Minimizing Equation 4.77 with respect to ii yields: (4.78) (4.79) Au(t) = [im o • • -}[gTRG + Q]-%g(w - f) (4.80) Chapter 4. Constrained Multivariable Model-Based Predictive Control 90 4.3 Recursive Parameter Estimation in Adaptive Control As mentioned previously, the process under this study possesses time-varying and non-linear characteristics. The best solution to compensate for the time-varying nature is to identify the changes and make an adjustment in the process model. This leads us to considering an adaptive version of a model predictive control incorporating a proper parameter estimation method. By selecting an adaptive scheme, the nonlinear problems can be treated as a linear one with a time-varying gain. Therefore, a more complex nonlinear control algorithm can be avoided. Parameter estimation incorporated with an adaptive control includes various recursive identification algorithms [95]. The most wildly used scheme of these algorithms is the Recursive Least Squares (RLS). The algorithm is: $t = yt- 4>J6t-i §t = dt^ + Ptfait Pt = A " 1 Pt-i (4-81) where <j>J is the vector of past inputs and outputs, 6t is the vector of the model parameters estimated, Pt is the square matrix proportional to the covariance of the parameters, et is the prediction error and A is the forgetting factor. The basic RLS is known to have optimal properties when the parameters are time invariant, but it is unsuitable for tracking time-varying parameters since the algorithm gain converges to zero. Modified versions of the algorithm have been developed under the effort of many researchers. A n improved algorithm is Exponential Forgetting and Resetting Algorithm (EFRA) algorithm proposed by Salgado et al. [96]. The algorithm is expressed as: et = yt~ <fIh-\ Chapter 4. Constrained Multivariable Model-Based Predictive Control 91 1 + <t>t Pt-x<Pt 1 CxPt^MjPt-l ^j. , p 2 A 1 + 0t -Pt-10t where the estimator's parameters a, (3,8 and A are constants. By making particular choices for a,/?,8 and A, the algorithm is reduced to many of the previous RLS based algorithms. Typical values of a, (3,8 and A suggested in [96] are a G [0.1,0.5], (3 G [0,0.01],8 G [0,0.01] and A G [0.9,0.99]. With these parameter ranges, the covariance matrix Pt is bounded to P m i n < Pt < P m ax- This algorithm incorporates both exponential forgetting old data to ensure continual parameter adaptation and exponential covariance resetting during periods of low excitation. Complete discussions on various parameter estimation algorithms for the linear regression type model can be found in Ljung [77] and Astrom and Wittenmark [97]. 4.4 Control under Constraints Constraints that are strictly enforced are referred to as hard constraints. Hard constraints are often taken into account when undesirable and unrealistic control inputs are computed because of process changes. These changes may happen during setpoint change or more generally when the process dynamic varies and leads to the control inputs which cannot be actually applied due to the reasons of safety, equipment physical limitations or product quality requirement. In addition, the economic operating condition of a typical process unit often lies at the intersection of constraints [98]. A successful industrial controller should, therefore, be able to maintain the process as close as possible to the constraints but without violating them. The unconstrained control criterion described earlier is to calculate a sequence of future control actions by minimizing a desired performance index. In the case when the Chapter 4. Constrained Multivariable Model-Based Predictive Control 92 controller which gives the minimum variance requires too large a control signal, it will result in the process operating out of its pre-specified operating limits. Although the control actions can be reduced by choosing relatively high control weights in the G P C algorithm, it may not be easy to find an accurate value to guarantee pre-specified variable constraints. This motivated the current study on the development of a constrained control strategy that would be suitable for the control of a T M P plant. The typical constraints for a SISO case may include amplitude/rate limits on the control signal and amplitude limits on the process output, i.e. " m i n < U(t + i - 1) < " m a x A « m i n < Au(t + i - 1) < A w m a x Vrnin < V(t + j) < ymax (4-83) i = l,---,NU j = N1,---,N2 where ( u m i n , u m a x ) '• the lower and upper bounds on the control amplitude. ( A u m i n , A u m a x ) : the lower and upper bounds on the control increment. (Vmin, 2/max) '• the lower and upper bounds on the process output amplitude. In GPC-based control development, the process is assumed to be linear. Hence any constraints on the future inputs and outputs (Equation 4.83) can be mapped to a linear combination of the future control increments according to: « m i n <u(t + i-l) < " m a x 4 (4-84) " m i n — u(t + l — 2) < Au(t + i — 1) < W m a x ~ u(t + i — 2) Chapter 4. Constrained Multivariable Model-Based Predictive Control 93 and 2/min 2/max 2/min < y(t + JVi + 1) < 2/max Ymax 2/min y(t + N2) Z/max (4.85) Ymin < y = G ' i i + f < y m a x However, it must be kept in mind that the mapping of process constraints depends greatly on the accuracy of the process model. Only an accurate model will result in the required effect as mapping of output constraints to the constraints on control increments u . Substituting Equations 4.84 and 4.85 into Equation 4.83 gives a more compact form as follows (see Appendix C for the details): A, INUXNU —INUXNU G' - G ' i i < Umax - A2u(t - 1) - u m i n + A2u(t - 1) A u m a x - A u m i n ymin f Ymax "t" f (4.86) where i i = [u(t) u(t + 1) • • • u(t + NU- l ) f (4.87) Chapter 4. Constrained Multivariable Model-Based Predictive Control 94 and A1 = 1 0 1 1 0 0 1 1 ••• 0 (4.88) A2 = [l 1 - - -1] T The elements in the vectors u m a x , umi n, Au m a x , Au m i n , y m a x and y m i n are associated values of the input and output constraints. If none of the constraints in Equation 4.86 is violated, the control actions derived by minimizing a desired performance index alone are the optimal solution. If however one or more constraints in Equation 4.86 are violated, the solution will no longer be optimal. The optimal solution should be obtained by minimizing the performance index with satisfaction of the input-output constraints, i.e. mina J = [G'u + f - w]Tii:[G'u + f - w] + uTQu subject to: (4.89) ;4i Umax - A2u(t - 1) - U m i n + A2u(t ~ 1) INUXNU u < A u m a x INUXNU - A U m m G' Ymin f - G ' Ymax "I" f The similar mathematical formula can be derived for M I M O cases. The conceptual struc-ture of the above constrained model-based predictive control strategy is also illustrated in Figure 4.29. Chapter 4. Constrained Multivariable Model-Based Predictive Control 95 objective J constraints disturbances setpoint + Optimizer Process output Model Update Output Prediction Figure 4.29: Basic structure of constrained model-based predictive control There is no analytical solution available for the general case of constrained predictive control [99]. Some analytical optimal solutions exist only for some special cases [100, 99]. The solution for a general case may be obtained with the help of some available commer-cial optimization software packages, but it turns out to be computationally demanding and may not be suitable for on-line implementations. The following presents the deriva-tion of an analytical solution for the constrained G P C for a simple case. A n optimal solution of a constrained G P C for the general case is derived by solving a mixed-weights least-squares (MWLS) problem. 4.5 Constrained Control via Quadratic Programming A n optimal solution of a constrained control problem was first developed by using linear programming, but the solution of the constrained problem was obtained using quadratic programming [101]. There are several algorithms that exist for solving QP problems. They can be roughly divided into two groups: those based on the computation of Kuhn Chapter 4. Constrained Multivariable Model-Based Predictive Control 96 Tucker multipliers [102] and Lagrange multipliers [103], and those reducing a QP problem to non-negative least squares [104] or to a linear complementary problem [105, 103]. Constrained SISO G P C for the case of NU < 2 was treated by Tsang and Clarke [100] using the Lawson and Hanson method [104], where the control signal for the rate or amplitude constraint (but not both) was calculated. Tsang and Clarke concluded in [100] that the Lawson and Hanson method is rather complex and computationally demanding, so it is not adequate for fast processes application. Camacho [106] transformed the G P C QP problem to a linear complementarity problem and solved it using the Lemke's algorithm. In Camacho's work, a simple algorithm to decrease the number of constraints and a modification to a standard algorithm for solving the linear complementary problem was used in order to reduce the amount of computation required. Recently, an analytical solution to the constrained G P C was derived by Mutha [99] through the computation of Kuhn Tucker multipliers for NU < 2 in a SISO G P C case and for NU = 1 in a M I M O G P C case. However, Mutha's work did not take output constraints into consideration in the controller development. In the following, an analytical solution to the constrained G P C subject to both input and output constraints is derived by modifying Mutha's method. 4.5.1 Quadratic Programming via Kuhn Tucker Multipliers Computation Kuhn Tucker multipliers method can be used to solve the minimization of a nonlinear quadratic cost function with linear constraints as shown below [102]: mmf(x) (4.90) subject to: hi(x) < 0 for i = 1 to p 9j(%) =2 for j — 1 to q (4.91) Chapter 4. Constrained Multivariable Model-Based Predictive Control 97 where f(x) is assumed to be continuous and continuously differentiable up to 1st (for the 1st order conditions) or 2nd order (for the 2nd order conditions) in the feasible region V which is defined by p + q constraints, x E Rm is said to be the "regular point" of the constraints if at this point all the constraints are independent. The above optimal problem can be rewritten as follows by including the constraints as an unconstrained problem: X,p) = fix) + i2Xj9j(x) + (4-92) j=i t=i where fa are Kuhn-Tucker multipliers and \j are Lagrangian multipliers. For $ to be minimum, the following 1st order and 2nd order conditions must be satisfied. The 1st order conditions are: ox M = <«3> ffh(x) = 0 The 2nd order conditions are: fa > 0 for i = 1 to p (4.94) | ^ is non-negative definite in the subspace of the constraints Note that the sign of the Kuhn-Tucker multipliers fa are subject to change with the nature of the problem (see Table 4.12 for different cases). The above conditions do not enable us to compute the values x directly since they are in general nonlinear and unsolvable analytically. Conventional solutions for the prob-lem may be via Iterative Search Methods [102], but it turns out to be computationally demanding. Hence the use of an analytical solution to solve the above problem is highly desirable as it requires reduced computation. Chapter 4. Constrained Multivariable Model-Based Predictive Control 98 Nature of the constraints Nature of the extremes Max f(x) Min f(x) subject to hi(x) < 0 subject to hi(x) > 0 Pi <o IM>0 Pi > o Pi <o Table 4.12: Kuhn-Tucker Multipliers for Different Cases 4.5.2 A n a l y t i c a l Solut ion for Constra ined SISO G P C Constrained SISO G P C for NU — 1 is described as below: m i n A u ( t ) J = [G' A u{t) + f - w ] T [G' A u(t) + f - w] + p A u2(t) subject to: «min <u(t) < "max A Umin < Au(t) < A u m a x 2/min Z/min <y = y(i + iV i + 1) < 2/max 2/min y(* + iv2) 2/max (4.95) where y ^ n = [Vmin, Vmm, ymm]T and y m a x = [ y m a x , y m a x , • • • , y m a x ] T are the vectors including the lower and higher bounds of the process output. In GPC-based control, as stated earlier, the process is assumed to be linear. Hence the input and output constraints in the above equation can be mapped as a linear combination of the control increment Au(t) by substituting u(t) — Au(t) + u(t — 1) and y = G ' A u(t) + f into the above equation, i.e.: "min - U(t - 1) < Au(t) < U m a x - u(t - 1) A t t m l n < Au(t) < A u m a x (4-96) G'T(ym i n-f) < A ( f ) < G'T(ym a ] [-f) G' TG' — ^"Vv — G' TG' As the optimal problem in Equation 4.95 involves only one variable Au(t), the mini-mization of the cost function subject to the constraints is solved by simply clipping the Chapter 4. Constrained Multivariable Model-Based Predictive Control 99 unconstrained solution Au*(t) to the appropriate bounds, i.e. a if Au*(t) < a Au(t) = I Au*(t) if a < Au*(i) < b b if 6 < A«*(t) where a and 6 are the scalars associated with the constraints, and given by: G ' T ( v • -(4.97) a = max < A i t m i n , i t m i n - u(t - 1 ) , G' TG min < A i t , max) "max (4.98) ^ G' T(ym a x-f)) - « ( * - ! ) , G' TG' / The analytical solution for NU = 2 in a SISO G P C case can be derived similarly as the derivation of the analytical solution for a constrained M I M O G P C for a NU = 1 case (see Section 4.5.3). 4.5.3 Analytical Solution for Constrained MIMO GPC Constrained G P C for NU = 1 in an m-input n-output case can be described as below: minA«(t) J = [QAu(t) + f - w]T R [SAu(t) + f - w)] + Au(t)TQ A u(t) subject to: Umin < U ( t ) < Umax Au(*)min < Au(t) < A u m a x (4.99) Y = •*• min ymin y(* + i ) Ymax Ymin < Y = y(* + 2) < Ymax : ymin _ y ( * + iV2) _ Ymax where u(t) = u2(t)---um(t)}T Au(t) = [Aui(t) A u2{t) ••• A um(t)]T y(t) = [yi(t)y2(t)---yn(t)} (4.100) Chapter 4. Constrained Multivariable Model-Based Predictive Control 100 Obviously, the solution obtained by simply clipping the control signals to their associated bounds is no longer optimal. To obtain an optimal solution, the quadratic programming is used based on the computation of the Kuhn-Tucker multipliers. After substituting u(t) = A u ( f ) + u ( t - l ) and Y = QAu(t)+f into Equation 4.99, the inequality constraints are mapped to a linear combination of the control increment Au(i) as follows: Umax iy (4.ioi) U m i n - U(t - 1) < Au(f) < U m a x - u(t - 1) and Y < Y < Y x min - x A max 4 (4-102) ( < ^ ) - i £ T ( Y m i n _ f) < Au(t) < (^ T ^) -^ T (Y m a x - f) Combining Equations 4.101 and 4.102 with A u ( t ) m i n < Au ( t ) < A u m a x gives: ax = Au(t) - a > 0 a2 = - A u(t) + 0 > 0 where a and j3 are the vectors of the dimension m, representing m constraints, and given by: a = max { Aumin, u m i n - u ( i - 1), ( ^ T ^ ) - 1 ^ T ( F m i n - f) j 1 J (4.104) P = min { A U m a x , Umax - U{t - 1), ( ^ ) - ^ T ( y m a x ~ f)} The vectors ax and a2 in Equation 4.103 are the linear combination of the constraints, defined by: On = [an a12 • • • aXr l T (4.103) lm (4.105) a2 = [a2i a22--- a2m]T After substituting Equation 4.103 into the objective function in Equation 4.99, the constrained control problem can then be solved as an unconstrained control problem as: m i n J = [SAu(t) + f - w] TR[<?Au(t) + f - w)] Au(t) Chapter 4. Constrained Multivariable Model-Based Predictive Control 101 + A u ( t ) T Q A u(t) + ^ T a x + y£a2 (4.106) The above equation can be modified as the following form: min J = A u T ( £ ) ( £ T R £ + Q)Au(t) + 2(f - w ) r R £ A u ( t ) + £ax + £a2 (4.107) Note that the term (f — w ) T ( f — w) was dropped from the above equation as it is inde-pendent of Au(t) . The associated Kuhn-Tucker multipliers for the constraints in Equa-tion 4.103 are: H2 = [A*21 A*22 • • • j"2m]T (4.108) where p and p2 are the vectors of dimension m x 1. To obtain the optimal solution Au(t) of Equation 4.107, the following 1st and 2nd order conditions have to be satisfied. 1st order conditions: d = 2(QTR g + Q) A u(t) + 2(f - w ) T R Q + £ l m + £(-In) = Q (4.109) 2nd order conditions: a i ' j f l i ' j = ° (4.110) * = 1,2, j = l , - - - , m such that an = 0 if an < 0 3 ^ J (4.111) aij > 0 if pij = 0 The optimal solution Au(t) in Equation 4.107 is obtained when the 1st order and 2nd order conditions (Equations 4.109 and 4.110) are satisfied. Note that the 1st order conditions in Equation 4.109 are a vectorized representation of m scalar equations. The steps involved in the calculation of Au(£) are based on Mutha's method [99]. These steps are: Chapter 4. Constrained Multivariable Model-Based Predictive Control 102 Step 1: Determine an unconstrained control solution Au*(t). Step 2: Test to determine if any of the constraints in Equation 4.103 are violated. If all dij > 0 (no constraints violation) then go to Step 5. Step 3: Saturate a constraint and verify if the solution satisfies the conditions of Kuhn-Tucker multipliers method, i.e. • assume any one or two (say and dji) of the constraints are satisfied by the constrained optimal solution, which gives a value Au(t). Therefore, all other /i's other than and fiji are zero. • calculate the Kuhn-Tucker multipliers and fj,ji with the help of Equa-tion 4.109, and verify if the Kuhn-Tucker multipliers u-ik and u-ji are < 0 If the Kuhn-Tucker multipliers are negative, then go to Step5. Step 4: Repeat Step 3 for the next (set of) constraint(s). Step 5: Implement the control action and repeat Steps 1-5 at next sample period. 4.6 Constrained Control via Mixed-Weights Least-Squares Algorithm The constrained control solution in the previous section is obtained by solving a quadratic programming problem. The problem stemmed from a QP-based optimal solution calcu-lation is that the shape of the set of feasible solutions is an irregular volume limited by hyperplanes defined by the constraints [107]. Furthermore, if the problem has no feasible solution, a separate procedure must ensure that a reasonable compromise is achieved be-tween the unsatisfied constraints [107]. The above problems can be avoided by applying an alternative mixed-weights least-squares (MWLS) algorithm. M W L S was proposed by Chapter 4. Constrained Multivariable Model-Based Predictive Control 103 Rossiter and Kouvaritakes [108] based on a simple modification of the Lawson's algorithm (the oo-norm minimization). A least-squares (LS) algorithm is a tool normally used for the approximation of a data set by a finite set of functions using the sum of squares of the errors as a measure of the goodness of the fit. LS has been used in almost every scientific field since it was first developed in the early 19th century by Gauss [109]. Its popularity is due in part to the fact that it can be easily understood and used. As a problem-solving technique, LS minimization is always convenient but not always meaningful or even desirable, because although the sum of the squared deviations may be small, there is no guarantee that one or more of the summands is small [110]. To solve the problem, Lawson [111] proposed a numerical procedure that consists of solving a weighted least-squares problem and recursively modifying the weights such that the larger errors receive higher weights. It was shown in Lawson's paper [111] that the algorithm always converges to a solution that minimizes the largest error (i.e. the oo-norm). Later, Lawson's algorithm was analyzed by Rice [112] and extended by Rice and Usow [113] to include the cases p G [2, oo]. Recently, Lawson's algorithm was modified by Rossiter and Kouvaritakis [108] to solve control problems with input-output constraints. Because Rossiter and Kouvaritakis method partitions the error vector into two components, one that takes the control objective into account, the other which is for the constraints, this method is referred to as the mixed-weights least-squares algorithm. Correspondingly, the method uses two different recursive relationships for the update of the weights. Very recently, the algorithm of Rice and Usow [113] was extended by Gendron [110, 107] to include the cases p G [l,oo]. Addition to QP based optimal solution calculation, the mixed-weights least-squares algorithm is used in the following to solve predictive control problems with input and output constraints. Using M W L S method, a control performance index defined by desired control strategy such as G P C can be easily augmented and the weights can be modified Chapter 4. Constrained Multivariable Model-Based Predictive Control 104 in a manner that encompasses both the requirement for the future control movements to lie inside the feasible region, and to minimize the control performance index. If the constrained optimization problem is unfeasible, the M W L S will converge to the point that minimizes the 'maximum constraint violation'. Exact theoretical optimality can be obtained only as the number of iterations within the M W L S tends to infinity [108], but the solutions that are optimal within the tolerances of practical implementations can be obtained after a few iterations only. Therefore, there is a significant reduction in the computational burden, especially when compared to the use of QP. A significant advantage of using the M W L S procedure over QP is that by minimizing an appropriate co-norm, it minimizes the likelihood of future constraint violations, and in certain cases overcomes unfeasibility problems altogether. Furthermore, the M W L S application is simple to understand and easy to implement. It can be used to solve constrained control problems for more general cases and with more flexibility. 4.6.1 Preliminary Linear p-norm problems Squared errors in solving an approximation problem is the special case of a p-norm problem. A p-norm problem is described as follows: (4.112) and its distance functional is thus: (4.113) where e/c = Vk ~ Vk (4.114) Chapter 4. Constrained Multivariable Model-Based Predictive Control 105 In Equation 4.113, || • | | p denotes the Euclidean p-norm. The vector Y — [yx y2 • • • V N ] T is the set of observations to be approximated. Y = [yi y2 • • • yN]T is the set of approximates to Y. If the approximation is considered as linear in parameters, the approximation may be expressed by: Y(e) = x{e) = Ylj^e (4.ii5) k=i where 6 is a M-dimensional vector of adjustable parameters. X — [xx x 2 • • • X J V ] t is a N x M matrix. || Y — Y \ \ in Equation 4.113 is a function that measures the distance between members of Y and members of Y. The approximation is judged according to the minimization of this measurement. Minimizing the p-norm of the distance in Equation 4.113 gives [110]: V « = E M P E N p - 2 e * V * e* = 0 (4.116) \fc=i / fe=i If assuming ^ 0, above equation becomes: E M p - 2 e f c V * e f c = 0 (4.H7) Least-squares algorithm - a linear 2-norm problem A least-squares problem is the special case of p-norms with p = 2. A least-squares or 2-norm problem can be expressed as follows by setting p = 2 in Equation 4.113: m i n l i y - y ^ H 2 (4.118) Substituting Equation 4.115 into the above equation and minimizing with respect to 6 gives: 0 = {XTX)-1XTY (4.119) or 0 = ( E x f c x f ) (4-120) \fc=i / k=i Chapter 4. Constrained Multivariable Model-Based Predictive Control 106 Weighted least-squares algorithm - a linear oo-norm problem Solving a 2-norm problem ( Equation 4.118) is to obtain the minimization of the sum of the squared approximation errors. But obviously, it does not guarantee that one or more of the summands is small. To improve this, Lawson [111] proposed weighted least-squares by attaching weights to the error vector in Equation 4.118, i.e., E = Y — Y(6). The associated weighted least-squares algorithm is described by: min\\W1/2(Y -Y(6)\\22 (4.121) 0 where W is a diagonal matrix with positive real diagonal elements, i.e., W = diag{wiW2 • • • n Computing the gradient of the above equation with respect to 6 and setting it to zero yield: § = {XTWX)-1XTWY (4.122) or alternatively / N \ ~ X N 6 = I E WkX-kX-l YI wkXkVk (4.123) \fc=i / fc=i Lawson suggested: with identity matrix, minimizing the objective function of Equa-tion 4.121 for the sequence of weighting matrices defined recursively by W t 1 ] = v n n J = l - N (4-124) Lk=i wkk \ek I where e ( i ) =y- x T 0 ( i ) (4.125) So, Equation 4.122 can be rewritten as: 6 = (XTW{i)X)-1XTW{i)Y (4.126) With increasing i, 6 will converge to the solution 6* [111, 108], which minimizes the largest member of the Y — Y(6) (hence the oo-norm) because the weights are chosen so Chapter 4. Constrained Multivariable Model-Based Predictive Control 107 as to place a greater penalty on the member of the error vector which has the largest error. 4.6.2 Constrained Control via Mixed-Weights Least-Squares Constraints In G P C based control, input-output constraints may be considered over the control hori-zon NU and predictive horizon N2. To apply the mixed-weights least squares algorithm, the constraints in the following are mapped in the form of a linear combination of control increments. Input amplitude constraints: "min < U(t) < "ma; "min < U(t + 1) < ttD (4.127) "min < U(t + NU~1)< "max where ("mim "max) are the lower and upper bounds on the control amplitudes. The above constraints can be modified in the following form: |"(^ ~\~ % 1) "centre | j* ^ radius ~ 1 i = 1 • • • NU (4.128) where "max "I" "min "'max "'min / . , nQ\ ^centre — n ^radius — ~ ^ 4 . 1 Z y j After substituting Au(t + i - 1) = u(t + i - 1) - u(t + i - 2) into Equation 4.128, the amplitude constraints on the control actions are mapped in the form of a linear Chapter 4. Constrained Multivariable Model-Based Predictive Control 108 combination of the future control increments: Au(f) Ucentre— u(t— 1) ^radius ^radius Au(t)+&u(t+l) Ucentre—u(t— 1) uradiuB uradiuB < 1 < 1 (4.130) Au(t)+Au(t+l)+---+Att(t+iVr/-l) _ u c e n t r e - u ( t - l ) "^radius uradius < 1 The above equation can be modified in a more compact form as follows: \Luu - Cu < 1 or expressed as an oo-norm problem: (4.131) \\Luu - C I U < 1 (4.132) where u = [Au(t) A u(t + 1) • • • A u(t + NU - 1)]T. Lu is an NU x iVC/ lower diagonal matrix whose non-zero entries are all equal to l/wradm«- Cu is an iV[/-diniensional vector whose entries are all equal to [ucentre — u(t — l)}/uradius. Input rate constraints: \Au(t)\ < A M |A«(* + 1)1 < Au 1 V n ~ (4.133) \Au{t + NU- 1)| < Au The above equation can be modified as: Au(t)/Au < 1 Au(t + 1 ) / A u < 1 (4.134) Au(t + NU -1)/Au < 1 The above equation can be easily rewritten in the convenient form as follows: L d u i i — Ojvuxil < 1 (4.135) Chapter 4. Constrained Multivariable Model-Based Predictive Control 109 or expressed as an oo-norm problem: | | L d u U — OjVf/xiHoo < 1 (4.136) where Ldu is an NU x NU diagonal matrix whose non-zero entries are all equal to 1/Au. Output amplitude constraints: (4.137) 2/min y(t + i) Umax 2/min < y(t + 2) < 2/max 2/min _ y(t + N2) _ 2/max The above output constraints can be modified as: \y(t + j)- ycentre\/yradius < 1 j = l,---,N2 where ycentre 2/max "I" ymin yradius y> max ymin (4.138) (4.139) 2 £7 T U U ' t ' u a ^ As derived earlier, the output prediction function for a SISO G P C case is expressed by: y = G'ii + f (4.140) To map the output constraints into the form of a linear combination of the future control increments, substituting Equation 4.140 into 4.138 gives: Lyu — Cy\ < 1 (4.141) or expressed as an oo-norm problem: \LyVt. — Cy oo < 1 (4.142) where Ly = G'/yradius is an N2 x NU matrix. Cy = ( y c - i)/yradius is an N2 x 1 matrix. y c is an N2 x 1 matrix whose entries are all equal to ycentre-Chapter 4. Constrained Multivariable Model-Based Predictive Control 110 Combining the constraints in Equations 4.132, 4.136 and 4.142 gives: where | | L u Lu Cu L = C = 0 Ly (4.143) (4.144) 4.6.3 C o n t r o l L a w The control problem for SISO G P C subject to the constraints in Equation 4.143 can be re-expressed as: m i n a J = [G'u + f - w] T [G 'u + f - w] + p u T u subject to: (4.145) | | L u - C||oo < 1 The above control performance index can be modified in the form of a 2-norm minimiza-tion as: m i n J = I I R u - V I where G ' w - f R = V = P1/2INU (4.146) (4.147) R is a matrix whose dimensions is are (N2+NU) x NU and V is a vector of the dimension N2 + NU. Chapter 4. Constrained Multivariable Model-Based Predictive Control 111 The minimization of the performance index (Equation 4.146) subject to the con-straints (Equation 4.143) is calculated by solving the following mixed-objective (or mixed-weights least-squares problem) function: 1/2 JMWLS = min W(») R L u V c (4.148) where is a positive scalar weight associated to the control performance index and is a positive definite diagonal matrix associated with the constraints, i.e., = diag {tun • • • WJJ • • • w$}. The weights are updated with: w(i) W (*+l) -2-fe=l wkk\ek W (*+l) 33 (*) (») 2^fc=l wkk ek (4.149) where iw^ are the entries of W^K The initial values of and are normally set to 33 identity. In the above equation, m is the total number of constraints (TO = 2NU + N2) and ej^ are the entries of the vector = [e^ e2^ • • • e^]T given by: e« = Lii ( i ) - C (4.150) The control movement in the above equation is obtained by minimizing the mixed objective function in Equation 4.148 with respect to ii, which is: (\ u (0 R L R L R L W(») V c (4.151) Repeat the procedure with increasing i until the algorithm converges. A criterion based on the relative variation of the solution from one iteration to the next is calculated to stop the iterations. The criterion is: , lluW 5' = (4.152) Chapter 4. Constrained Multivariable Model-Based Predictive Control 112 > s* continue if sl { (4.153) [ < s* stop where s* is some threshold. By recursively using the procedure with increasing i, the algorithm can only converge to the vector i i which minimizes the mixed objective function in Equation 4.148 and it always converges to either [108] depending on whether a feasible solution exists or not. The proof of this convergence and a discussion of its associated properties is given in [108]. Because the algorithm uses two different recursive relationships for the update of the weights, it is referred to as the mixed-weights least-squares algorithm. 4.7 Simulat ions The purpose of simulations was to demonstrate the effectiveness of the proposed con-strained G P C (described earlier) with the application of the quadratic programming and mixed-weights least-squares algorithms. A single-input single-output linear plant was first selected for implementing the control algorithm. A 2 x 2 linear plant was then selected for further investigation of the proposed control algorithm. Example 4.1: Consider the single-input single-output second-order linear plant: In the following simulations, GPC-based controller tuning parameters were set to Ni = 1,N2 = 10, NU = 1 and p = 0. A l l the initial parameter estimates were 0.1. The initial covariance matrices for the RLS algorithm were 1000 x / , and all the forgetting factors UQP or m in | |Lu — C | | u oo (4.154) y(t) - 1.5y(t - 1) + 0.6y(t - 2) = 0.006u(t - 1) + 0.06«(t - 2) (4.155) Chapter 4. Constrained Multivariable Model-Based Predictive Control 113 det u __ 1 1 V V 1 I i " i \ li! I 1 1 0 50 100 150 200 250 300 350 400 Sampling Interval Figure 4.30: Constrained G P C via analytical QP (t > 300 sampling intervals: -1.5 < ut < 1.5, -2.0< Aut <2.0) Chapter 4. Constrained Multivariable Model-Based Predictive Control 114 det u 1 .... ! I V V i 0 50 100 150 200 250 300 350 400 Sampling Interval Figure 4.31: Constrained G P C via M W L S ( t > 300 sampling intervals: -1.5< ut < 1.5, -2.0 < Aut < 2.0) Chapter 4. Constrained Multivariable Model-Based Predictive Control 115 were 0.99. Simulations were performed with and without constraints. The output setpoint was square waves of amplitudes 1 unit. Input Amplitude and Rate Constraints: Figures 4.30 and 4.31 show the behav-iors of the constrained control through solving the quadratic programming and mixed-weights least-squares problems, respectively. To make a comparison, the standard uncon-strained GPC was simulated for the first 300 sampling intervals. The constrained control algorithm was then turned on at sampling interval 300. The maximum and the minimum amplitude constraints on the control input were set to +1.5 and —1.5. The maximum and the minimum rate constraints on the control input were set to +2.0 and —2.0. As can be seen, the proposed control algorithms based on both quadratic programming and mixed-weights least-squares worked well since the rate and amplitude constraints on the input were successfully limited to their specified ranges without a sacrifice of setpoint tracking. For the constrained control based on the mixed-weights least-squares algorithm, we set the threshold for the stopping of the algorithm to 0.001 (i.e. s* = 0.001 in Equation 2.24). After about 20 iterations, the solution of minimizing the mixed objective function was obtained and the iterations stopped. Example 4.2: Consider the 2x2 second-order linear plant: yi(t) = 1.57yi(t-l)-0.615j/i(t- 2) + 0.002ui(t-l) -0.00213«i(t - 2) + 0.072w2(t - 1) - 0.054u2(t - 2) y2(t) = 1 .57y 2 (t- l ) -0 .615y2(*- 2) + 0 . 0 7 0 2 « i ( t - l ) -0.053ui(t - 2) + 0.03w2(i - 1) - 0.0246u2(t - 2) (4.156) For all the simulations in the following, GPC-based controller tuning parameters were set to Ni = 1, iV2 = 10, NU = 1 and px = p2 = 0. Simulations were performed with and without the constraints. The output setpoints were square waves of amplitudes 1 unit. Chapter 4. Constrained Multivariable Model-Based Predictive Control 116 1n 1 1 1 1 r 0.9- I 0.8- ; I 0.7- : j 10.6- ; Iteration Figure 4.32: Criterion for stopping the iteration Input A m p l i t u d e and Rate Constrained: Figures 4.33 and 4.34 show the be-haviors of the constrained G P C with input amplitude/rate constraints using quadratic programming and mixed-weights least-squares algorithms, respectively. In Figures 4.33 and 4.34, the standard unconstrained G P C was used for the first 300 sampling intervals. The setpoint tracking of this period simulation was ideal. But some control moves were relatively larger in amplitude and rate. To reduce the amplitude and rate of the control inputs, the proposed constrained G P C was used for the next 100 sampling intervals with -4.0 < «i < 4.0, -4.0 < At*i < 4.0, -4.0 < u2 < 4.0 and —4.0 < A u 2 < 4.0. As can be seen in Figures 4.33 and 4.34, the input amplitude and rate were limited to their specified ranges for the sampling interval > 300 without a sacrifice of setpoint tracking. For the constrained control based on M W L S , we set the threshold for the stopping of the algorithm to 0.001 (i.e. s* = 0.001 in Equation 2.24). After about 50 iterations, the solution of minimizing the mixed objective function was obtained and the iterations stopped. Figure 4.32 shows the feasible value of s* as a function of the iteration number. Chapter 4. Constrained Multivariable Model-Based Predictive Control 117 Outpu t Constraints : Figures 4.35 and 4.36 show the behaviors of the output con-strained G P C using quadratic programming and mixed-weights least-squares algorithm, respectively. The standard unconstrained G P C was used for t < 300 sampling intervals and the setpoint tracking of the simulations was ideal. To reduce the small overshoots on the first output yi, the constrained G P C was on at t = 300 sampling interval and yi was limited to —1.0 < yi < 1.0. Figures 4.35 shows that the constrained G P C via quadratic programming was successful since the output yi was limited to its desired range. To further investigate the constrained control algorithm by solving the M W L S , the maxi-mum and minimum amplitudes on both yi and y2 were set to 1.0 and —1.0. Figure 4.36 shows the time-trajectory plot of the simulation of the constrained G P C using the mixed weights least-squares algorithm. Chapter 4. Constrained Multivariable Model-Based Predictive Control 118 10 5 = o - 5 - 1 0 10 5 <3 0 - 5 ^- k r f ( 100 200 300 Sampling Interval 400 -10 100 200 300 Sampling Interval 400 •-r 100 200 300 Sampling Interval 400 10 CM -10 1 100 - i — n f r -1 200 300 Sampling Interval Figure 4.33: Constrained G P C via analytical QP (t > 300 sampling intervals < «i < 4.0, -4.0< u2 < 4.0, | A « i | < 4.0, | A u2\ < 4.0) 400 -4.0 Chapter 4. Constrained Multivariable Model-Based Predictive Control 119 100 200 300 Sampling Interval 400 100 200 300 Sampling Interval 400 10 5 ^ 0 - 5 k L i ( ( 100 200 300 Sampling Interval 400 -10'— 100 200 300 Sampling Interval 400 rr 100 200 300 Sampling Interval 400 10 CM 3 CD T 3 - 5 -10 I 100 200 300 Sampling Interval 400 Figure 4.34: Constrained G P C via M W L S (t > 300 sampling intervals: -4.0 < ui < 4.0, -4.0< u2 < 4.0, | A «i | < 4.0, | A u2\ < 4.0) Chapter 4. Constrained Multivariable Model-Based Predictive Control 120 10 8 6 4 2 0 - 2 - 4 - 6 - 8 - 1 0 ^ V T J 100 200 300 Sampl ing Interval 400 10 8 6 4 2 0 - 2 - 4 - 6 - 8 - 1 0 I 1 V V V f r I if 100 200 300 Sampl ing Interval 400 Figure 4.35: Output constrained G P C via analytical QP (t > 300 sampling intervals: -1.0 < V l < 1.0) Chapter 4. Constrained Multivariable Model-Based Predictive Control -2 100 200 300 Sampl ing Interval 400 2 1.5 1 0.5 0 -0.5 -1 -1.5 - 2 100 200 300 Sampl ing Interval 400 100 200 300 400 100 200 300 400 Sampl ing Interval Sampl ing Interval Figure 4.36: Output constrained G P C via M W L S (t > 300 sampling intervals: -1.0< y i 1.0, -1.0 < y 2 < 1.0) Chapter 4. Constrained Multivariable Model-Based Predictive Control 122 4.8 Summary and Conclusions A n analytical solution exists for the unconstrained generalized predictive control (GPC). However, no analytical solution exists for the general case of constrained G P C . The so-lution for the general case of constrained G P C may be obtained by using commercial optimization software packages, but it turns out to be computationally demanding as compared with an analytical solution. This chapter describes the derivation of a con-strained M I M O G P C subject to input and/or output constraints. A n analytical solution of the constrained M I M O G P C for the case of NU = 1 was obtained through solv-ing a quadratic programming (QP) problem by computing Kuhn Tucker multipliers. A problem resulting from a QP-based optimal solution calculation is that there will be no feasible solutions if one or more constraints are overly stringent. To avoid the problem, the mixed-weights least-squares (MWLS) algorithm was introduced to the constrained M I M O G P C to solve an optimization problem. A significant advantage of using the M W L S is that if the constrained optimization problem is unfeasible, the M W L S will converge to the point that minimizes the maximum constraint violation. In addition, the MWLS-based optimal solution is able to handle the general cases of constrained G P C with more flexibility. The proposed control schemes were tested on SISO and 2 x 2 linear models, showing the success of handling input and/or output constraints. The proposed control schemes will be used for the control of the T M P process as given in Chapter 6. Chapter 5 Laguerre Function-Based Adaptive-Predictive Control 5.1 Introduction In dealing with a model-based control approach, parametric representations, e.g. trans-fer functions such as an A R M A X or A R I M A model, are commonly used. A possible reason for this scheme is that the model is simple and the theory based on this model is well developed [114]. However, a linear transfer function-based control approach will work well only when the nonlinearity of a plant is not severe and the actual plant is represented properly by the structure of the model. In the case of severe nonlinearity, the linear approximation-based control may not be satisfactory. It may result in, in some cases, significant static errors and in other cases, oscillation or even instability [114]. Also, the linear transfer function-based approaches generally require accurate knowledge of the model order and time-delay, and can be very sensitive to a model mismatch. A l -ternatively, a non-parametric representation like the impulse response requires no other information than the linearity of the plant [115]. A disadvantage of the impulse response representation is that a large number of parameters are needed. Considered in this chapter is the use of the orthonormal functions or filter representa-tions to present the dynamics of an actual plant. By the right choice of an orthonormal filter, plant dynamics can be represented properly without the need of assumptions about plant order and time delay. A n interesting set of orthonormal filters is the set of Laguerre functions. To overcome the problems caused by severe nonlinearities, the use of nonlinear 123 Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 124 Laguerre functions representations is the nature choice. The nonlinear Laguerre func-tion is a Wiener type model, which is composed of a linear dynamic block followed by a memoryless nonlinear function. The nonlinear Laguerre function approximation-based control will be very beneficial for the control of wood chip refining when the nonlinear-ity between refiner motor load and plate gap becomes severe, i.e., the incremental gain between refiner motor load and plate gap subject to a sudden change in the sign (Figure 6.50). This is particularly critical when a refiner is operating close to the maximum load point. The linear approximation-based control for this case may lead to unstable oper-ation. Alternatively, the nonlinear Laguerre functions are used to approximate process dynamics and static nonlinearity. The Laguerre function is expressed in a state-space form, so any state-space control design techniques may be used for a state-feedback con-trol. However, the generalized predictive control algorithm is used for the controller development because of its simplicity and ease of use. The optimal control actions are calculated by minimizing future output errors and control movements. The future out-puts are predicted based on the Laguerre model under the influence of future control actions. For an adaptive control implementation, the nonlinear Laguerre model can be arranged in a linear form in the parameters so that the powerful but simple Recursive Least Squares (RLS) algorithm can be used for the parameter estimation. The use of orthonormal functions for obtaining approximations goes back as far as the development of the Fourier series [116]. However, an appropriate and well-known orthogonal basis set is the set of Laguerre filters [117]. The studies and applications of Laguerre approximations have been carried out recently by many researchers. Linear Laguerre model-based adaptive control theory was first proposed by Dumont and Zervos [118]. Their theory was then applied to pH control by Dumont, et al. [119], and indus-trial application results showed that the Laguerre model-based adaptive pH control was Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 125 superior to those using other types of models [119]. Laguerre filter-based nonlinear dy-namic approximation is based on Volterra functional series expansion [120, 114]. Early work on the Volterra model directly applies the nonlinear memoryless function to the plant input signal and its delays, resulting in a very large number of model parameters even for a simple nonlinear process [114]. However, using Laguerre filter-based nonlinear approximation the number of model parameters can be significantly reduced [120]. In a paper by Fu and Dumont [51], a Laguerre filter-based nonlinear dynamic model was used for the adaptive control of the motor load of a wood chip refiner used in the pulp and paper industry. Fu and Dumont's work was further extended by Du et al. [52] into a 2 x 2 control system based on mixed linear and nonlinear Laguerre functions for the control of both the motor load and the outlet consistency. Other studies on the theories and applications of the Laguerre approximation can be found in [116, 121, 114]. Studies and discussions of time domain properties of Laguerre approximation are detailed in [121]. Laguerre function-based representation is often used because of its flexible structure so that the model complexity can easily be changed on-line. Also, Laguerre function-based representation is very similar to transient signals due to its block-oriented feature and thus can be understood easily. First, linear and nonlinear Laguerre function-based approximations are introduced in Section 5.2. A linear Laguerre function-based G P C algorithm is derived in Section 5.3. In the same section, the limitation of a linear model-based control is described and a nonlinear Laguerre function-based control is then presented. Section 5.4 presents the technique for the estimation of Laguerre function parameters. Summary and conclusions are given in Section 5.5. Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 126 5.2 The Laguerre Functions The Laguerre functions, a complete orthonormal set in L 2[0 ,oo) , have often been used because of their convenient network realizations [117] and their similarity to transient signals [122]. In the following, a linear Laguerre model is a linear representation where the model output is a linear combination of the Laguerre filter outputs, while a nonlinear Laguerre model is a nonlinear representation where the model output is a nonlinear combination of the Laguerre filter outputs. 5.2.1 The Linear Laguerre Functions The continuous-time Laguerre functions The continuous-time linear Laguerre model has been discussed and used for the control purposes in [116, 119, 121]. In the s-domain, the i th order Laguerre function JPJ(S) is: as follows: = v^(!PK- (5-157) Any transfer function of a stable system can be described by this Laguerre ladder network using infinite number of filters as: CO G(s) = YciFi(s) (5.158) i = i In practice, the first iV filters are used to approximate the system. This approximation can be adjusted by the convergence results given in [121]. p is the Laguerre filter pole. The output of the iVth order Laguerre model is: y(t) = E c i I i ( t ) (5.159) i=l where li(t) is the i th Laguerre function output. The Laguerre ladder network is also illustrated in Figure 5.37(a) where the first block is a low-pass filter and the remaining Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 127 u(t) /2p li(t) s-p s+p s+p l2(t) IN (t) Laguerre Network Parameters (a) Laguerre model in s domain u(t) / 1 - a 2 z-a Laguerre Network Parameters 1-az z-a l2(t) 1-az z-a IN (t) C N Summing Circuit y(t) (b) Laguerre model in z domain Figure 5.37: Laguerre ladder networks blocks are all-pass niters. As can be seen from Figure 5.37(a), the simplicity of a phase-shift chain, convenient network realization and the similarity to transient signals are the significant advantages of the Laguerre function representations. The discrete-time Laguerre functions The transformation of the continuous-time Laguerre functions to the discrete ones has been discussed in [122, 123]. A discrete-time Laguerre ladder network, similar to the continuous one, is illustrated in Figure 5.37(b). In the z-domain, the i th order Laguerre filter is: -(1 -azf-1 Ftiz) = V(l-«2)- {z - ay Any stable transfer function G(z) can be expanded as: (5.160) G(z) = Y,CiFi{z) (5.161) Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 128 If the first N filters are used to approximate a process, the output of the iVth order discrete-time Laguerre model is: N y(t) = 5>fc(*) i = l (5.162) where li(t) is the ith. Laguerre filter output and Cj are the Laguerre function coefficients. The state space formula of the discrete-time Laguerre ladder network can be easily derived from Figure 5.37(b) as: L(t) = AL(t - 1) + Bu(t - 1) y(t) = CTL(t) where C = [ci c2 • • • cN}T L(t) = [h(t) l2(t)---lN(t)}T and (-a , ) ( l-a?) 0 1 -a? ( - a ^ - 2 ( l - a f ) ( - a ^ " 3 ( l - a 2 ) 0 • •• 0 0 • •• 0 •• 0 •• 0 . . . . 1 ( -Oi) ( - « i ) 2 (5.163) (5.164) (5.165) (-a*)"- 1 As can be seen from Equation 5.165, the matrix A and the vector B depend on the Laguerre filter pole a which is a free parameter. The choice of a Laguerre filter pole was first discussed in [124] and more recently in [125]. Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 129 5.2.2 The Nonlinear Laguerre Functions The nonlinear Laguerre model is a special case of the general Wiener-Volterra series representation [126], where the Wiener-Volterra functional series is truncated and the Volterra kernels are approximated using the Laguerre functions [120] with the assumption that all the kernels are in L 2 [0, oo). The continuous-time nonlinear Laguerre functions The Wiener-Volterra series representation for a nonlinear and time invariant system is: 0 0 POO POO y(t) = h0 + V / • • • / hn(rx, rn)u(t - TI) • • • u(t - Tn)drx • • • drn (5.166) where u(t) is the input and y(t) is the output. The functions hn(ri • • • r„) are the Volterra kernels which represent the nonlinear dynamic of the system. In practice, the system is often approximated by a finite Wiener-Volterra series. To simplify the notation, the above series is truncated to have only the first two Volterra kernels /ii(ri) and h,2(ri,T2). Thus we have: POO POO POO y(t) = h0+ h1(T1)u(t-Tl)dTl+ I I h2(Ti,T2)u(t - r^uit - r2)dTidT2 (5.167) Jo Jo Jo Let <f>i(t) be the ith-order Laguerre function. The i th order Laguerre filter output is: poo k(t) = / <f>i(T)u(t - r)dr (5.168) Jo Since all the Volterra kernels form a complete orthonormal set in L 2[0, oo), the truncated Volterra kernels using N Laguerre filters are expressed as [120]: N N N M r i ) = E c f e ^ ( r i ) h2{rUT2) = X ] E c n m </» n ( r i )0 m ( r 2 ) (5.169) fc=l n= l m= l where c*. and cnm are constant coefficients. Substituting the above equation into Equation 5.167 and considering the orthonormality of the Laguerre functions give: N N N y(t) = C 0 + E CMt) + E E Cnmln(t)lm(t) (5-170) fc=l n= l m = l Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 130 (5.171) The continuous-time state-space formula of the 2nd order nonlinear Laguerre model is then expressed as: L(t) = AL(t) + Bu(t) y(t) = c 0 + CTL(t) + LT(t)DL(t) where C = [ci,--- ,cN]T L(t) = [h(t),---,lN(t)}T D = C l l C12 C21 C 2 2 ClN C2N (5.172) Cm CjV2 • • • CpfM Since the Volterra kernels are symmetric, i.e., cmn — c n m , the number of unknown coeffi-cients in Equation 5.170 is (N + 1)(N + 2)/2. The matrix A and the vector B in Equation 5.171 depend only on the selected Laguerre filter pole. The discrete-time nonlinear Laguerre functions Similarly to the continuous-time nonlinear Laguerre formula, the discrete-time non-linear Laguerre model is expressed as follows by extending the continuous-time formula using the discrete-time Laguerre filters: (5.173) L(t) = AL(t - 1) + Bu(t - 1) y(t) = c0 + CTL(t) + LT(t)DL{t) + ••• where the scalar c 0, the vector C and the matrix D are the Laguerre function coefficients to be estimated. Under the condition of neglecting the higher-order Volterra kernels, the above equation is then an approximation of Wiener-Volterra series representation for a nonlinear dynamic system. The discrete-time nonlinear Laguerre network is also illustrated in Figure 5.38, where the linear dynamic part, in occurrence as the Laguerre Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 131 Laguerre Network u(t) z-a li(t) 1-az z-a l2(D 1-az z-a ls(t) ; 1st order terms 2nd order terms [1,ll(t),l2(t),l3(t),...] [Il(t),l,(t)l2(t) l2(t),...J3(t),. c r y(t) — * • Memoryless Nonlinearity Figure 5.38: The discrete-time nonlinear Laguerre network network, is followed by a memoryless nonlinear function which is a nonlinear combination of the Laguerre filter outputs. 5.3 SISO Laguerre Function-Based Predictive Control 5.3.1 Linear Model-Based Control Model output prediction Consider a stochastic process described by the following linear Laguerre model: L(t) = AL(t - 1) + Bu(t - d) y{t) = CTL(t) + m/& (5.174) (5.175) where u(t) is the input and y(t) is the output. £(£) is uncorrelated, normally distributed zero-mean random noise. The difference operator A = 1 — z - 1 in the denominator of the noise term forces an integrator into the controller in order to eliminate any offset between the measured output and its setpoint. The matrix A and the vector B are dependent on Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 132 only the selected Laguerre filter poles (Equation 5.165). If Equation 5.174 is multiplied by A with the assumption that the time-delay d = 1 (without losing the generality by adding states) and the future control movements Au(t + j) = 0, then we have: AL(t + j) = Aj A L(t) + (A*-1 -\ h I)B A u(t) (5.176) To derive a j-step ahead predictor y(t + j), consider the following identity: y(t + j) = y(t) + Ay( t + l ) + . . - + Ay(t + j) = y(t) + J2Ay(t + i) (5.177) 8 = 1 Multiplying Equation 5.175 by Azj and combining with Equations 5.176 and 5.177 gives: y(t + j) = y(t) + CTAi A L(t) + J2 ^(A*-1 H vI)B A u(t) + + *) (5-178) i = l i = l i = l As the noise components in the above equation are all in the future, the output predictor, given measured output and control action up to time t, is: y(t + j) = J2 ^{A*-1 H h I)B A u(t) + J2 c T A i A L(t) + V(t) (5-179) i=i i=i The output predictor y(t + j) consists of two parts: the first term depending on the control action yet to be determined, and the rest depending on the past control action and the measured output up to time t. Considering the output predictor y(t + j) over the predictive horizon N2 gives: y(t + l) = CTA°B A u(t) + CTA A L(t) + y(t) y(t + 2) = E 2 =i C T ( - 4 i _ 1 H r- I)B A u(t) + £ 2 = i CTA{ A L(t) + y(t) y(t + N2) = Z?IiCT(Ai-1 + --- + I)BAu(t) + z¥1CTAiAL(t)+y(t) Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 133 The above equation can be rewritten in a more compact form as: y = G A t i ( t ) + f (5.181) where G = CTA°B CT(Ai~1 + ... + I)B (5.182) and fl CTA A L{t) + y(t) h J22i=1CTAiAL(t) + y(t) f = — / A T 2 E i l a i C T ^ A L ( t ) + y(t) (5.183) Predictive control law The target of GPC-based predictive control is to minimize the following performance index: { N2 NU E [y(* + j)-w(t+j)? + E P A U ( H J - I ) 2 i=JVi j=i The above performance index can be rewritten in a more compact form: J = E {(y - w ) T ( y - w) + pu T u} (5.184) (5.185) where (5.186) y = [y(t + Nr), y(t + N1 + l),---,y(t + N2)]T w = [w(t + JVi), w(t + N1 + l),---,w{t + N2)]T u = [Au(t), Au(t + 2), • • •, Au(t + NU- 1)] T Substituting Equation 5.181 into Equation 5.185 for the case N\ = 1, NU — 1 and minimizing J yields: Au(t) = ( G T G + p ) " 1 G r ( w - f) (5.187) Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 134 Estimation of Laguerre model parameters The output of the iV-filter Laguerre network in Equation 5.175 is simply the linear combination of all filter outputs, i.e. N y(t) = 5 > * « ( t ) + £ ( * ) / A (5-188) The above equation can then be rearranged in the form: Ay(t)=6?4>t + t(t) (5.189) where the model parameters to be estimated are: 6j = [C l c 2 . . . C j v ] (5.190) and the regression vector is: 4>Tt = [Ah(t) A l2(t) • • • A lN(t)] (5.191) Since Equation 5.189 is a linear form in the parameters, the standard RLS parameter estimation algorithm can be used for the adaptive control purpose. Limitations of linear model-based control Control of a wood chip refiner in a T M P process is sometimes difficult due to the nonlinear relationship between the refiner motor load and plate gap with a sudden change in the sign of incremental gain, illustrated by the dashed line in Figure 5.39. For the nonlinearity represented by the solid curve in Figure 5.39, the linear model 1 based control will be satisfactory as long as the output setpoint yr is reachable. For this case, the process will settle at the operating point A. However, when the nonlinearity is represented by the dashed line, the linear model 2 based on the operating point B will predict that point C is the desired position. The control action will move the process to the position D. The linear model based control will then predict that the desired position Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 135 Linear model 2 Linear model 1 Linear model 3 Output Set point y r Control Input Figure 5.39: Linear model-based control will lead to oscillation when the nonlinearity is represented by the dashed curve [1] is F , but resulted control action leads the process back to B point. If no heuristic is used, the process will be in a limited cycle and can never be stabilized at the extreme point E because at that point the corresponding process gain is zero, resulting in an infinite control signal. One solution to the problem is the use of a nonlinear model-based control approach. 5.3.2 Nonlinear Model-Based Control Model output prediction Consider a 2nd order nonlinear Laguerre model using N filters: L{t) = AL(t - 1) + Bu(t - d) y(t) = c 0 + CTL(t) + LTDL(t) (5.192) Recursive use of Equation 5.192 gives: L(t + j) = AiL(t) + (A3'1 + ••• + I)Buit) = A>L(t) + A3Bu{t) (5.193) Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 136 where Aj = Aj~1 + --- + I (5.194) For simplifying the above derivation, we assumed that the time delay is one (i.e. d = 1) and the future control movements are zeros (i.e. u(t) = u(t + 1) = • • • = u(t + j)). Predictive control law Predictive control action is obtained by minimizing the performance index: N2 J = £ [w(t + j) - y(t + j)}2 + pA u2{t) (5.195) The minimization of the above performance index is done by solving: d J -2 f: [w(t + j) - y(t + j)} 9y{* + j ) + 2p A u(t) = 0 (5.196) du(t) jttr J' "K J,i du(t) or P A «(t) = E Mt+j) - vd+j)} dj^L (5-197) Substituting the output equation in Equation 5.192 into the above equation gives: p A u(t) = E + J) " Vi* + MCT + 2LT(t + j)D]AjB (5.198) The above equation will involve solving a nonlinear equation and the analytical solution for a general case is difficult to obtain. A n analytical solution for the special case of Ni = N2 = d and p — 0 is derived in [51] which is: [a u2{t) + b u{t) + c][2a u{t) + 6] = 0 (5.199) The control solution is then: -fefcvV-4oc i f b 2 _ 4 a c > 0 u(t) = { 2 a ~ (5.200) - £ i f f e 2 - 4 a c < 0 Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 137 where a - BTAjDAdB b = CTAdB + 2LT(t) (Ad)TDAdB (5.201) c = co + CTAdL(t) + LT(t)(Ad)TDAdL(t) - w(t) and w(t) = w(t + j),j > 1, i.e., the future setpoints are considered to be equal to the present one. If it is the case of b2 — 4ac > 0, one way to select the sign in Equation 5.200 is to look at how the future steps will be influenced, as suggested in [127]. However, if it is the case of b2 — Aac < 0, the solution in Equation 5.200 is corresponds to the extreme point in Figure 5.39. Estimation of Laguerre model parameters The Laguerre model output given in Equation 5.192 can be rearranged as: y(t) = c 0 + ciZi(t) + • • • + cNlN(t) + Clll2(t) + c12lx(t)l2(t) + • • • +c1Nh(t)lN(t) + c22l2(t) + ••• + cNNl2N(t) (5.202) = 0T<t>t where the vector of the model parameters to be estimated is: = [c0 ci • • • cN cu C12 • • • c22 ••• cNN] (5.203) and the regression vector is: tf = [lh ••• l N l l h h ••• l\ ••• l%] (5.204) As can be seen, the nonlinear Laguerre model output is simply in a linear form in the model parameters. Therefore, the standard RLS parameter estimation algorithm can be used for the model parameter estimation. Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 138 5.4 MIMO Laguerre Function-Based Predictive Control In the following, we use refiner control as an example to describe a multivariable Laguerre function-based predictive control. 5.4.1 Mixed Linear-Nonlinear Laguerre Model Consider a 2 x 2 mixed linear-nonlinear Laguerre functions using the first N filters for modelling a wood chip refiner: Li,t — i4i^i,t- i + BiUift_dl Hi = A2L2^X + B 2 u 2 ^ d 2 ( 5 2 Q 5 ) Vi,t = Co + CjLht + LjjD^t + ••• + fciz/2,t-i + Nltt V2,t = C2TL2,t + fayi,t-i + N2tt where u\ is the pulse sent to the hydraulic cylinder to position the gap, u2 is the dilution flow rate to control the consistency in the refining zone, y\ is the motor load and y2 the outlet consistency, di and d2 present time delays in terms of sampling unit for each input-output pair. Ni and N2 are the unmeasured disturbance and measurement noises acted on the outputs. kiy2 and k2yi present the load disturbances induced due to the changes in the motor load and outlet consistency. Matrix Ai and vector B{ depend on Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 139 only the Laguerre filter poles a;: Oj and A = 1 - a 2 ( - 0 0 ( 1 - 0 ? ) 0 1 -a? ( - a ^ - 2 ( l - a 2) ( - a ^ - 3 ( l - a 2 ) 0 • •• 0 0 • •• 0 Oi • •• 0 .. . •• 0 .. . •• Oi 1 ( - O i ) ( - O i ) 2 (~Oi) N-l Li — [In ll2 • • • ^ liv]T C\ = [Cu Cu • • • CIJV] 1/2 = [^ 21 ^22 " • " ^ 2Jv]T C*2 = [C21 C 22 • • • C2jv] (5.206) (5.207) d l x d 1 2 • • • <iiAr d 2 i d22 • • • d2jv (5.208) djvi djv2 • • • djvjv The discussions on the selection of Laguerre poles was first given in [124] and recently in [114]. The Laguerre model outputs in Equation 5.205 can be rearranged in a linear in parameter form: yi,t = 9i,t<f>i,t V2,t = 02,t<fa,t (5.209) Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 140 where Q\ = [Co,Cn,Ci2 ,dii,di2,d22 , fci] (5.210) 0i,< = [1) hi, h2i hih2, l\2iy2,t-i]T and * = [.*,<*.*] ( 5 2 n ) 4>2,t = [ h i , k 2 , y i , t - i ] T By the above arrangement, the simple but powerful RLS is readily implemented. This estimation algorithm is: @i,t = @i,t-i + Pi,t<t>i,tei,t ei,t = Vi,t - <Pj,A,t-i where Pitt are the error covariance matrix, ei<t are the predictive error and A; are the forgetting factors. 5.4.2 MIMO Laguerre Function-Based Predictive Control Law The future control actions are determined by minimizing a designed trade-off between control performance and constraints on the control movements as: f N2 J = E { E [(Vi,t+i - wht+j)2 + (y2,t+j ~ w2it+j)2] (5.213) U=JVi +Pi A u\tt + p2 A u\t) where W{ are the setpoints of the outputs and Aui>t are the control increments (Auitt = Ui,t - « i , t - l ) -Assuming time-delay di = 1 (without losing the generality by adding states) and Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 141 recursively substituting Equation 5.205 itself yields: Lltt+j = A{Lu + (A{-l + ---+I)BlUl<t 1 ' _ ' (5.214) L2,t+3 = A{L2>t + (Af1 + • • • + I)B2u2tt = A32L2tt + A2B2u2tt It was assumed in the above equation that the future control movements Auitt+j = 0, j > 0 (or u(t) — u(t + 1) = • • • = u(t + j)). To derive a j-step ahead predictor of yitt+j, the output functions in Equation 5.205 are rewritten as: Vi,t+j = c 0 + CxTLlit+j + Ljj+jDtLxj+j H h kiy2,t+j-i + £i,t+j ^ V2,t+j = C2 L2,t+j + k2yu+j-i + 6 , t + j / A where are the independent, normally distributed (white) noise sequence with variance p|. The noise term N^t acted on yitt in Equation 5.205 is assumed to be white noise to simplify the derivation. To derive a j step ahead predictor of y2,t+j, consider the following identify: V2,t+j = V2,t + Aj/2,t+i + • • • + Ay2,t+j = yu + J2Ayu+i (5.216) i=i Substituting Equation 5.214 into Equation 5.215 and combining the resultant with Equa-tion 5.216 gives: Vi,t+j = c0 + C1T[A{Lht + A1B1ulit\ + + [A{Lltt + A^u^Y DX[A{LU + I i B i « M ] + k^t+j-i (5-217) y2,t+i = y2,t + £ C 2 T [ 4 A L 2 > t + i 2 B 2 A u 2 > t ] + + E fc2 A y1>t+i-i i=l Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 142 The optimum control actions can then be determined by substituting the above equation into the performance index in Equation 5.213 and minimizing it, i.e.: — = 2 £ [ylit+j - wltt+j][C? + 2Llt+jD1}A1Bl + 2Pl A t t M = 0 (5.218) dJ N* j — = 2 £ (y2tt+j - w2tt+j) E Cj(A\-x + ••• + I)B2 + 2p2 A u2it = 0 ou2,t j = N l i = 1 Let Ni = N2 = d, the above equation can then be simplified as: pi A u M = (a «J > t + b uht + c)(2 a « 1 > t + 6) P2 A Ma,* = 7 l A «2,t + 70 where (5.219) a = B?A'(D1A1B1 b = C ^ S x + 2 L f ) t ( A f ) T D 1 I 1 S 1 c = co + C?AiLltt + Llt{Af)TDlAdL1,t + kmt+d-i - «>i 7i = - [ E t 1 C 2 T ( A 2 - 1 + --- + / ) B 2 ] 2 7o = [w2 - V2,t - E t i C 2 T 4 A L 2 , t - E t x *2 A y,^} E f = 1 C ^ " 1 + • • • + J ) B 2 (5.220) If px is set to zero, i.e., pi = 0, the above equation can be solved analytically as follows: f -6±Vy-4ac i f f t 2 _ 4 a c > 0 «i,t = 2 a " (5.221) [ - £ i f 6 2 - 4 a c < 0 v>2,t = —— \-u2,t-i (5.222) Pi - 7 i In the case of b2 — Aac > 0, one way to select the sign in Equation 5.221 is to look at how the future steps will be influenced, as suggested by Wittenmark [127]. However, for the control of the refiner load, minus sign in the above equation is used since it corresponds to the plate position in the normal operating region [51]. When b2 — Aac < 0, the only solution is u\it — ~ ^  which corresponds to the extreme point. Chapter 5. Laguerre Function-Based Adaptive-Predictive Control 143 5.5 Summary and Conclusions A linear transfer function-based control will work well only when the nonlinearity of a plant is not severe and the model structure error is not significant. In the case of severe nonlinearity and unknown dynamics of a process, the linear transfer function based control may lead to remarkable static errors or unstable operations. To avoid the problems, the nonlinear Laguerre functions-a type of orthonormal functions are used to represent an actual plant. By the right choice of an orthonormal filter, the plant can be represented properly without the need for assumptions about plant order and time delay, i.e., accurate assumptions about these true values are not necessary. The Laguerre function-based approximation is very similar to transient signals due to its block-oriented feature and thus can be understood easily. The Laguerre functions are expressed in a state-space form, so any state-space control design techniques may be used for a state-feedback controller. However, the G P C algorithm is used for the controller design because of its simplicity and ease of use. The optimal control actions are calculated by minimizing future output errors and control movements. The future outputs are predicted based on a mixed linear-nonlinear Laguerre model under the influence of future control actions. For an adaptive control application, the Laguerre model can be arranged in a linear form in model parameters so that the powerful but simple RLS algorithm can be used for the parameter estimation. The proposed Laguerre model-based control approach, an alternative to a linear transfer function-based control scheme, will be used to handle the severe nonlinearity of a T M P process, i.e., the incremental gain between refiner motor load and plate gap is subject to sudden changes in the sign. The is particularly critical when a refiner is operating close to the maximum load point (see Chapter 6). Chapter 6 Constrained Multivariable Control Of A TMP Plant 6.1 Introduction A thermo-mechanical pulping (TMP) process is a rather complex process. Good control will improve the production, efficiency and pulp quality, while poor control can lead to operating difficulties and poor physical and papermaking properties. The development of a control strategy suitable for the process is a challenging task due to the multivari-able nature, variables' interactions, nonlinearity, time varying, time delay and stochastic disturbances of the process. For the reasons of equipment physical limitations, stable operation and product quality requirement, the plant has to be operated under a number of constraints, e.g. the maximum load attainable by each refiner, the maximum hydraulic closing pressure of a hydraulic loading system, and the desired ranges for the dilution flow and the refining consistency. The controller for this particular process must thus be able to handle the problems efficiently together with the robustness to modelling errors. In the past, various control strategies were proposed towards a control of the process and the product quality. However, many of them were based on the quality-energy rela-tionship, i.e., by adjusting specific energy alone to control pulp quality such as freeness. This is not satisfactory. The earlier research work shows that at least two property vari-ables should be used to predict pulp properties rather than just using freeness alone [66]. Also, a given specific energy may produce fibres with substantially different characteris-tics. Figure 6.40 shows that the given specific energy, Ei = E2, may produce the fibres 144 Chapter 6. Constrained Multivariable Control Of A TMP Plant 145 with substantially different characteristics, where E\ results in fibre cutting and E2 leads to fibre fibrilization. Actually, pulp quality is mainly a function of two variables: the number of refining bar impacts N per unit of pulp and the specific energy per bar impact or refining intensity / . The right combination of these two variables will completely de-termine the pulp quality. Another drawback associated with the past control approaches I Cutting Fibrillation N Figure 6.40: Specific energy is the area under curve of N plotted against I is that, in many cases, the control calculations were made and implemented under the assumption that the resulted control actions were able to be implemented. This is not satisfactory because the inability to implement the control signal exists with the con-cerns of equipment physical limitations, large variations on variables and product quality requirement. Implementing control actions without considering process constraints may lead to degraded process performance and even closed-loop instability. The limitations associated with the past control methods motivated this study towards a better control of a T M P plant. This chapter presents a novel control strategy for the control of a T M P plant by applying the control theories that have been developed and presented in the thesis. The control objective is to achieve specified pulp quality by selecting correct setpoints of the Chapter 6. Constrained Multivariable Control Of A TMP Plant 146 operating variables: specific energy, refining intensity, production rate, motor load and refining consistency. The setpoints of the operating variables are achieved by automat-ically manipulating the closing pressure of a hydraulic loading system, the chip screw speed and the dilution flow rate, without violation of process constraints. The chapter is structured as follows. Section 6.2 describes a choice of variables and the proposed control strategy for the two-stage T M P plant. Section 6.3 presents the control strategy for the control of T M P refiners and simulation results. Laguerre function-based control strategy is developed in Section 6.4 to overcome severe nonlinearity of the process. Summary and conclusions are given in Section 6.5. 6.2 Choice of Variables and Overall Control Strategy 6.2.1 Choice of Variables As mentioned earlier, many pulp properties in terms of papermaking like freeness, fibre length distribution, fibre specific surface, shive content, coarseness, flexibility and so on can be used to describe pulp quality. However, most researchers agree that at least two property variables should be used to define pulp quality. Under this study, the following pulp properties: • Freeness, CSF • Long fiber content, LF • Shive content, SC are chosen to predict pulp quality. The choice of these variables is because (i) these variables can adequately describe the handsheet, strength and drainage properties of mechanical pulp [70], and (ii) they can be reasonably measured on-line at the present Chapter 6. Constrained Multivariable Control Of A TMP Plant 147 time. The targets of pulp properties, normally given by plant operators or from higher level plant optimization systems, are generally expressed using a quality window as: Pulp Properties Lower Bound Desired Value Upper Bound Freeness (ml) Long fibre content(%) Shive content(%) The subscripts min, max refer to the lower and upper bounds in the pulp properties. For a given raw material and refiner design, the pulp properties are strongly affected by the following operating variables at each stage: • Production rate, Fc • Specific energy, E • Refining intensity, J or specific refining power, e • Motor load, P • Refining consistency, C To achieve the setpoints of the above operating variables according to pulp quality re-quirement, the following input variables are manipulated at each stage: • Closing pressure of a hydraulic load system to position the gap, Pc • Chip transfer screw speed to achieve the desired production rate, Sc • Dilution water flow rate to control the refining consistency, Fj, CSFmm < CSF* < CSFt LF^ < LF* < LF, iSCmln ^ SC* ^ SC, Chapter 6. Constrained Multivariable Control Of A TMP Plant 148 6.2.2 Constraints A number of constraints need to be considered during automatic operation of a T M P plant. Amplitude Constraints. The maximum load attainable and the maximum hydraulic pressure are bounded to avoid plate clashes and for the reason of product quality require-ment. The outlet consistency is limited with concerns of both pulp quality requirement and stable operation of a refiner because operation instability may occur in a condition of high refining consistency. The dilution flow rate at each stage is also bounded to avoid valve saturation and a plate clash because plate clashes may happen when the di-lution flow rate is increased up to a certain level [73]. The described possible amplitude constraints are summarized as follows: Amplitude Constraints Variables Lower Bound Typical Value Upper Bound Motor Load Consistency Cmin Closing Pressure Dilution flow Fdmi, Rate Constraints. For the reasons of large variations on the variables and physical limita-tions associated with actuators or instrument, the rate constraints in the closing pressure, dilution flow rate and chip screw speed are considered as follows: Rate Constraints Variables Lower Bound Typical Value Upper Bound Closing Pressure A P c m i n < A P c < Dilution flow A F t i m i n < AFd < Chip Screw Speed A 5 c m ; n < A 5c < P < P •*• _ •* max < c < c„ Pc < PCr, < Fd< Fdn APc™ A F c U AiSCma Chapter 6. Constrained Multivariable Control Of A TMP Plant 149 A controller for this particular process should not be allowed to drive the process variables outside the specified rate and amplitude constraints by taking those constraints into account in the controller design. 6.2.3 Control Strategy Over Two Stages As can be seen from Equation 2.40, pulp quality is expressed as a function of the specific energy E, refining intensity / and specific refining power e. However, specific refining power e can be determined from refining intensity / as given in Equation 2.36. Therefore, Prod. Rate constraints objective J (R:,Sc , F d ) Pulp Quality (CSF,LF,SC) Quality Model El, II i E2,12 N. L. Mech. Model setpoints (Se,P1,C1,P2,C2> Constrained GPCU-T Primary Refiner constraints objective J Constrained G P C U u:Pc ,Fd Secondary Refiner J pulp pulp Figure 6.41: Block diagram for the control of a two-stage T M P plant Chapter 6. Constrained Multivariable Control Of A TMP Plant 150 pulp quality is mainly a function of specific energy E and refining intensity / . The right combination of E and / will completely determine the pulp quality for a given raw material. Figure 6.41 illustrates the proposed control strategy for the two-stage T M P plant i.e., the desired pulp quality is achieved by selecting the setpoints of E and / for each stage. If there were on-line sensors for E and / , there would have automatic control of E and / according to pulp quality requirement. As E is a function of production rate Fc and motor load P (see Equation 2.23), and / is function of refining consistency C for given specific energy and disc rotational speed (see Equations 2.31, 2.32 and 2.33), E and / for each stage are controlled through selecting the setpoints of production rate which is set by chip the screw speed Sc, motor load P and outlet consistency C. The manipulated variables at each stage are the chip screw speed 5 C (for the 1st stage), closing pressure Pc and dilution flow rate Fd- The details regarding the controller design for each stage are given in Section 6.3. To calculate the total specific energy Et and refining intensity It for the two-stage refining, the following equations can be used: Et = E\ + Eo (6.223) It = Et/Nt where Nt = Ni + N2 (6.224) The subscripts 1,2 refer to the first and second stage refining respectively, and the sub-script, t, refers to the total stage refining. N is the number of refiner bar impact received by a unit pulp (Equation 2.33). Chapter 6. Constrained Multivariable Control Of A TMP Plant 151 6.3 Control of T M P Refiners For control purposes, a refiner is considered as a unit operation with three inputs; (i) the wet mass flow rate of chips/pulp, (ii) the hydraulic closing pressure forcing the plates together and (iii) the dilution flow rate to control refining zone consistency. The outputs of a refiner are the specific energy, pulp production rate, motor load, outlet consistency, refining intensity and pulp quality. The process is highly coupled, as most of the inputs simultaneously affect most outputs. The motor load can be easily measured. Consistency measurements become available with the help of some consistency sensor developed lately. Pulp property measurements can be obtained using quality sensors available on the market. In many mills, the production rate is predicted from chip screw speed. Because of the lack of associated sensors, the specific energy and the refining intensity will be predicted from the production rate, motor load and refining consistency based on mechanistic models. The controller for a refiner is developed to achieve the desired specific energy and re-fining intensity through control of the production rate, motor load and outlet consistency by manipulating the closing pressure, chip screw speed and dilution flow rate. 6.3.1 Process Model Any model-based control technique depends upon the existence of a mathematical model to represent the input-output relationship. The dynamic relationships between refiner inputs and outputs about some operating point have been identified as first order plus time delay as follows: P W = i T ipc(t-di)+i T 1Sc(t-d2) + Nl(t) (6.225) 1 — OiiZ 1 1 — 021-Z C 0 = i T iSc(t-ds)+, T 1 F d ( t - d 4 ) + ^ - 1 ) P ( t - l ) + iV 2(t) 1 — O 3 1 Z 1 1 — OA\Z 1 Chapter 6. Constrained Multivariable Control Of A TMP Plant 152 where di axe the time delay associated with each input-output pair. The load disturbance d ( z - 1 ) P ( £ — 1) from the motor load to the consistency is also included to enhance the control performance instead of relying on the load disturbance properties of the controller. Ni(t) represent the joint effect of all unobsesrved disturbances acting within the process and may be modelled by: W) = 1~6£Z~1k(t) (6-226) where & are independent, normally distributed (white) noise sequences with variance tr|.. The dynamic-stochastic model in Equation 6.225 may not be appropriate for on-line parameter estimation, but it can be placed in the form of a C A R I M A model by multiplying both sides by the common denominator polynomial, i.e. A1(z-l)y1,t = Bi + B2(z-l)u2,t_d2 + ^ - \ , « (6-227) A ^ z ' ^ t = £ 3 ( * _ 1 K t - * + £ 4 ( * _ 1 K t - d 4 + £ ( * _ 1 ) i / M - i + ^2,t where Aiiz"1) = 1 - aiXz~x - ai2z~2 Bjiz-1) = bj0 + bjlZ-1 Ci(z-l) = l-ci^-ci2z-*-ci3z-> ( 6 2 2 g ) Diz'1) =do + d1z~1 + ---i = 1,2 j = l , - - - , 4 To simplify the notation in Equation 6.227, the inputs and outputs are re-defined as: yi : the motor load, P y2 : the outlet consistency, C ui : the hydraulic closing pressure, Pc (6.229) u2 : the chip transfer screw speed, Sc u3 : the dilution flow rate, Fd Chapter 6. Constrained Multivariable Control Of A TMP Plant 153 Since many mills use the chip transfer screw to set the desired production rate, the actual screw speed u2 is composed of two terms: the target screw speed u2 associated with the desired production rate and its offset \/u2 manipulated by the controller to compensate for output deviations, i.e. u2,t = V«2 , t + «5,t (6-230) Substituting Equation 6.230 into 6.227 yields: A^z'^y^t = B i ( 2 r - 1 ) M M _ i + B 2 («~ 1 ) v « 2 , t - i + - B 2 ( ^ _ 1 K , t - i + C1(z-l)Cu/A A2{z-l)y2,t = Bz(z'1) v « 2 , t - i + 5 4 ( ^ 1 ) « 3 , t - i + - B 3 ( - 2 " 1 K t - i (6-231) + £>(z - 1 )y M _ i + C 2 ( z - 1 ) 6 , t / A Note that d; — 1 leading zeros have been added to the J3j(z _ 1) polynomial in the above equation to simplify the notation, i.e., deg(-Bj) = n&; + dj — 1. To derive a j-step ahead predictor y^t+j based on Equation 6.231, consider the Dio-phantine identity: diz-1) = E^z-^Mz-1) A +z->FiJ(z-1) (6.232) where E i ) i ( z _ 1 ) (deg(F i j ) = j - 1) and Fj j (z _ 1 ) (deg(F i J ) = max(nai,ncj - j)) are deter-mined by given C^z'1), Ai(z~l) and j. The C j (z _ 1 ) polynomials used for the controller design might differ from the noise terms identified from the process identification stage if the nature of process disturbances or setpoint changes expected in the future is different from that obtained during identification [128]. Multiplying Equation 6.231 by Eitj A z3 and then combining with Equation 6.232 give: Vi,t+j = E^Bi A ult+j^ + EhjB'2 v 4,t+j-i + EUB2 A Chapter 6. Constrained Multivariable Control Of A TMP Plant 154 + F^ylt + Eutw (6.233) V2,t+j = E2jjB'3 v 4,t+j-i + E2jB4 A i4,t+j-i + E2,jB3 A vT{it+i_x + E2JD A ylt+j_x + F2tjyf2tt + E 2 ^ t + j In the above equation, B'^z'1) = A B ^ z - 1 ) effectively eliminates integral action from «2,t- This allows the chip screw speed to be operated about the target production rate. Superscript ' denotes filtering by l/C^z-1). After eliminating future unknown noise terms Eitj£itt+j, the output predictors can be modified in the form: Vi,t+j = Gld A « { i t + i _ i + G2tj v «2,t+j-i + G h A u*{t+j-i + Fijy[ti V2,t+i = Gaj^^j-i + GAjAui^^ + G'^Au'i^^ + Fij^ (6.234) + E2ijD A 2/iit+j_i where = E\jB\ = EhjB2 = E2,3BZ = E2jB4 < ^ = EitjB2 = E2jB3 (6.235) To separate the past known filtered control actions from the current and future un-filtered control actions yet to be determined, consider the identity: G^z-1) = C^iz-^Ciiz-1) + z-'T^z-1) (6.236) where G'ktj(z-1)(deg(Gkj(z-1) = j - l ) and rA. ) j(2:-1)(deg(r fc i i(2:-1) = max(n6fc - 1, nCi -1)) are defined by Gkjiz'1), C^z"1) and j. Substituting Equation 6.236 into 6.234 gives: y M + i = G'ld A u i , t + J - i + G'2J V «2,*+i-i + TU A «i,t-i + r 2 , i V «2 , t - i Chapter 6. Constrained Multivariable Control Of A TMP Plant 155 + G'lj A u%+j_, + F l j y l (6.237) m,t+j = G'ZJ v u2,t+j-i + G'4J A u 3, t+j-i + r 3 i j v " ( t - i + r 4 j A The first two terms on the right side of Equation 6.237 include the effect of the future control inputs yet to be determined: future closing pressure changes (Auitt+j-i), future chip screw speed deviations from the target production rate (vw2,t+j-i) and future dilu-tion flow rate changes (Au3<t+j-i). The rest comprises process free response (open-loop predictions), including the past known controls, measured outputs, load disturbance and future changes in the production rate. Since the G"TJ polynomials in Equation 6.237 are in the order of j — 1 4- nfej, i.e. GUz-1) = 9'io + 9'ixz-1 + ••• + 9h-^-{J~X) + • • • (6 23g) Gltiz-1) = 9lo + 9'^z-1 + ••• + flj-iZ-V-V + ••• all terms except for the jth-term in G'l^z'1) A U g f t + j - i involve the past and future changes in the production rate. The production rate can be assumed constant over the predictive horizon if its future changes are unknown. After eliminating the past and future changes in the production rate, Equation 6.237 is then written as: yltt+j = G'HJ A «i,t+i-i + G'2J v «2,t+j-i + Ti,,- A + T2<J y u{t-i + g'2\j-i A u*f2t + Fldy(t (6.239) m,t+j = G'3J v «2,t+j-i + G\4 A M 3 ) t+ i - i + r 3 ) j v « 2 ) t - i + r 4 j A u{t_x + g'l^ A u% + E2JD A yf1<t+j + F2Jyf2ii The above equation can be modified in a more compact form over the predictive and control horizons: yi = Gi*i i + G ' 2 u 2 + fi ( 6 2 4 o ) y2 = G3U2 + G^u3 + f2 Chapter 6. Constrained Multivariable Control Of A TMP Plant 156 where and fi = y i = bi,t+i yi,t+2 • • • 2/i,t+iv2]T y 2 = [y2,t+i V2,t+2 • • • y2,t+Ni\T u i = [A«i ) t A ui,t+i • • • A uitt+NU-i]T U 2 = [Vw2,t V "2,t+l • • • V u2,t+NU-l}T u 3 = [Au 3 , t A w3)t+i • • • A •u 3 i t + j V r 7 - i ] : r r u A «{ft_! + r2,i v u{t-i + FhivL + 0 2 , 0 A U*L Tl.2 A u{ i t _ x + r2j2 V « 2 , t - l + * l , 2 V l , t + 02,1 A U*2,t rliJVa A «{it_! + r2)JV2 v « 2 , t - i + ^i,jv2j/i,t + tftM-iA (6.241) (6.242) f2 = r3,i V « 2 , t - i + r 4 , i A u( t _! + F 2 ) 1 y ( t + p£0 A u*f2,t r3,2 V « 2 , t - i + T4 , 2 A «J)t_! + F 2 | 2 y ( t + A u*£ t r3,Af2 v ui,t-i + r4,iv2 A u( t _! + F2 , iv 2y 2 , t + s ^ - i A u * l t J (6.243) 6.3.2 Controller Development Simply stated, the control objective is to maintain the production rate, the motor load and the outlet consistency at their respective setpoints according to the desired specific energy and refining intensity through manipulating the chip screw speed, the closing pressure and the dilution flow rate. Since the number of inputs, in this application, is one greater than the number of outputs, one of the inputs must be constrained about some steady-state value [128]. Many mills use the chip screw to set the desired production rate and, therefore the controller should constrain the chip screw speed about the target production rate while still leaving some room for manipulation. The described control Chapter 6. Constrained Multivariable Control Of A TMP Plant 157 strategy can be achieved by designing a controller to minimize the following quadratic cost function: f N a min J = E l J2 \Qi(Vi,t+j - Wi,t+j? + "2(2/2,*+; - w2,t+j)2 A-ui,v«2,Au 3 I j = N l L NU I + £ [pi A u h + j - i + P2 V 4 t + i - i + Pa A ult+j_x] J (6.244) where Oj are the output weighting factors, pi are the control weight factors and Wi are the output setpoints. Sju2 is the difference between the actual chip meter speed u2 and the target production rate u\ (u2%t = w2 + V u 2 ) - ' • The above quadratic cost function can be modified and put in a more compact form: min J = 4i(yi - w i ) T ( y i - wi ) + q2(y2 - w 2 ) T ( y 2 - w a ) U l U2 Us + p i U j i i i + p2U2 u 2 + P3uju3 (6.245) where W{ are the vectors of the future setpoints, i.e., = [witt+i Wijt+2 • • • wi<t+N2]T• Substituting Equation 6.240 into the above equation gives: m i n U l a 3 f i , J = ft[Giui + G ' 2 u 2 + f i - w 1 ] T [ G ' 1 u 1 + G 2 u 2 + f 1 - w x ]+ +9 2 [G^u 3 + G ^ u 3 + f2 - w 2 ] T [ G 3 u 2 + G^iis + f2 - w 2 ] (6.246) +p1uju1 + p 2 u 2 U 2 + /53U3U3 As stated earlier, an industrial refiner needs to be operated under a number of con-straints. The necessary constraints given earlier are rewritten as: | A u l i t + i _ i | < At t i Wlmin < Ml,t+t-l < W i m a x 2/lmin < Vl,t+j < J/lmax |At*2,t+i-l| < A « 2 «3min < M3,t+i-l < «3max 2/2 min < 2/2,i+j < 2/2 max . n A n . | A « 3 , t + i - i | < A « 3 (6-247) Z = 1 , - - - , JV[ / Chapter 6. Constrained Multivariable Control Of A TMP Plant 158 If none of the constraints is violated, the control actions u* = [uj uJ]T derived by minimizing Equation 6.246 are optimal. If, however, one or more constraints are violated, u* will no longer be optimal. The optimal solution should be derived by minimizing the DESIGN constraints objective J P setpoint -C setpoint -i i i •* Constrained „t 1 |CP G P C +A 1 TS^ Output Prediction Model Update DW PRIMARY REFINER PR setpoint P C P = Motor Load C = Outlet Consistency PR = Production Rate CP = Hydraulic Closing Pressure TS = Transfer Screw Speed DW = Dilution Water Flowrate Figure 6.42: Constrained multivariable adaptive-predictive control of the primary refiner quadratic cost function in Equation 6.246 with the satisfaction of the constraints in Equation 6.247. The following is the derivation of the optimal solution through solv-ing a quadratic programming (QP) problem and a mixed-weight least-squares (MWLS) problem. Optimal solution via analytical QP The cost function in Equation 6.246 subject to the constraints in Equation 6.247 can be modified and put in the form of: minfc J = u Hu + 2 c T i i subject to: (6.248) A(u) > 0 Chapter 6. Constrained Multivariable Control Of A TMP Plant 159 where u = [uJ uJ u3]T (6.249) hn hn. 0 H = h2\ h22 h2z 0 h32 q1G\TG>1+p1 and lTn,T 2 ^ 1 g i G 0 Cl c = c 2 = c 3 g i G ' / G ' / 0 l 2T G 2 + g 2 G 3 T G 3 + P2 g 2 G ^ T G 3 0 o2 G 3 G 4 g 2 G ^ T G ^ + p 3 (6.250) (6.251) ftGi^fx-W!) q1G'2T(f1 - w x ) + g 2 G 3 T ( f 2 - w 2 ) 9 2 G ^ T ( f 2 - w 2 ) Matrix A(u) = [ai a 2 • • • ae]T is a linear combination of the constraints and the control actions yet to be determined. Generally, the quadratic cost function in Equation 6.244 solves both for the current input and NU — 1 future inputs. However, here we use NU = 1 for all three inputs since this tends to simplify the computations. For the case of NU = 1, the constraints in Equation 6.247 can be modified in the form (see Appendix B.I.): ai = a2 = a 3 = a 4 = a 5 = a 6 = A«i , t A u i , t V«2 , t V « 2 , t A u 3 | t A i t 3 , t -cti > 0 A > 0 -a 2 > 0 ft > 0 - a 3 > 0 & > 0 (6.252) Chapter 6. Constrained Multivariable Control Of A TMP Plant 160 where cti = m a x ( - A H 1 , t i i m i n - t i i , t - i , [1 0 0] • P • hL) p\ = min( A u n « l m a x - « i , t - i , [1 0 0] • P • hH) ct2 = max( -Au 2 - Au*2<t + V « 2 , t - i , [0 1 0] • P • hL) p\ = min( A u 2 - A u ^ + V«2,*-i> [0 1 0] • P • hH) a 3 = m a x ( - A u 3 , u3mxn - w3)<_i, [0 0 1] • P • hL) p\ = min( A u 3 , U 3 max - « 3 , t - i , [0 0 1] • P • /in) (6.253) and P = G' x G' 2 0 0 G(, G^ G i G' a 0 0 G 8 G 4 -1 G i G' 2 0 0 G'« G^ T 3^ 1 m i n - f l ' y i m a x - f l " y2 m i n - f 2 . y 2 m a x - f 2 . (6.254) (6.255) The constrained control problem in Equation 6.248 can be further expressed as an un-constrained control problem by including the Kuhn-Tucker multipliers as: min J = uTHu + 2cTu + pTA (6.256) u where A is a vector of defined constraints: A = [di a2 a 3 0,4 a 5 a6]T (6.257) and p is a vector of Kuhn-Tucker multiplier associated with the constraints: P = [Pi Pi P3 Pi P5 P&]T (6.258) Substituting Equations 6.250, 6.251, 6.257 and 6.258 into Equation 6.256 yields: mina J = h n A w 2 t + (h21 + J112) A u 1 | t y «2,t + h22 V ul,t + (h32 + h23) v "2,t A u 3 > t + h33 A u\t + 2ci A u u + 2c2 v U2,t + 2C3 A u3<t + /i iOi + / i 2 a 2 + / i 3 Q 3 + P4CL4 + P5&5 + P60-6 (6.259) Chapter 6. Constrained Multivariable Control Of A TMP Plant 161 The Kuhn-Tucker multipliers method requires that for u to be optimal, it should satisfy the following conditions: 1. The first order conditions are: dJ/du = 0 (6.260) 2. The second order conditions are: aj/ij = 0 for all i = 1 to 6 such that a* = 0 if m<0 (6.261) or a, > 0 if /ij = 0 The optimal solution u of Equation 6.256 will be obtained when the 1st and 2nd order conditions are satisfied. The steps involved for the calculation of u are given below: Step 1: Determine an unconstrained control solution u*. Step 2: Test to determine if any of the constraints in Equation 6.252 are violated. If all aj > 0 (no constraints are violated) then go to Step 5. Step 3: Saturate the solution on one or two constraints and verify if the solution satisfies the above 1st and 2nd order conditions: • Assume any one or two (say am and an) of the constraints are satisfied by the constrained optimal solution, which gives a value of i i . Therefore, all ^'s other than u.m and (j,n are zero. • Calculate the Kuhn-Tucker multipliers /x m and //„ with the help of Equa-tion 6.260. Verify if the Kuhn-Tucker multipliers u.m, /xn < 0 Chapter 6. Constrained Multivariable Control Of A TMP Plant 162 If the Kuhn-Tucker multipliers are negative, then go to Step5. Step 4: Repeat Step 3 for the next (set of) constraint(s). Step 5: Implement the control action and repeat Steps 1-5 at next sample interval. A total of 18 sets of constraints are possible in Step 3 of the above algorithm. They are listed in Appendix B.2. Optimal solution via the mixed-weights least-squares algorithm The optimal solution for the general case of constrained control problems can also be obtained through solving a mixed-weights least-squares problem. To do so, the quadratic cost function in Equation 6.245 is re-arranged as a 2-norm minimization problem as: Y i - w i y 2 - w 2 PIINU ®NU ®NU \ ®NU ®NU PSINU } where u = [uJ uJ u3]T. Substituting Equation 6.240 into the above equation and mod-ifying yields: 11 (6-263) min J u \ 1/2 U (6.262) minJ= ll itf i-VI 1 2 where R = G i 0 PIINU ONU L \ °NU G' 2 G 3 P2^NU 0 G 4 PZINU J 1/2 (6.264) and Chapter 6. Constrained Multivariable Control Of A TMP Plant 163 W i - f i y = w 2 — f 2 0jv 2xi The input amplitude constraints in Equation 6.247 can be modified as: |L„u - Cu < 1 (6.265) (6.266) or expressed as an oo-norm problem as: Luii - Cu lloo^ 1 (6.267) where Lu — and .radius A/u3 ,radius Cu — {Ultcentre ~ ^ l , t - l ) - ' N U x l / / u l , r a d i u s OjVuxl .centre ~ u 3 , t - l ) l N U x l / u 3 , radius "lrnax + « l i Ml radius max - U i min (6.268) (6.269) _ ^3max + ^3min _ M3max ~ ^3min / „ 0 7 n \ u3,centre — ^ u3,radius — ^ (0.2/UJ Matrix A is an NU x NU lower diagonal matrix whose non-zero entries are all equal to 1. The output amplitude constraints in Equation 6.247 are modified in the form of: Lyu — Cy\ < 1 (6.271) or expressed as an oo-norm problem as: LyXL — Cy | 1 (6.272) Chapter 6. Constrained Multivariable Control Of A TMP Plant 164 where and L„ — G'l/y\,radius Gr^/ 'y i , radius ^N2xNU ON2XNU .radius *-*4/V2,radius Cy — radius (y2,centrelN2xl ~ ^^2 ,radius yi,centre 2/l max + 2/1 min yi,radius 2/l max 2/l min (6.273) (6.274) 2/2 max + 2/2 min 2/2 max 2/3 min y2,centre — 2~ 2/2,radiua ~~ 2~ T/ie m/?u£ rate constraints in Equation 6.247 are re-arranged in the form: |LduH ~ Cdu| 1 or expressed as an oo-norm problem as: where r INU/Altj 0NU 0NU 'du= I 0jv(7 ^ 4 I / A M 2 Ojvry Oivc; Ojvrj / w / / A u 3 and 1 0 0 - 1 1 0 ••• 0 0 - 1 1 0 0 ••• - 1 1 A2 = ®NUxl Cdu — A2/Au2 OjVUxl (V«2,t-1 - A«5, t)/A1L 0 0 (6.275) (6.276) (6.277) (6.278) (6.279) Chapter 6. Constrained Multivariable Control Of A TMP Plant 165 Combining Equations 6.267, 6.272 and 6.277 gives: I I L u - C H o o ^ l where (6.280) Cu L = Ly C = Cy Ldu Cdu (6.281) The constrained control problem is then solved by minimizing the cost function in Equa-tion 6.263 with the satisfaction of the constraints in Equation 6.280, i.e. m i n f i J = | | i 2 u - V\\l subject to: (6.282) I  L u - C H o o ^ l The optimal solution of the above constrained problem can be obtained through mini-mizing the following mixed-weights least-squares or mixed-objective function: 1/2 mm W W R L u — V C (6.283) where is a positive real scalar and is a diagonal matrix with positive entries. Both and are calculated iteratively as: W Z^fc=l wkk ek\ W V V 33 CJ 33 Z^fe=l wkk ek where and e W = L u ( i ) - C (6.284) (6.285) (6.286) WJf are the entries of the and m is the number of the constraints (m < 2NU + N2). Chapter 6. Constrained Multivariable Control Of A TMP Plant 166 6.3.3 Simulation Results and Discussions The proposed control strategy described previously was tested on a combination of empir-ical and mechanistic models (see Figure 6.43), to approximate the features of an industrial situation. The dynamics from inputs to outputs was based on the identification results performed on industrial refiners. The static relationship between inputs and outputs was simulated based on the mechanistic models. The effect of operating variables on refining disturbances Pc '—*\ 1 st order + delay N.L Mech. Model | i ... St g*: 1 st order + delay Fd j—*1 1 st order + delay - • F c - * -P c the primary refiner Pc: closing pressure N PQ E : specific energy Fc: production rate P : motor load St: chip transfer screw speed Fd: dilution flowrate c : o u t l e t consistency x : residence time I : refining intensity N : number of bar impacts PQ: pulp quality Figure 6.43: Input-output dynamic and static relationship for the primary refiner consistency was modelled using mass and energy balances. Production rate control was operated at steady state and it was reasonably assumed to be a linear function of the transfer screw speed. The model used for the controller design was given in Equation 6.240. In the following simulations, a nominal set of tuning parameters qi = q2 = 1, iVi — 1, N2 = 15, NU = 1 and Ai = A 2 = 0.98 were used. * Response to Rate Constraints Chapter 6. Constrained Multivariable Control Of A TMP Plant 167 Due to equipment physical limits, manipulated variables should be bounded before their implementation. Figure 6.44 shows constrained control via M W L S with the rate constraints on manipulated variables. The hard limits on the closing pressure, screw speed and dilution flow were set to | A i t l | < 50, | Au2\ < 0.2 and | Aw3| < 0.5| respectively. The controller parameters were p2 = 1.0, P\ = p 3 = 0.0001. No constraints on the amplitudes of the manipulated variables were implemented. To compare constrained control with unconstrained one, the simulation was run without hard constraints for the first 200 intervals and constrained control was turned on at interval 200. Figure 6.44 shows the movements of the manipulated variables were effectively reduced into their desired limits. Response to Amplitude Constraints Amplitudes of manipulated variables should also be limited to their desirable ranges for the reasons of safety, stable operation and equipment physical limitation. This is particularly important in a T M P process, e.g. maximum closing pressure applied for each refiner to avoid plate clashes, maximum and minimum dilution flow for the reasons of stable operation and dilution valve's physical limit. Figure 6.45 shows the simulation performed without constraints for the first 200 sec. and with hard constraints for the rest of time. The amplitude constraints on the closing pressure « i and dilution flow M 3 were set to 820 < M i < 1150 and 3.4 < M 3 < 5.2. Figure 6.45 demonstrates that the setpoints were achieved by manipulating closing pressure and dilution flow within their desirable ranges while more room for manipulation on the screw speed to compromise the outputs deviations from the targets. Response to Deterministic Disturbances and Changes in Process Dynamics Figure 6.46 shows the response to a number of deterministic load disturbances and changes in the process dynamics. Six time series plots are shown in the figure: (i) the motor load (yl) and its setpoint change, (ii) the consistency (y2) and its setpoint, (iii) Chapter 6. Constrained Multivariable Control Of A TMP Plant 168 CM 15 10 I 50 40 1 1 V : \f IV r \ 1 1 0 1500 ^ 1000 500 i 30.1 ^ 30 29.9 co " = 4 2 I 100 0 -100 I 0.2 -0.2 50 — i — 50 —\— 50 50 50 50 50 100 — i — 150 — i — 200 250 — i — 300 — i — 100 — I — 150 — i — 200 250 300 100 150 200 250 300 100 150 200 250 300 100 150 200 250 300 100 150 200 250 300 100 150 200 250 300 SV-~!— 50 100 150 200 250 Time (sec.) 350 — i — 350 350 350 350 350 350 — i — 300 350 400 400 400 A I A \ I*-Y f V 400 400 1 | | 1 ; : A n r ; ^ I: I i i i 1 i i i i 400 I { r T r p i—"t_>— V 400 400 Figure 6.44: Constrained control via M W L S , t> 200 sec: | A « 3 | <0.5 (ql=qa = l,p2 = 1.0, P l = p3 = 0.0001) A u i | <50, | A t*21 <0.2, Chapter 6. Constrained Multivariable Control Of A TMP Plant 169 y1: motor load (MW) 16 14 12 10 0 100 200 300 400 u1: closing pressure (psi) 1500 1000 0 100 200 300 400 u2: screw speed (rpm) 30.1 30 29.9 \ — r 0 100 200 300 400 Time (sec.) y2: consistency (%) 0 100 200 300 400 u3: dilution flow (kg/sec) i 1 — r " • ? 1 r~—I 1 \ r "0 100 200 300 400 production rate (t/d) 326 325 324 r 0 100 200 300 400 Time (sec.) Figure 6.45: Constrained control via analytical QP, t>200 sec: 820 < Ui < 1150, 3.4 <u3< 5.2 (qi = q2 = l,p2 = 1.0, px = p3 = 0.001) Chapter 6. Constrained Multivariable Control Of A TMP Plant 170 the closing pressure (wl), (iv) the target and actual screw speed (u2), (vi) the dilution flow (u3), and (vii) the target and actual production rate. A set of tuning parameters P2 = 10, pi = P3 — 0.0001 were used with the objective of maintaining the production rate to its target while manipulating the closing pressure and dilution to reach the setpoints. Load disturbances both from the motor load to the consistency and the dilution to the motor load were simulated in Figure 6.46. A n increase (or decrease) in the motor load results in an increase (or decrease) in the consistency. This is because increased (or decreased) motor load results in more (or less) water converted into steam and therefore, an increase (or decrease) in the consistency. To compensate for the motor load effect on the consistency, the dilution flow was adjusted to maintain the target of the consistency. Due to process nonlinearity, an adaptive scheme was employed to estimate time-varying gain on-line. Figure 6.46 shows that proposed controller performed well with process nonlinearity and time-varying gain. Tuning to Reduce Chip Screw Manipulation Figure 6.47 shows a simulation performed with the objective of decreasing the amount of chip screw manipulation while maintaining the targets of the motor load and consis-tency. This may be considered a realistic control objective because excessive screw speed manipulation can result in a large variation in production rate and thereby, varying the specific energy applied and pulp quality. To reduce screw manipulation, p 2 was set to p2 = 10. Figure 6.47 shows substantial reduction in the screw speed by increasing associated weighting factor, without sacrificing output performance. The pulp quality response to setpoint changes Figures 6.48 and 6.49 show a set of simulation runs undertaken with the changes in the motor load setpoint. The initial motor load setpoint was determined from given total Chapter 6. Constrained Multivariable Control Of A TMP Plant y1: motor load (MW) 100 200 300 400 u1: closing pressure (psi) 0 100 200 300 400 u2: screw s p e e d (rpm) r 0 100 200 300 400 T i m e (sec.) y2: consistency (%) 50 45 40 35 i 8 6 4 2 .IL Tr-i m 0 100 200 300 400 u3: dilution flow (kg/sec) 1 r 1 i h h > n i If : i U \ 0 100 200 300 400 production rate (t/d) 325.5 325 324.5 \ 1 T f 1 i : [ I -1 0 100 200 300 400 T i m e (sec.) Figure 6.46: Unconstrained controller (qx = q2 = l,p2 = 1.0,/?i = pz = 0.0001) Chapter 6. Constrained Multivariable Control Of A TMP Plant y1: motor load (MW) y2: consistency (%) 16 14 12 10 i 1500 1000 500 i 30.01 0 100 200 300 400 u1: closing pressure (psi) 0 100 200 300 400 u2: screw s p e e d (rpm) 30 29.99 r-f\ 0 100 200 300 400 T i m e (sec.) 0 100 200 300 400 u3: dilution flow (kg/sec) 8 6 4 2 i 325.5 325 324.5 r f n / |f -. . . . (... .. H ; If : • i u f "0 100 200 300 400 production rate (t/d) Tt 0 100 200 300 400 T i m e (sec.) Figure 6.47: Unconstrained controller (cji = q2 = l,p2 = 10, pi = pz = 0.0001) Chapter 6. Constrained Multivariable Control Of A TMP Plant 173 specific energy, power split and production rate as follows: Total specific energy = 6.0(MJ/kg) Power split = 57/43 Production rate = 325(t/d) The consistency setpoint was set to 45%. These setpoints are given according to the requirement of product quality. Perturbation signals were added at each control output to obtain persistent excitation due to an adaptive scheme applied. The tuning parameters were set to qx — q2 = l,p2 = l , P i = Pz = 0.0001. A n increase in the motor load setpoint was made at t = 150sec. and a decrease in the consistency setpoint was made at t = 300sec. Figure 6.48 shows that the motor load setpoint was achieved through manipulating the closing pressure while the screw speed was adjusted around desired production rate. In practice, an increase in the motor load will increase the consistency. Figure 6.48 shows that consistency variations resulting from the load setpoint change was compromised through automatically increasing the dilution flow. Figure 6.49 shows that at the given production rate, an increase in the motor load increased the specific energy. As a result, the inlet consistency decreased for a given out-let consistency because with higher specific energy, there will be more water converted to steam and therefore, less water in the refining zoom. To keep the consistency in the refining zone constant, inlet dilution was increased and as a result the inlet consistency reduced. The residency time is proportional to the specific energy and the inlet con-sistency. Increased residence time indicated that an increase in the specific energy was relatively significant than a decrease in the inlet consistency. Increases in the refining intensity and the specific power were due to a relatively significant increase in the specific energy. A significant decrease in the freeness was resulted from a decrease in the specific energy. Increased refining intensity and specific power due to an increase in the specific energy resulted in decreases in the long fibre and shive contents. Chapter 6. Constrained Multivariable Control Of A TMP Plant 174 Figure 6.48 shows that a decrease in the consistency setpoint was made at t = 300sec. Wi th the proposed control, the consistency setpoint was achieved by increasing the dilu-tion flow and slightly decreasing the chip feed due to a tight constraint on the production rate. Figure 6.49 shows that for a given specific energy, decreased consistency resulted in decreases in the inlet consistency and the residence time whereas increases in the refining intensity and the specific power. The explanation is as follows. Increased dilution flow at a given specific energy brought about a higher centrifugal force on the pulp flow. As a result, the residence time for a pulp to pass through the refining zone was reduced and thereby, the total number of bar impacts received by the pulp was reduced. At a given specific energy, reduced number of bar impacts increased the average specific energy per impact, i.e., the refining intensity. Increased refining intensity and specific power caused more fibre cutting, resulting in less long fibre content and shive content. Due to a strong relationship with the specific energy, the freeness was almost unchanged due to a constant specific energy. Chapter 6. Constrained Multivariable Control Of A TMP Plant 175 u2: screw s p e e d (rpm) production rate (t/d) 0 100 200 300 400 T ime (sec.) 0 100 200 300 400 T ime (sec.) Figure 6.48: Unconstrained controller (p2 = l , P i = Pz = 0.0001) Chapter 6. Constrained Multivariable Control Of A TMP Plant specific energy (MJ/kg) 3.5 O 100 200 300 400 residence time (sec.) 0.35 0.3 0.25 O 100 200 300 400 specif ic power (MJ/kg/sec.) O 100 200 300 400 30 25 20 inlet consistency (%) 6 5.5 5 4.5 550 500 450 l*A*Mi 0 100 200 300 400 x - | o~ 4 r e f i n i n 9 i n t e n s i t y (MJ/kg) i AJU it. An 0 100 200 300 400 C S F (ml) W i /I O 100 200 300 400 55 50 45 40 L F (%) 100 200 300 T ime (sec.) 400 S C (%) ' ^ y ^ W v ^ ^ 100 200 300 T i m e (sec.) 400 Figure 6.49: Unconstrained controller (p2 — 1, p\ = p3 — 0.0001) Chapter 6. Constrained Multivariable Control Of A TMP Plant 177 Extreme Point Plate G a p (x 0.1 mm) Figure 6.50: The gain between motor load and plate gap is subject to slow drift and a sudden change in the sign from negative to positive 6.4 Laguerre Function-Based Control of TMP Refiners Control of a wood chip refiner in TMP manufacturing is sometimes difficult due to the incremental gain between refiner motor load and plate gap subject to a sudden change in the sign. This is particularly critical when the refiner is operating close to the maximum load point (see Figure 6.50). The linear approximation-based control for this case may not be satisfactory and may lead to unstable operation. Alternatively, the Laguerre functions presentation is used in the following to approximate process dynamics and static nonlinearity. A generalized predictive control algorithm is then derived based on the nonlinear Laguerre function representation. 6.4.1 Nonlinearity Description It has been known for a long time that the incremental gain (see Figure 6.50) between refiner motor load and plate gap is subject to a slow drift due to plate wear and to Chapter 6. Constrained Multivariable Control Of A TMP Plant 178 sudden changes in sign due to collapse of the pulp pad in the refining zone. Plate wear generally occurs gradually over a period of several hundred hours, but can also occur suddenly in the case of plate clash, i.e., metal to metal contact between two plates. The maximum load corresponds to the minimum gap below which the pulp pad in the refining zone cannot sustain the pressure and then collapses, resulting in plate clashes and the incremental gain between the load and the gap changing in sign from negative to positive. A plate clash is highly undesirable because it will disrupt production and damage the plates. To avoid a plate clash, the gap has to be opened beyond the point where the collapse occurred, i.e., the refiner has to be operating in the region where closing the gap increases the load. However, the gap at which the pad collapses is unpredictable, and it depends on wood species, wood chip quality, refining consistency, plate condition, etc. These make refiner control more difficult. 6.4.2 Past Control Solutions A first attempt at applying adaptive control to the problem was made by Dumont [30]. In his solution, Fortescue's variable forgetting factor based recursive least squares esti-mator was used to identify the gain, and a set of rules were set to back out the refiner plates and restore the pulp pad when the gain changed sign. Some success on an in-dustrial chip refiner was reported, but the choice of the variable forgetting factor design parameters to achieve the required trade-off between tracking slow and fast gain changes proved to be difficult in practice [33]. Dumont and Astrom [32] then investigated several alternatives including dual control and proposed use of a nonlinear performance index which explicitly penalizes operation in the pad collapse region. In their proposed control strategy, additional terms reflecting the nonlinear nature of the process were introduced into a performance index, instead of using an index function reflecting the output error alone. Using this control criterion, penalty is given when the estimated gain is small. In Chapter 6. Constrained Multivariable Control Of A TMP Plant 179 the case of an unreachable motor load setpoint, the controller will tend to keep the gain negative instead of eliminating the output error. When the motor load is lower than the setpoint, such control will lead the control action to close the gap when the estimated gain is negative with relatively large amplitude and open the gap when gain tends to be positive. Their idea was further developed and tested on an industrial refiner by Allison et a/. [48]. Preliminary industrial results showed the general success and the potential of industrial control implementation. However, one drawback is the computational diffi-culty attached to solving the dual control problem on-line. A new way of representing a refiner was proposed by by Dumont and Fu [114] for the control of a wood chip refiner. In their control scheme, the Laguerre function expansion via Volterra kernel was used to model process nonlinearity and dynamics, and the clashes of the refiner plate was avoided without adding extra measurements and computational difficulty which are associated with dual control problem. The control approaches given above were based on single-variable control philosophy (i.e., the pulp quality is determined by the energy input), so they were limited by their control performance. This is because pulp quality is mainly a function of two variables: the specific energy and the specific energy per bar impact or refining intensity. For a given plate configuration and a constant disc rotational speed, the refining intensity is, with some approximation, largely dependent on refining consistency. To obtain desired pulp quality, the specific energy as well as the refining consistency need to be controlled. The following discussion presents a multivariable Laguerre function-based control scheme in which the specific energy and the refining intensity are controlled through control of the motor load and refining consistency by manipulating the hydraulic closing pressure to position the gap and the dilution flow rate. Chapter 6. Constrained Multivariable Control Of A TMP Plant 180 6.4.3 2 x 2 Mixed Linear-Nonlinear Laguerre Function-Based Control Of T M P Refiners A 2 x 2 Laguerre function-based predictive control scheme developed in Section 5.4 is applied towards stable control of a T M P refiner when the refiner is operating close to the maximum load. 1st order terms 2nd order terms 1st order terms In l - a , z z—a, D...1...-1 [l2 1 1 ••• l 2 1 [>3.. l - a 2 ; z - a , Figure 6.51: Structure of a 2x2 Laguerre model for a T M P refiner Laguerre Model Figure 6.51 shows the 2 x 2 Laguerre model structure using the first N filters for modelling a wood chip refiner. The nonlinear dynamic relationship between the motor load and plate gap is modelled using the nonlinear Laguerre functions representation, whereas the dynamic relationship between the outlet consistency and dilution flow rate is approximated by the linear Laguerre functions representation. The mathematical formula of this Laguerre model structure was given in Equation 5.205. Adaptive-Predictive Control based on Laguerre Model Chapter 6. Constrained Multivariable Control Of A TMP Plant 181 A generalized predictive control algorithm is then derived based on the mixed linear-nonlinear Laguerre model. The control objective is to maintain the motor load y i and refining consistency y2 to their perspective setpoints according to the desired specific en-ergy and refining intensity by adjusting the plate gap through manipulating the hydraulic closing pressure Ui, and the dilution flow rate u2 to the refiner. The control actions are calculated by minimizing the performance index given in Equation 5.213. The Laguerre model given in Equations 5.205 is used to predict process output trajectories. The dif-ferences between the predictions and actual process measurements are used to estimate the model parameters on-line to handle the time-varying nature of the process. 6.4.4 Simulation Results and Discussion The described adaptive-predictive control scheme based on the mixed linear-nonlinear Laguerre model has been tested extensively in simulations. The simulation results show its potential applicability and illustrate its practicality as follows. The simulations were performed on a combination of empirical and mechanistic models given in Equations 6.287 ~ 6.290, to approximate the features of an industrial situation. The dynamics from inputs to outputs were based on the identification results performed on industrial refiners. The static relationship between the plate gap and the motor load was simulated as shown in Figure 6.50. The effect of operating variables on the consistency was simulated using the results based on the mass and energy balances. The above described simulation models may be presented as: X i t t = anXltt-x + & n U i , t - 2 + & 1 2 « l , t - 3 (6.287) X2,t = a 2 i X 2 , t - l + &21«2,t-2 + &21«2,t-3 (6.288) (6.289) Chapter 6. Constrained Multivariable Control Of A TMP Plant 182 + N2,t (6.290) V2,t = + x2,tL - rj y1)t where y\ is the motor load, y2 is the outlet consistency, U i is the pulse sent to the hydraulic cylinder to position the gap x\, and u2 is the dilution flow rate used to control the consistency. iVj in the above equation account for the joint effect of all external disturbances and measurement noise on an output. In the following simulations, the linear dynamic parameters in the above equations are: an = 0.75, bn = 0.25,612 = 0, a 2 i = 0.8187,62i = 0.1813,622 = 0. The nonlinear relationship given in Equation 6.289 is simulated with gi = 7600, g2 = 800, a = 1 in a normal state, and a = 0.77 in a collapsed state when x < 0.9. Other operating variables are: the pulp production rate Fc = 3.7616 kg/sec, the latency heat L = 2.5 MJ /kg , the inlet solid content s — 0.5 and the power efficiency coefficient n = 0.9. However, those parameters and nonlinear functions in the above equations are supposed unknown. There are white noises added to the motor load and consistency outputs, simulating measurement noises. The predictive controller tuning parameters, in the following simulations, are set to Ni = 1, N2 — 15, pi = 0.0 and p2 = 0.0001. The amplitude and rate constraints on the plate gap are applied, which are 0.7 < Xij < 5.0 and | A xltt\ < 0.1. To start with simple Laguerre functions representation, two Laguerre filters are used with the Laguerre filter poles ax = a2 = 0.7. The severe nonlinear relationship between the plate gap and the motor load is modelled by the 2nd order nonlinear Laguerre function. The RLS estimation algorithm is used with the forgetting factors: Ai = A 2 = 0.98. Motor Load Setpoint Unreachable. Figures 6.52 shows a simulation of the proposed controller when the load setpoint is unreachable. For the period up to about t = 80 sec, the refiner is operating in the normal state, i.e., closing the gap increases the motor load. The controller for this case is able to adjust the plate gap through manipulating the closing pressure to minimize the error between the load and its setpoint. The load Chapter 6. Constrained Multivariable Control Of A TMP Plant 183 setpoint then begins to increase at t = 80 sec from 6380 k W to 6580 kW. The controller, in the beginning, closes the plate gap in order to reach the load setpoint. As the gap decreases, the extreme point (see Figure 6.50) will be reached, because of the setpoint is unreachable in this case. The controller based on the estimated nonlinear Laguerre model is able to identify the extreme point and evaluate the control input close to the true extreme point 0.9. Correspondingly, the motor load output is then stabilized below its setpoint. Plate Clash Avoided. Figure 6.53 demonstrates how a plate clash can be avoided by the proposed control strategy. As can be seen at about t = 350 sec in Figure 6.53, when a refiner is operating near the maximum load point, just process disturbances can produce a pulp pad collapse. As a result, the refiner will be operating on the left side of the extreme point (see Figure 6.50), i.e., closing the gap decreases the motor load. For this case, the nonlinear Laguerre model based controller is able to detect this collapse and quickly open the gap, allowing the pulp pad to rebuild. The refiner then goes back to the normal operation state and a plate clash is avoided. Changes in Consistency Setpoint. Figure 6.54 shows a simulation of the proposed control scheme when the consistency setpoint changes. As can be seen, a step change in the consistency setpoint is made at t = 220 sec and the proposed controller is able to adjust the dilution flow rate to meet the setpoint. Although the effect of operating variables on the consistency is nonlinear (Equation 6.290), the proposed controller based on the estimated linear Laguerre approximation is able to identify the nonlinearity as a gain change and evaluates the control input accordingly. Response to Load Disturbances. In practice, there is a coupling between the motor load and refining consistency, i.e., an increase in the motor load will increase the consistency and vice versa. Figure 6.52 shows that to maintain the outlet consistency constant when the motor load increases resulting from an increase in its setpoint, the dilution flow rate is Chapter 6. Constrained Multivariable Control Of A TMP Plant 184 increased. A decrease in the motor load will decrease the outlet consistency. To maintain a constant consistency, the proposed controller automatically reduces the dilution flow to the refiner. Figure 6.54 shows an increase in the outlet consistency will result in an increase in the load. To keep the load constant, some control effort goes into decreasing the load which is done by opening the plate gap. Simulations of the proposed control scheme demonstrated that it is more robust than the conventional model-based approaches in the sense of no need for assumption of the model order and time-delay, and is able to handle an industrial process with severe nonlin-earity, unknown dynamics and external disturbances. In the case of a pulp pad collapse, the proposed controller was able to identify the collapse and quickly open the gap for the pad to rebuild, whereas the linear model-based control schemes for this case may cause some control problems such as process oscillation, pulp pad collapse acceleration and plate clashes. Using Laguerre functions approximation, plant nonlinearity and dynamics can be modelled properly with a flexible model structure and less model parameters. 6.5 Summary and Conclusions A new control strategy was developed towards better control of a two-stage T M P plant, where pulp quality was controlled through selecting the setpoints of the specific energy and the refining intensity at each stage. The targets of the specific energy and the refining intensity were achieved through control of the motor load, production rate and outlet consistency by manipulating the closing pressure, chip feed and dilution flow rate. The proposed control strategy is advantageous because it defines chip refining more accurately than conventional quality-energy control schemes, i.e., chip refining was defined by the specific energy alone. Therefore, the proposed control scheme provides more degrees of freedom in the control of the process. A constrained adaptive-predictive controller based Chapter 6. Constrained Multivariable Control Of A TMP Plant 185 Motor Load (kW) Plate Gap (x 0.1 mm) o > o > o ) o > o > a > o > o > a > Figure 6.52: Motor load setpoint unreachable Chapter 6. Constrained Multivariable Control Of A TMP Plant 186 Motor Load (kW) Plate Gap (x 0.1 mm) c » c » c » c » o o > a > o > o > Dilution Flow (kg/sec) Consistency (%) ui b> bo co co o to -p> cn oo o Figure 6.53: A plate clash was avoided by the proposed control scheme Chapter 6. Constrained Multivariable Control Of A TMP Plant 187 Figure 6.54: Response to the changes in consistency setpoint Chapter 6. Constrained Multivariable Control Of A TMP Plant 188 on non-square M I M O C A R I M A model was developed for the control of T M P refiners in the process. The control algorithm was based on the G P C because of its simplicity and ease of use along with its ability to handle problems in one algorithm. A n optimal solution of the constrained control problem was obtained by solving a QP problem. For the case of unfeasible constrained optimization problem, the M W L S was used as an alternative to calculate an optimal solution of the constrained control problem. Simulations of the proposed control schemes showed the successes in handling problems such as coupling, variable interactions, time-varying dynamics, time delay and external disturbances in one algorithm, without violations of input and output constraints. Control of a wood chip refiner in a T M P process is sometimes difficult due to the severe nonlinearity of the process, i.e., the incremental gain between refiner motor load and plate gap is subject to a sudden change in the sign. This is particularly critical when the refiner is operating close to the maximum load point. The linear approximation-based control for this case may lead to unstable operation. To avoid the problem, the nonlinear Laguerre functions were used to approximate process static nonlinearity and dynamics. The multivariable generalized predictive controller is then derived based on a mixed linear-nonlinear Laguerre model. The nonlinear Laguerre functions were used to represent the severe nonlinearity between refiner motor load and plate gap, while the linear Laguerre functions were used to approximate nonlinearities that are not severe. For an adaptive control implementation, the nonlinear Laguerre model can be arranged in a linear form in the parameters so that the powerful but simple Recursive Least Squares (RLS) algorithm can be used for the parameter estimation. Simulations of the proposed control scheme showed its successes in setpoint tracking as well as load disturbance reject. In the case that the incremental gain between refiner motor load and plate gain was subject to a sudden change in the sign due to a pulp pad collapse, the proposed controller had the ability to detect the collapse and quickly open the gap for the pad to Chapter 6. Constrained Multivariable Control Of A TMP Plant rebuild in order to avoid a plate clash. Chapter 7 Conclusions and Future Work 7.1 Conclusions This thesis has presented the development of a constrained multivariable predictive con-trol approach to the control of a T M P plant. Control of this particular process is difficult because of the multivariable nature, process constraints, nonlinearity, time-varying dy-namics and time delay along with stochastic disturbances. The objective of this work, therefore, has been to develop a control strategy suitable for handling these problems and for an industrial implementation towards an automated control of T M P manufacturing. The process was modelled using a combination of mechanistic and empirical methods. Static relationships between inputs and outputs were obtained based on first principles such as mass and energy balances to represent the quantitative relationships and complex interactions between variables. To avoid the complexity of the development of differential equations, the dynamics of the process were obtained through process identification. M i l l experiments were performed in Canadian mills from which a large number of on-line operational data was collected. Identification results showed that process dynamics can be reasonably represented by first order with time delay and that process disturbances feature nonstationary characteristics with first-order dynamics. The developed model was used for process pre-analysis, controller design and proposed controller investigation. Instead of the usual quality-energy control approaches, pulp quality was considered as a function of two variables: the energy per unit mass of pulp or specific energy and the 190 Chapter 7. Conclusions and Future Work 191 specific energy per bar impact or refining intensity. Therefore, pulp quality was controlled through selecting the setpoints of specific energy and refining intensity for each stage. The targets of specific energy and refining intensity at each stage were achieved through the control of motor load, production rate and refining consistency. The manipulated variables for obtaining the desired operating conditions were chip feed, dilution flow rate and closing pressure at the inlet of each refiner. In T M P manufacturing, the constraints on different flows, as well as constraints in-volved by equipment physical limitations and product quality requirements, all play a crucial role in control design. Applications of control action without consideration of the constraints may lead to a degraded process performance and, in some cases, closed-loop instability. To deal with the problems created by constraints, a constrained multivariable adaptive-predictive control scheme was developed for the control of a wood chip refiner of the T M P process. The control algorithm was based on G P C because it is easily under-stood and can handle the problems in one algorithm along with its ability to incorporate constraints explicitly. Future control actions were determined at each sampling time by minimizing the predicted errors subject to operating constraints. On-line estimates of a M I M O C A R I M A model were used for output prediction. The model structure was based on the identification results obtained from the process identification stage. A n analytical solution of the constrained M I M O G P C for NU = 1 was derived by solving a quadratic programming problem. In the case of unfeasible constraints, the mixed-weights least-squares algorithm was used to solve the constrained optimization problem. The MWLS-based solution will always converge to the point that minimizes the maximum constraint violation. The proposed control schemes were investigated on a combination of empirical and mechanistic models obtained from the modelling stage to approximate an industrial situation. Simulation results showed that the proposed schemes were able to handle the problems in one algorithm without violation of process constraints along Chapter 7. Conclusions and Future Work 192 with the potential of optimizing process operation and improving product quality. Control of a wood chip refiner is challenging due to the problem associated with the severe nonlinearity between refiner motor load and plate gap, i.e., the incremental gain between refiner motor load and plate gap is subject to a sudden change in the sign due to pulp pad collapse. This problem is particularly critical when the refiner is operating close to the maximum load point. To deal with the problem, a new approach was taken to represent the process, i.e., a M I M O Laguerre model was used to approximate process dynamics as well as static nonlinearity. Generalized predictive control solutions were then derived based on a mixed linear-nonlinear Laguerre function representation. The nonlinear Laguerre functions were used to approximate the severe nonlinearity between the refiner motor load and plate gap, whereas the linear Laguerre functions were used to represent nonlinearities that are not severe. The Laguerre function-based control scheme is advantageous because of its simplicity and flexibility while eliminating the need for assumptions about model order and time delay. Simulations of the proposed control scheme demonstrated that the nonlinearity and dynamics of the process can be represented properly even using lower order Laguerre functions. When the sign of the incremental gain between load and gap started to change, the proposed controller was able to detected the change and quickly open the gap to avoid plate clashes. The Laguerre model-based control scheme was considered as an alternative to the transfer function-based control when a wood chip refiner is operating close to the maximum load point. 7.2 Future Work It was difficult to address all the issues that arose during this research. A number of topics still need to be dealt with in the future. Chapter 7. Conclusions and Future Work 193 The proposed control schemes have been tested on a combination of empirical and mechanistic models to approximate an industrial plant. Further research could be under-taken towards in-depth investigation of the behaviour of the proposed control schemes over two stages in series. Throughout the study and simulation in this thesis, desired pulp quality was achieved through control of specific energy and refining intensity by adjusting the setpoints of motor load, throughput and outlet consistency. Further work would be undertaken by including quality feedback control. A n outer loop based on measurements of pulp proper-ties such as freeness, long fibre content and shive content would be included as a cascade control in providing the setpoints for inner loop control of specific energy and refining intensity at each stage. In this thesis, the chip screw speed was constrained about the desired production rate while still leaving some room for manipulation. A n extension of this study would consider maximizing the production rate for a given energy input to meet the different requirements of the process. Bibliography [1] G.A. Dumont, Y . Fu, and G. Lu. Nonlinear adaptive generalized predictive con-trol and applications. In Advances in Model-Based Predictive Control, volume 1, Department of Engineering Science, Oxford University, 1994. [2] G. A . Smook. Handbook for Pulp and Paper Technologists. Joint Textbook Com-mittee of the Paper Industry, 1992. [3] R. Franzen and R. Sweitzer. Refining forces in high stock concentration pulping. Appita, 36(2):116-121, 1982. [4] R. A . Leask. Mechanical pulping. In Pulp and Paper Manufacture, volume 2. the Joint Textbook Committee of the Paper Industry, 1987. [5] M.I . Stationwala, D. Atack, J.R. Wood, D.J . Wild, and A . Karnis. The effect of control variables on refining zone conditions and pulp properties. In 1979 Interna-tional mechanical pulping conference, Toronto, 1979. [6] K . B . Miles, H.R. Dana, and W.D. May. The flow of steam in chip refiners. In Inter-national Symposium on Fundamental Concepts of Refining, pages 30-42, Appleton, Wisconsin, Sept. 1980. [7] C. Earl Libby, editor. Pulp and Paper Science and Technology, volume 1. The Joint Textbook Committee of the Paper Industry, 1962. [8] J .C .W. Evans. T M P survey. Pulp and Paper, pages 99-110, June 1978. [9] R. J . Kerekes. Characterization of pulp refiners by a C-factor. Nordic Pulp and Paper Research Journal, 5(l):3-8, 1990. [10] K . B . Miles. Refining intensity and pulp quality in high-consistency refining. Paperi ja Puu-Paper and Timber, 72(5):508-519, 1990. [11] K . B . Miles, W.D. May, and A . Karnis. Refining intensity, energy consumption, and pulp quality in two stage chip refining. Tappi Journal, pages 221-229, March 1991. [12] J.D. Hoffmann and L . V . Welch. The Escher Wyss refiner a tool for pulp evaluation. In 1992 CPPA Pacific Coast Branch Spring Conference, Parksviille, B C , Canada, Apri l 1992. 194 Bibliography 195 [13] P.J. Leider and A . H . Nissan. Understanding the disk refiner-the mechanical treat-ment of the fibers. TAPPI Journal, 60(10), 1997. [14] B . C . Strand and N . Hartler. Modelling and optimization of full scale chip refining. In International Mechanical Pulping Conference, pages 46-54, Stockholm, Sweden, May 1985. [15] W.D. May, M.R. M c R A E , K . B . Miles, and W . E Lunan. A n approach to the mea-surement of pulp residence time in a chip refiner. Journal of Pulp and Paper Science, 14(3):J47-J53, 1988. [16] K . B . Miles. A simplified method for calculating the residence time and refining intensity in a chip refiner. Paperija Puu, 73(9):852-857, 1991. [17] W . E . Lunan, K . B . Miles, and W.D. May. The effect of differential pressure on energy consumption in thermomechanical pulping. Journal of Pulp and Paper Science, 11(5):J129-J135, 1985. [18] J.R. Rodarmel, M . J . Sferrazza, and P.E. Sharpe. Multiple speed refining evalua-tions, single/double disc refiners. In 76th CPPA Annual Meeting, Montreal, Feb. 1990. [19] W . C. Strand, O. Ferritsius, and A . V . Mokvist. Use of simulation models in the on-line control and optimization of the refining process. Tappi Journal, pages 103-108, Nov. 1991. [20] D. Church. Current and projected pulp and paper industry problems in process control and process modelling. In AIChE Symp. Series, volume 72, 1976. [21] R . B . Michaelson. A summary of thermomechanical pulping plant advanced control applications. In ISA PUPID/PMCD Symp., pages 65-78, Portland, Oregon, 1978. [22] J .W.T. Battershill and J .H. Rogers. On-fine computer installations: a survey. Pulp and Paper Canada, 81(8):44-52, 1980. [23] G.A. Dumont. Application of advanced control methods in the pulp and paper industry-survey. Automatica, 22(2):143-153, 1986. [24] J .H. Rogers, SC. Vespa, C. Mills, D.J . Butler, and N.D. Legault. Automatic control of chip refining. Pulp and Paper Canada, 81(10):T277-T282, Oct. 1980. [25] J .A . Nordin, J . Hi l l , and L. Nelvig. The SCS control system for the control of the mechanical pulping processes. In CPPA Process Control Conference, Halifax, 1980. Bibliography 196 [26] W.P. Alsip. Digital control of a T M P operation. Pulp and Paper Canada, 82(3):T76-T79, March 1981. [27] M . G . Horner and S. Korhonen. Computer control of thermomechanical pulping. Pulp and Paper Canada, 82(3):T72-T75, March 1981. [28] R .E . Jones and A . Pila. Process control for thermomechanical pulping. In IFAC In-strumentation and Automation in the Paper, Rubber, Plastics and Polymerization Industries, Antwerp, Belgium, 1983. [29] H.N. Koivo, P. Sammalisto, M . Perkola, and J . Renfors. A microprocessor-based control system for a T M P refiner. In the Mini and Microcomputers Conference, Sant Feliu, 1985. [30] G. A . Dumont. Self-tuning control of a chip refiner motor load. Automatic, 18(3):307-314, 1982. [31] S. Gendron, S. Banerjee, and L . Gerencser. Adaptive control using a Dahlin con-troller with application to wood chip refining. In IASTED International Symposium on Adaptive Control and Signal Processing, 1990. [32] G. A . Dumont and K . J . Astrom. Wood chip refiner control. IEEE Control System Magazine, 8(4):38-43, 1988. [33] B . J . Allison. Dual Adaptive Control of Chip Refiner Motor Load. PhD thesis, Department of Chemical Engineering, University of British Columbia, Vancouver, Canada, 1994. [34] W . R Cluett, J . Guan, and T .A . Duever. Control and optimization of T M P refiners. Control Systems'92, pages 1-4, 1992. [35] S.J. McQueen. Process identification and control of a T M P wood chip refiner. Mas-ter's thesis, Department of Chemical Engineering, University of Toronto, Toronto, 1995. [36] G.R. Fralic. Feedrate/motor load and dilution/consistency control of a T M P re-finer. Master's thesis, McGi l l University, Montreal, 1996. [37] A . Roche, J . Owen, K . Miles, and R. Harrison. A practical approach to the control of T M P refiners. In Control Systems'96, Halifax, Canada, 1996. [38] J . Hi l l , K . Saarinen, and R. Stenros. On the control of chip refining systems. In 1991 International Mechanical Pulping Conference, pages 235-241, Minneapolis, M N , 1991. Bibliography 197 A . Savonjousi, R. Stenros, and K . Saarinen. Automatic control of the consistency and density of thermo-mechanical pulp (TMP) . In Control systems' 90, pages 173-184. R. Evans, K . Saarinen, and R. Sutinen. Developments in sensors and controls strategies for refiners. In 78th CPPA Annual Meeting, pages A309-A314, Montreal, 1992. K . K . Hilden. On-line refiner consistency control increases production, efficiency and pulp quality. In 1994 Spring Conference, CPPA Pacific Coast and Western Branches, Jasper, Alberta, May 1994. K . Partanen and H . Koivo. On-line freeness control at a thermomechanical pulping plant. ISA, 1984. G. A . Dumont, N.D. Legault, J.S. Jack, and J .H. Rogers. Computer control of a T M P plant. Pulp and Paper Canada, 83(8):54-59, 1982. F . Cameron, H . N . Koivo, and K . Partanen. Self-tuning control of freeness in a thermomechanical pulping plant. In IAESTED Int. Symposium on Applied Control and Identification, ACI'83, Copenhagen, Denmark, 1983. H . T. Toivonen and J . Tamminen. Minimax robust L Q control of a thermomechan-ical pulping plant. Automatica, 26(2):347-351, 1990. G. A . Dumont. Self-tuning regulator: principles, present status and significance. Pulp and Paper Canada, 82(4):T117-T120, Apri l 1981. W . L . Bialkowski. Process control technology and practice: past, present and future. Pulp and Paper Canada, 91(1):T34-T41, Jan. 1990. B. J . Allison, J .E. Ciarniello, P.J-C. Tessier, and G.A. Dumont. Dual adaptive control of chip refiner motor load: Industrial results. Pulp and Paper Canada, 96(3):T73-T79, 1995. S.B.L. Kooi. Adaptive inferential control of wood chip refiner. Master's thesis, Concordia University, Montreal, 1991. S.B.L. Kooi and K . Khorasani. Control of wood chip refiner using neural networks. Tappi Journal, pages 156-162, June 1992. Y . Fu and G.A. Dumont. Chip refiner motor load adaptive control using a nonlinear laguerre model. In Second IEEE Conference on Control Applications, pages 371-376, Vancouver, Canada, 1993. Bibliography 198 [52] H . Du, G.A. Dumont, and Y . Fu. Nonlinear control of a wood chip refiner. In 4th IEEE Conference on Control Applications, Albany, New York, 1995. [53] H . Du and G.A. Dumont. Constrained multivariable control of a wood chip refiner. In 13th World Congress of IF AC International Federation of Automation Control, San Francisco, 1996. [54] M . Fournier, H . Ma, P . M . Shallhorn, and A . A . Roche. Control of chip refiner operation. Journal of Pulp and Paper Science, 18(5):J182-J187, Sept. 1992. [55] K . B . Miles and W.D. May. The flow of pulp in chip refiners. Journal of Pulp and Paper Science, 16(2):J62-J72, March 1990. [56] J .B. Cort and M . M . Barrette. Installation of blowline consistency control improves mill operation. In 81st CPPA Annual Meeting, Montreal, 1995. [57] J . Hi l l , H . Westin, and R. Bergsto. Monitoring pulp quality for process control. In 1979 International Mechanical Pulping Conference, pages 111-125, Toronto, 1979. [58] S. Rydefalk, T. Pettersson, E . Jung, and I. Lundqvist. The STFI optical fibre classifier. In 1981 International Mechanical Pulping Conference, Oslo, June 1981. [59] J . A . Nordin, M . Jackson, and H . Skold. Continuous on-line measurement of pulp quality. In 77th CPPA Annual Meeting, Montreal, 1991. [60] K . B . Miles and W.D. May. Predicting the performance of a chip refiner: a constitu-tive approach. In 1991 International Mechanical Pulping Conference, Minneapolis, M N , June 1991. [61] W . G . Mihelich, D.J . Wild, S.B. Beaulieu, and L.R. Beath. Single-stage chip refining - some major operating parameters and their effects on pulp quality. In 1971 International Mechanical Pulping Conference, Ottawa, June 1971. [62] B . Falk, M . Jackson, and O. Danielsson. The use of single and double disc re-finer configurations for the production of T M P for filled super calendered and light weight coated grades. In 1987 Mechanical Pulping Conference, pages 13-14, Van-couver, 1987. [63] O. Ferritisius and G. Gamte. Single and double disc refining at Stora Karnsveden. In 1989 Mechanical Pulping Conference, Helsinki, 1989. [64] M.I . Stationwala, K . B . Miles, and A . Karnis. The effect of first stage refining condi-tions on pulp properties and energy consumption. In 1991 International mechanical pulping conference, Minneapolis, M N , 1991. Bibliography 199 G.A. Smook. Handbook for Pulp and Paper Terminology. Angus Wilde Publica-tions, Vancouver, 1990. O.L. Forgacs. The characterization of mechanical pulps. Pulp and Paper Magazine of Canada, pages T89-T118, 1963. N . Ryt i and H . Manner. A control strategy for a thermomechanical pulping process. Paperi ja Puu, 10:640-651, 1977. G. Stephanopoulos. Chemical Process Control. Prentice-Hall, Inc., New Jersey, 1984. D.D. Ruscio. Topics in Model Based Control with Application to the Thermo Me-chanical Pulping Process. PhD thesis, The Dept. of Engineering Cybernetics, The Norwegian Institute of Technology, 1993. X . Qian and P.J.-C. Tessier. Modelling of a chip refining process and pulp quality prediction. In Post-Graduate Research Laboratory Report. Pulp and Paper Research Institute of Canada, Pointe Claire, Quebec, Oct. 1993. D. Atack and M.I . Stationwala. On the measurement of temperature and pressure in the refining zone of an open discharge refiner. Transactions of the Technical Section, l(3):71-76, Sept. 1975. B . Beaudry. Mechanical pulping. In Pulp and Paper Manufacture, volume 2, page 130. M . Redding. Effect of refiner dilution water flow rate on fibre characteristics and pulp properties. Master's thesis, Chemical Engineering, the University of British Columbia, Vancouver, 1992. A . J . Astrom and P. EykhofT. Introduction to stochastic control theory. In Mathe-matics in Science and Engineering, volume 70. G.E.P. Box and G . M . Jenkins. Time Series Analysis: forecasting and control Holden Day, San Francisco, 1970. B . J . Allison. Time series process model identification: Theory and application. Miscellaneous report (MR282), Pulp and Paper Research Institute of Canada, Van-couver, May 1994. L . Ljung. System Identification. Prentice-Hall, Inc., 1987. D.W. Clarke, editor. Advances in Model-Based Predictive Control. Oxford Science Publications, 1994. Bibliography 200 J . Richalet, A . Rault, J .L. Testud, and J . Papon. Model predictive heuristic control: applications to industrial processes. Automatica, 14:413-428, 1978. C R . Cutler and B .L . Ramaker. Dynamic matrix control - a computer control algorithm. In 1980 Joint Automatic Control Conference, volume 1, pages W P 5 - B , San Francisco, Aug. 1980. J .B. Froisy. Model predictive control: Past, present and future. ISA Trans., 33:235-243, 1994. D.W. Clarke and L . Zhang. Does long-range prediction work ? In IEE Control'85, Cambridge, 1985. M . Morari, C. E . Garcia, and D. M . Prett. Model predictive control: Theory and practice. In T .J . M c A V O Y , Y . Arkun, and E . Zafiriou, editors, Model Based Process Control, 1988. D. W. Clarke. Advances in model-based predictive control. In Advances in Model-Based Predictive Control, volume 1, University of Oxford, 1993. C E . Garcia, D . M . Prett, and M . Morari. Model predictive control: Theory and practice - a survey. Automatica, 25(3):335-348, 1989. J.R. Bosley, T .F . Edgar, A . A . Patwardhan, and G.T. Wright. Model-based control: A survey. In IFAC Advanced Control of Chemical Process, Toulouse, France, 1991. K . R . Muske and J.B. Rawlings. Model predictive control with linear models. AIChE J., 39(2), 1993. D. W. Clarke, C Mohtadi, and P.S. Tuffs. Generalized predictive control-part I. the basic algorithm. Automatica, 23(2): 137-148, 1987. P.S. Tuffs and D.W. Clarke. Self-tuning control of offset: a unified approach. Proc. IEE, 132D:100-110, 1985. D. W. Clarke, C. Mohtadi, and P.S. Tuffs. Generalized predictive control: Part II.-extensions and interpretations. Automatica, 23(2):149-160, 1987. D. W. Clarke and C. Mohtadi. Properties of generalized predictive control. Auto-matica, 25(6):859-875, 1989. A .R . Mcintosh, S.L. Shah, and D.G. Fisher. Analysis and tuning of adaptive generalized predictive control. The Canadian Journal of Chemical Engineering, 69(1):97-110, 1991. Bibliography 201 [93] S.L. Shah, C. Mohtadi, and D.W. Clarke. Multivariable adaptive control without a prior knowledge of the delay matrix. In Systems & Control Letter 9, pages 295-306. Elsevier Science Pulisher B .V . , North-Holland, 1987. [94] C. Mohtadi, S.L. Shah, and D.W. Clarke. Generalized predictive control of mul-tivariable systems. In Proc. of 5th Yale Workshop in Adaptive Systems Theory, pages 54-59, 1987. [95] T. Soderstrom, L . Ljung, and I. Gustavsson. A theoretical analysis of recursive identification methods. Automatica, 14:231, 1978. [96] M . E . Salgado, G.C. Goodwin, and R .H . Middleton. Modified least squares algo-rithm incorporating exponential resetting and forgetting. International Journal of Control, 47(2):477-491, 1988. [97] K . J . Astrom and B. Wittenmark. Adaptive Control. Addison-Wesley, 1989. [98] D . M . Prett and R.D. Gillette. Optimization and constrained multivariable control of a catalytic cracking unit. In 1980 Joint Automatic Control Conference, San Francisco, 1980. [99] R . K . Mutha. Constrained long range predictive control. Master's thesis, University of Alberta, Canada, 1990. [100] T .T .C . Tsang and W.D. Clarke. Generalized predictive control with input con-straints. IEE Proceedings, Pt.D, 135(6):451-460, Nov. 1988. [101] J . M . Dion, L . Dugard, A . Franco, N . M . Tri, and D. Rey. M I M O adaptive con-strained predictive control case study: A n environmental test chamber. Automat-ica, 27(4):611-626, 1991. [102] M . G . Singh and A . Ti t l i . Systems Decomposition, Optimization and Control. Perg-amon Press, 1978. [103] R. Fletcher. Practical Methods of Optimization. John Wiley and Sons, Inc., 1990. [104] C L . Lawson. Solving least squares problems. Prentice Hall, 1974. [105] M . S. Bazaraa and C. M . Shetty. Nonlinear Programming: Theory and Algorithms. Wiley, New York, 1979. [106] E .F . Camacho. Constrained generalized predictive control. IEEE Transactions on Automatic Control, 38(2):327-332, 1993. Bibliography 202 S. Gendron and M . Perrier. Constrained predictive control and model weighting estimation of low order plants. In 13th World Congress of IFAC International Federation of Automatic Control, San Francisco, June 1996. J . Rossiter and B . Kouvaritakis. Constrained stable generalized predictive control. Proceedings of the IEE-D, 140(4):243-254, 1993. L . Ljung and T. Soderstrom. Theory and practice of recursive identification. MIT Press, Cambridge, MA, 1983. S. Gendron. A class of algorithms for solving CDcontrol problems: convergence re-sults and simulated examples. In International CD Symposium, Tampere, Finland, June 1997. C L . Lawson. Contributions to the theory of linear least maximum approximations. Master's thesis, The University of California, Los Angeles, 1961. J.R. Rice. The approximation of functions: Nonlinear and multivariate theory, volume 2. Addison Wesley, 1969. J.R. Rice and K . H . Usow. The Lawson algorithm and extensions. Mathematics of Computation, 22:118-127, 1968. G.A. Dumont and Y . Fu. Nonlinear adaptive control via Laguerre expansion of Volterra kernels. International Journal of Adaptive Control and Signal Processing, 7(5):367-382, 1993. P. Eykhoff. System Identification. London, Wiley, 1974. C C Zervos and G.A. Dumont. Deterministic adaptive control based on Laguerre series representation. International Journal of Control, 48(6):2333-2359, 1988. Y . W . Lee. Statistical Theory of Communication. John Wiley, N Y , 1960. G.A. Dumont and C C Zervos. Adaptive control using Laguerre functions. In IFAC Workshop on Adaptive Control and Signal Processing, Lund, Sweden, 1986. G.A. Dumont, C C Zervos, and G. Pageau. Laguerre-based adaptive control of pH in an industrial plant extraction stage. Automatica, 26:781-787, 1990. M . Schetzen. The Volterra and Wiener Theories of Nonlinear Systems. John Wiley, N Y , 1980. P . M . Makila. Laguerre series approximation of infinite dimensional systems. Au-tomatica, 26(6):985-995, 1990. Bibliography 203 [122] R .E . King and N . Paraskevopoulos. Digital Laguerre filters. Circuit Theory and Applications, 5:81-91, 1977. [123] C.C. Zervos. Adaptive Control Based on Orthonormal Series Representation. PhD thesis, Electrical Engineering, the University of British Columbia, Vancouver, 1988. [124] G.J . Clowes. Choice of the time scaling factor for linear system approximations using orthonormal Laguerre functions. IEEE Trans, on Automatic Control, 10:487-489, 1965. [125] Y . Fu and G.A. Dumont. On determination of Laguerre filter pole through step or impulse response data. In IFAC 1993 World Congress, volume 5, pages 303-307, Sydney, Australia, 1993. [126] H . Nijmeijer and A . J . Van Der Schaft. Nonlinear Dynamic Control Systems. Springer Berlin, 1990. [127] B . Wittenmark. Adaptive control of a stochastic nonlinear system: an example. In 2nd Workshop on Adaptive Control: applications to nonlinear systems and robotics, Cancum, Mexico, 1992. [128] T . J . Harris and J.F. MacGregor. Design of multivariable linear-quadratic con-trollers using transfer functions. AIChE Journal, 33(9): 1481-1495, Sept. 1987. Appendix A Time Series Analysis Tools A . l Auto-Correlation Functions A stationary Gaussian time series yt can be uniquely described by its mean y and auto-covariance cyy(k), or equivalently by its mean y, variance o~y and auto-correlation function ryy(k). The auto-correlation function ryy(k) of yt at lag k is defined as: »•«,(*) = (A-291) v where yy(k) = jf - y)(yt+k - v) ^y = i^li(yt-y)2 (A.292) y = N S t l i yt where cyy{k) is auto-covariance of yt at lag A;. A; = 0,1, 2, • • •, K is the number of discrete sampling intervals between the observations of yt(t = 1,2, • • •, N). ay is the special case of cyy(k) at A: = 0 (i.e., a2. = cyy(Q)) and it presents the variance of the time series yt. It can be seen from above equation that cyy(k) = cyy(—k) and ryy(0) = 1 (i.e., auto-correlation is always unity at k = 0). Therefore, the auto-correlation ryy(k) is often plotted for the value of the lag k — 0,1, - • • ,K. The auto-correlation coefficients ryy(k) generally tend to zero as k increases. To have a 204 Appendix A. Time Series Analysis Tools 205 crude check on whether auto-correlation is effectively zero beyond a certain lag, the auto-correlations are usually plotted out with their 95% confidence intervals. A n approximate expression for the variance of estimated ryy(k) is defined in the following [75] on the assumption that the time series yt is completely random (uncorrelated or white). yav[ryy(k)} ~ 1 (A.293) Thus standard error (S.E.) of the estimated partial auto-correlation is S.E.[ryy(k)] ~ -j= (A.294) Thus 95 % confidence intervals are given by 2 95% Confidence Intervals ± 2S.E. = (A.295) Any auto-correlation coefficient ryy(k) which falls in 95% confidence intervals is insignif-icant. A.2 Partial Auto-Correlation Functions Partial auto-correlation function is a tool used to identify the order of an Auto-Regressive (AR) function which is fitted by an observed time series. The partial auto-correlation function of a time series yt at lag k may be estimated by fitting successively A R processes of order 1, 2, • • •, k shown in the following: ryyU) = <t>kylyryy(j - 1) + • • • + <ffi-»rm(j - k + 1) + <f>kykyryy(j - k) ^ j = l,2,---,k Appendix A. Time Series Analysis Tools 206 Solving above equations for k = 1,2,3, ••• successively by applying a simple recursive method and picking out the estimates (pyy, 4>yy, • • •, 4>yy. For an A R process of order n, the partial auto-correlation coefficient <j)yy will be: # * ^ 0 if k<n y y (A.297) $J = 0 if k>n In other words, the partial auto-correlation function of a nth-order A R process has a cutoff after lag k = n. The partial auto-correlation coefficients 4>yy(k) generally tend to zero as k increases. To have a crude check on whether the partial auto-correlation function is effectively zero beyond a certain lag, the partial auto-correlation coefficients are usually plotted out with their approximate 95% confidence intervals. The approximate 95% confidence intervals on the partial auto-correlations are defined as the same as in Equation A.295 under the assumption that the series yt is white. A.3 Cross-Correlation Functions The cross-correlation function ruy[k) between input ut and output yt is a tool used to model the input-output dynamics of the process by fitting the time series pair (ut,yt) to a model. The cross-correlation function ruy(k) is estimated by ruy(k) = ^ ® (A.298) 0~u(Jy where Cuy(^) £ Ei=i («t - u)(yt-h - y), k - 0,1,2, • • • [ jj E&*(y t - V)(ut-k -*),k = 0, - 1 , - 2 , • (A.299) Appendix A. Time Series Analysis Tools 207 cuy(k) is the cross-covariance coefficients between ut and yt (t — 1,2, • • •, JV) at lag k. u and y are the means of the series ut and yt respectively. o\ and o2 are the variances of the series ut and yt respectively, a2, is the special case of CuU(k) at k = 0 (i.e., a\ = cuu(0)), while a2 is the special case of cyy(k) at A; = 0 (i.e., cr2, = cyy(Q)). Since rUJ/(A;) is not generally equal to ruy(—k), the cross-correlation function is not symmetric about A; = 0, but in fact, the cross-correlation coefficients generally tend to zero in the range (—oo to i) or (i to + co). To have a crude check on whether the cross-correlations ruy(k) are effectively zero beyond a certain lag, the cross-correlations are usually plotted out with their 95% con-fidence intervals, calculated based on the assumption that ut and yt are uncorrelated. Upon this assumption, the variance of the estimated ruy(k) is given by [75] var[ruy(k)} ~ — - £ r t t u ( 0 r r o ( 0 (A.300) l=—oo Appendix B Optimal Solution via Analytical Quadratic Programming B.l Constraints For the special case of NU = 1, the constraints given in Equation 6.247 can be rewritten as follows: I A u 1 ) t | < Aui Wlmin < Ul,t < "1 max 2/lmin < Vl,t+j < 2/lmax I A u 2 , t | < A « 2 min max V2 min < V2,t+j < V2 max v / I A w 3 ) t | < A u 3 j = 1,---,JV 2 According to A u t = ut — ut-i, the amplitude constraints on the control signals can be modified as follows: Wlmin < Uitt = AuU + Uiit-1 < ^ lmax (B.302J W3min < «3,t = A « 3 , t + « 3 , t - l < «3max or «1 min - Ul,t-1 < A « M < « lmax ~ « l , t - l (13.303) M3min ~ « 3 , t - l < A U 3 ) t < « 3 m a x ~ « 3 , t - l The rate constraints on the control signals in Equation B.301 can be modified in a form as: —Aui < Auij < Aui -Au~2 < A u 2 ; t < Sua" (B.304) - A M 3 < A u 3 | t < A « 3 Since Au2,t = «2,t - « 2 , t - l = («2,t + V«2,t) - + V " 2 , t - l ) (B.305) 208 Appendix B. Optimal Solution via Analytical Quadratic Programming 209 the rate constraint on u2 is modified in a form of the control movement yet to be deter-mined: - A « 2 - Au5 > t + V « 2 , t - i < VU2,t < A M 2 - Aultt + V"2 , t -1 (B.306) The amplitude constraints on the outputs in Equation B.301 can be written over the predictive horizon iV2as follows: y i i and y2 min 2/l min l/l,t+l 2/1 max 2/l min < 2/l,t+2 < 2/l max — y i max 2/l min _ Ul,t+N2 2/l max 2/2 min V2,t+1 2/2 max 2/2 min < V2,t+2 < 2/2 max — y2 max V'2 min V2,t+N2 2/2 max (B.307) (B.308) Since (B.309) where y gives: y i = G'jUx + G' 2 u 2 + f x y 2 = G 3 u 2 + G 4 U 3 + f2 [2/t+i Vt+2 • • • Vt+mY• Substituting Equation B.309 into B.307 and B.308 y i m i n - f l < G i l i i + G ' 2 U 2 < yimax - f l y2min — f*2 < G ' 3 U 2 + G 4 U 3 < y 2 m a x _ $2 (B.310) Appendix B. Optimal Solution via Analytical Quadratic Programming 210 The above equation can be modified as follows: V l m i n - f l < [Gi G ' 2 0] y 2 m m - fa < [0 G 3 Gi] Auu V«2 , t V«2 , t A« 3,t < y i : f l < yai The above equation can be modified in a more compact form as A u i , t y i m i n — fl < G i G' 2 0 y2 m i n — ?2 0 G' 3 G i V « 2 , i A« 3,t Combining Equations B.303, B.304 and B.306 yields: ai < A « i i t < ft « 2 < V«2 , t < ft a 3 < A u 3 i t < ft where < y i max fl y2 max — ?2 (B.311) (B.312) (B.313) (B.314) ctx = m a x ( - A M 1 , « i m i n - «i,t-i) ft = min( At*!, « i m a x - «i,t-i) a 2 = - A « 2 - A i t 2 ) t + V u 2 , t - i ft = A w 2 - A«5 i t + VW2 , t - i a 3 = m a x ( - A u 3 , u 3 m i n - n 3,t-i) ft = min( A u 3 , w 3 m a x - u 3 , t - i ) B.2 Optimal Solution via Analytical QP The following is the derivation of the optimal solutions for the constrained control of the primary refiner for NU = 1. The 1st and 2nd order conditions to be satisfied are Appendix B. Optimal Solution via Analytical Quadratic Programming 211 rewritten in the follows. 1. The first order conditions are: dJ/d A uitt = 2hu A t t M + (h12 + h2i) y u2,t + 2cx + px - p2 = 0 dJ/d v «2,t = {h12 + h21) A u u + 2h22 y u2,t + (^23 + ^32) A u 3 ) t + 2c2 + p5 - Pe = 0 dJ/d A M 3 > T = (/i23 + ^32) y u2,t + 2h33 A w 3 j t + 2c 3 + p3 - p4 = 0 (B.315) 2. The second order conditions are: ciiPi = 0 for all i = 1 to 6 such that Oj = 0 if /ij < 0 (B.316) or a, > 0 if /ij = 0 The possible cases of rate and amplitude constraints are given as follows: Case 1: if ai is violated, then 2h33[(h12 + h21)ax + 2c2] - 2c 3(fr 2 3 + h32) {h23 + h32f - Ah22h33 {h23 + h32) c 3 ^ V ~ h~ ^^ 33 ft33 V«2,t A«3 , t Sufficient conditions are: «i > 0, i = 2 to 6 (B.317) Pi <o Case 2: if a 2 is violated, then A « i , t 2fe3s[(fei2 + M A + 2c2] - 2c3(h23 + fesa) (^ 23 + h32)2 - 4/l22^33 V«2,t A«3 , t Appendix B. Optimal Solution via Analytical Quadratic Programming 212 Sufficient conditions are: <k > 0 , i = 1,3,4,5,6 t*2 < 0 Case 3: if a3 is violated, then A (^12 + h2i)a2 + 2c x 2hn Vu2,t = Oi2 A _ (^23 + h32)a2 + 2C 3 ^ « 3 , t T T 2fl33 Sufficient conditions are: <k > 0,i = 1,2,4,5,6 A*3 < 0 Case 4: if a 4 is violated, then * _ (fci2 + M f t + 2C l i - i i t i j — —• 2hn V«2 ,t = ft A (/»23 + M f t + 2C3 2h33 Sufficient conditions are: ai > 0,* = 1,2,3,5,6 /z4 < 0 Case 5 : if a 5 is violated, then (h12 + h2i)[(h23 + h32)a3 + 2c2] - 4c i / i 2 2 Atti-4/ln^22 - (/ll2 + ^ 2 l ) 2 2/iu A 2ci V « 2 , t = - T — - r - A M M / l l2 + ^21 ' hl2 + ^21 A u 3 i t = a 3 (B.318) (B.319) (B.320) Appendix B. Optimal Solution via Analytical Quadratic Programming 213 Sufficient conditions are: <k > 0,i = 1,2,3,4,6 PB < 0 Case 6: if a 6 is violated, then (h12 + h2i)[(h23 + h32)(33 + 2c2] - 4c i / i 2 2 AMI,* = Ahnh22 - (h12 + h2i)2 2hn 2ci v«2,t = - y — T T - A U M -hi2 + h2i ' hi2 + h2i Au 3 , t = fh Sufficient conditions are: a>i > 0, i = 1 to 5 A*a < 0 Case 7: if ai and a3 are violated, then (B.321) (B.322) Aui,* = a i V w2,t = <*2 . /l23 + ^32 C 3 A u 3 l t = — " 2 - 7 — A I I 3 3 n33 Sufficient conditions are: <k > 0,i = 2,4,5,6 PuPs < 0 Case 8: if ai and a4 are violated, then A«i , t = «i V«2 , t = A A ^23 + ^32 a C3 A« 3 , t = JJ7 ^"33 "33 (B.323) Appendix B. Optimal Solution via Analytical Quadratic Programming 214 Sufficient conditions are: &i > 0, i = 2,3,5,6 Pi,p4 < 0 Case 9: if a x and a 5 are violated, then Auht - a i (hi2 + h21)a1 + (h23 + h32)ot3 + 2c2 2n22 Sufficient conditions are: <k > 0,i = 2,3,4,6 Pi,p5 < 0 Case 10: if a\ and a 6 are violated, then A«i,t = a i (hi2 + /i 2i)o;i + (h23 + h32)/33 + 2c2 V**< = 2 / ^ A«3,t = ft Sufficient conditions are: a* > 0,« = 2,3,4,5 Case 11: if a2 and a 3 are violated, then A« i > t = ft V«2 , t = "2 / l 2 3 + ^ 32 C 3 A« 3 , t = — a2 - — ^ " 3 3 f i 3 3 (B.324) (B.325) (B.326) Appendix B. Optimal Solution via Analytical Quadratic Programming 215 Sufficient conditions are: aj > 0, i = 1,4, 5,6 P2,P3 < 0 Case 12: if a2 and a 4 are violated, then A u l i t = p\ V«2 , t = r% = _ h 2 1 ± h B c^ ^ 3 3 ^33 Sufficient conditions are: <k >0,i= 1,3,5,6 p2,p4 < 0 Case 13: if a2 and are violated, then (B.327) (B.328) (/112 + h2i)Pi + {h23 + h32)a3 + 2c2 V ^ = 2 ^ A« 3 , i = 0=3 Sufficient conditions are: ^ > 0,i = 2,3,4,5 P2,P5 < 0 Case 14: if a2 and a 6 are violated, then (B.329) Appendix B. Optimal Solution via Analytical Quadratic Programming 216 A«i , t = ft V«2 , t = foi2 + h2i)l3i + (fea3 + M f t + 2c2 2/l22 A u 3 > t = ft Sufficient conditions are: Oj > 0, i = 1,3 )"2,A*6 < 0 Case 15: if a 3 and as are violated, then V«2 , t = "2 A u 3 > t = « 3 Sufficient conditions are: Oi > 0,i = 1,2 Case 16: if a 3 and are violated, then A«i,t A«i , t = V«2 , t = « 2 Au 3 , t = ft Sufficient conditions are: Appendix B. Optimal Solution via Analytical Quadratic Programming 217 Oj > 0, i = 1,2,4,5 (B.332) Case 17: if a 4 and 0,5 are violated, then Vu2,t = A A« 3 , t = « 3 Sufficient conditions are: <k > 0,t = 1,2 ^4 ,^5 < 0 Case 18: if a 4 and a6 are violated, then A«i , t A«i , t = -/^ 12 + •21 2hn V«2 , t = ft AMS,* = ft Sufficient conditions are: a, > 0, i = 1 to 5 A*4,A*e < 0 (B.334) Appendix C Constraint Mapping The following describes how to map the typical constraints for a SISO system to the expression in the form of future control increments. The amplitude constraints on the control signal over the control horizon Nu: (C.335) Umax U m i n < u(t + 1) < Umax Umin < U ( t + NU-1)< Umax According to Au(i) = u(i) — u(t — 1), above equation can be modified as: «min < U(t) = Au(t) + u(t - 1) < Wmax «min < U(t + 1) = Au(t + 1) + Au(t) + u{t - 1) < Wmax «min < U(t + NU-1) = Au(t + NU~1) + Au(t + NU- 1) + , • • • , + A U(t) < Umax (C.336) or 218 Appendix C. Constraint Mapping 219 Wmin - U(t - 1) < Au(t) < U m a x - u(t - 1) Wmin — U(t — 1) < Au(t) + Au(t + 1) < « m a x ~ u(t — 1) U m i n - u(t - 1) < Au(t) + Au(t + 1) + , ••• , + Au(t + Nu — 1) < Umax ~ u(t - 1) (C.337) Above equation can be rewritten in a more compact form as: or where Umin - A2u(t - 1) < AiU < U m a x ~ A2u(t - 1) AIVL < u m a x - A2u(t - 1) AxVL < - U m i n + A2u(t - 1) (C.338) (C.339) Ax = 1 0 ••• 0 1 1 ••• 0 1 1 ••• 0 (C.340) A2 = [l l..-l]T and u = [u(t) u(t + 1) • • • u(t + Nu - 1)]T, u m a x = [ u m a x u m a x • • • u m a x ] T and Umir [^min U m j n • • • U m j n ] • The rate constraints on the control signal over the control horizon Nu: Appendix C. Constraint Mapping 220 A u m i n < Au(t) < Aumax Aumia < Au(t + 1) < A w m a x A « m i n < Au(t + Nu - 1) < A u n Above equation can be written in a more compact form as: A u ^ < INUXNUH < A u m a x or (C.341) (C.342) -INUXNUH < -Au^ where A u m i n = [ A w m i n A t t m i n • • • A t t m i n ] T and A u m a x = [Aumax A u. (C.343) max A t i m a x T h e ampl i tude constraints on the output over the predictive horizon N2: ymin <y{t + N\) < J max ymin < y(t + JV1 + 1) < J max (C.344) ymin < y ( t + iV2)< j max Since y = G ' u + f where y = [y(t+l) y(t + 2) • • • y(t + N2)]T, above inequality equations can be modified to: 2/min 2/min leqG'u + f < 2/max 2/max' (C.345) 2/mir Appendix C. Constraint Mapping Above equation can be rewritten in a more compact form as: Ymin < G ' u + f < y m a x or G ' u < y m a x - f - G ' u < - y m a x + f where y m i n = [ymin 2/min • • • y m i n ] T and y m a x = [ymax ymax • • • Umax Combining Equation C.339, C.343 and C.347 gives: 221 At Umax - A2u(t - 1) -Ax - U m i n + A2u(t - 1) INUXNU u < A u m a x INUXNU - A u m i n G ' ymin f - G ' ymax ~T" f (C.346) (C.347) (C.348) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items