Applied Science, Faculty of
Electrical and Computer Engineering, Department of
DSpace
UBCV
Puwakkatiya-Kankanamge, Eranda Harinath
Doctor of Philosophy - PhD
University of British Columbia
The focus of this thesis is to develop an advanced control system for existing and to-be developed Thermo-Mechanical Pulping (TMP) refining processes. Therefore, the thesis has two parts. The first part developed two different but yet complementary closed-loop identification methods to update the models used in the low-level control layer of the current advanced control systems for TMP refiners. These newly developed identification methods are specific to closed-loop systems controlled by Model Predictive Control (MPC) techniques, and successfully tailored to the presence of the MPC. By updating the existing process models, the current advanced control systems can be operated to give better performance. However, these control systems may not be able to provide optimal performance due to process disturbances. Hence, a novel economically oriented advanced control system is developed for the existing two-stage TMP refiner processes for their optimal operation when disturbances are present. This novel technique dynamically optimizes the TMP processes as opposed to conventional two-layer methods which perform process optimization at steady-state, and has shown a potential economical benefit through reduction of total specific energy of a two-stage TMP process in a simulation study. In the second part, a novel TMP process with multi-stages of Low Consistency (LC) refining is studied for its optimal operation. The proposed Nonlinear Model Predictive Control (NMPC) technique dynamically optimizes the process and provides better performance when disturbances are present. This economically oriented NMPC (econNMPC) minimizes the total specific energy consumption of the process while respecting all the process constraints and achieving final desired pulp quality. In simulation studies, a TMP process with multiple stages of LC refining was able to save significant specific energy consumption with setpoint tracking control when disturbances are present. Moreover, further reduction of specific energy has been achieved with the developed econNMPC technique.
Identification, Control and Optimization Strategies for Thermo-Mechanical Pulping (TMP) Processes by Eranda Harinath Puwakkatiya-Kankanamge B.Sc., University of Moratuwa, 2003 M. Eng., Memorial University of Newfoundland, 2007 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) June 2012 c© Eranda Harinath Puwakkatiya-Kankanamge 2012 Abstract The focus of this thesis is to develop an advanced control system for existing and to- be developed Thermo-Mechanical Pulping (TMP) refining processes. Therefore, the thesis has two parts. The first part developed two different but yet complemen- tary closed-loop identification methods to update the models used in the low-level control layer of the current advanced control systems for TMP refiners. These newly developed identification methods are specific to closed-loop systems con- trolled by Model Predictive Control (MPC) techniques, and successfully tailored to the presence of the MPC. By updating the existing process models, the current advanced control systems can be operated to give better performance. However, these control systems may not be able to provide optimal performance due to process disturbances. Hence, a novel economically oriented advanced control system is developed for the existing two-stage TMP refiner processes for their optimal operation when disturbances are present. This novel technique dynamically optimizes the TMP processes as opposed to conventional two-layer methods which perform process optimization at steady-state, and has shown a potential economical benefit through reduction of total specific energy of a two-stage TMP process in a simulation study. In the second part, a novel TMP process with multi-stages of Low Consis- tency (LC) refining is studied for its optimal operation. The proposed Nonlinear Model Predictive Control (NMPC) technique dynamically optimizes the process and provides better performance when disturbances are present. This economi- cally oriented NMPC (econNMPC) minimizes the total specific energy consump- tion of the process while respecting all the process constraints and achieving final desired pulp quality. In simulation studies, a TMP process with multiple stages of LC refining was able to save significant specific energy consumption with set- point tracking control when disturbances are present. Moreover, further reduction ii Abstract of specific energy has been achieved with the developed econNMPC technique. iii Preface Chapter 3 is based on two published papers: 1. Rui Huang, Eranda Harinath, and Lorenz T. Biegler. Lyapunov stability of economically oriented NMPC for cyclic processes. Journal of Process Control, 21(4):501-509, 2011, and 2. Rui Huang, Lorenz T. Biegler, and Eranda Harinath. Robust stability of economically oriented infinite horizon NMPC that include cyclic processes. Journal of Process Control, 22(1):51-59, 2012. This research was initiated by Dr. Huang and Dr. Biegler, Department of Chemical Engineering, Carnegie Mellon University (CMU). Dr. Huang was a former Ph. D. student at CMU and worked with Dr. Biegler. Dr. Huang is currently working in the United Technologies Research Center, USA. I contributed to the analysis of infinite horizon NMPC and wrote sections of the manuscripts. The rest of the analysis and paper writing were done by Dr. Huang and Dr. Biegler. A part of Chapter 4 is based on a published paper, 1. Jun Yan, Eranda Harinath, Guy A. Dumont. Closed-loop identification for model predictive control: Direct method. Proceeding of the 48th IEEE Con- ference on Decision and Control, 2592-2597, Shanghai, China, Dec 16-18, 2009. The idea of use of explicit solutions of MPC stemmed from a discussion between myself and Dr. Yan, who was a post doctoral fellow at the Department of Electri- cal and Computer Engineering, University of British Columbia (UBC), and worked with Dr. Dumont. Dr. Yan is currently working in Hirain Technologies, China. I contributed to the analysis, section writing, and simulation studies. Dr. Dumont iv Preface kindly reviewed the theoretical proofs. Versions of Chapter 5 are published in two papers: 1. Eranda Harinath, Lorenz T. Biegler, and Guy A. Dumont. Control and op- timization strategies for thermo-mechanical pulping processes: Nonlinear model predictive control. Journal of Process Control. 21(4):519-528, 2011, and 2. Eranda Harinath, Lorenz T. Biegler, and Guy A. Dumont. Advanced Step Nonlinear Model Predictive Control for Two-stage Thermo Mechanical Pulp- ing Process. the 18th IFAC world congress, Milan, Italy, Aug. 28 - Sept. 02, 2011. I performed the complete analysis, simulation study, and paper writing for these research works. Dr. Biegler assisted to formulate the optimization problems and provided necessary softwares to solve the problems. Dr. Dumont and Dr. Biegler gave useful feedback during the research study and kindly reviewed the papers. A version of Chapter 6 is submitted for publication. 1. Eranda Harinath, Guy A. Dumont, and Lorenz T. Biegler. Predictive opti- mal control for thermo-mechanical pulping processes with multi-stage low consistency refining, 2012. I performed the complete analysis, simulation study, and paper writing for this research work. Dr. Biegler provided necessary softwares to solve the problems. Dr. Dumont and Dr. Biegler gave useful feedback during the research study and kindly reviewed the paper. The other papers which are published or submitted for publication during my Ph. D. study but not included in this thesis are: 1. Eranda Harinath, Jun Yan, and Guy A. Dumont. Advanced control of a TMP refiner line: A simulation study. the 95th Annual Meeting of PAPTAC, Montreal, Quebec, Canada, Feb 2-3, 2010. v Preface In this research work, I performed the complete analysis, simulation study, and paper writing. Dr. Dumont and Dr. Yan gave useful feedback during the research study and kindly reviewed the paper. 2. Jun Yan, Eranda Harinath, and Guy A. Dumont. Tuning and Identification for Model Predictive Control: An Iterative Approach. In Proceedings of IFAC Symposium on System Identification, St. Malo, France, July 6-8, 2009. This research was initiated by Dr. Yan. I contributed to the analysis and simulation studies. Dr. Dumont kindly reviewed theoretical proofs and gave useful feedback during the research study. 3. X. Shi, Eranda Harinath, D. M. Maijer and Guy A. Dumont, Nonlinear Model Predictive Control in a Low Pressure Die Casting Process, submit- ted for publication, 2012. This research was initiated by X. Shi, Dr. Maijer and Dr. Dumont. Dr. Maijer is with the Department of Materials Engineering at the University of British Columbia, and X. Shi is a Ph. D. student of Dr. Dumont and Dr. Maijer. In this research, I contributed to the analysis of system identification and designing of nonlinear controller. vi Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Two-stage TMP Process . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Research Motivation . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Research Objective . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Summary of Contributions . . . . . . . . . . . . . . . . . . . . . 7 1.5 Organization of the Thesis . . . . . . . . . . . . . . . . . . . . . 7 2 Model Predictive Control (MPC) . . . . . . . . . . . . . . . . . . . 9 2.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Receding Horizon Control . . . . . . . . . . . . . . . . . . . . . 9 2.3 Basic NMPC Problem Set-up . . . . . . . . . . . . . . . . . . . 10 vii Table of Contents 2.4 Nominal Stability of NMPC . . . . . . . . . . . . . . . . . . . . 12 2.5 Robust Stability of NMPC . . . . . . . . . . . . . . . . . . . . . 17 2.5.1 Inherent Robustness of NMPC . . . . . . . . . . . . . . 17 2.5.2 Robust Design of NMPC . . . . . . . . . . . . . . . . . 19 2.6 Solution of NMPC Optimization Problem . . . . . . . . . . . . 20 2.6.1 IPOPT . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.6.2 asNMPC . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Economically Oriented NMPC . . . . . . . . . . . . . . . . . . . . . 25 3.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Infinite Horizon economically oriented NMPC (econNMPC) . . 27 3.4 Nominal Stability of the econNMPC . . . . . . . . . . . . . . . . 29 3.5 Robust Stability of the econNMPC . . . . . . . . . . . . . . . . 32 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4 Closed-loop Identification with Constrained MPC . . . . . . . . . . 36 4.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3 Output Feedback MPC . . . . . . . . . . . . . . . . . . . . . . . 38 4.3.1 Output Prediction . . . . . . . . . . . . . . . . . . . . . 39 4.3.2 Explicit Solution of MPC . . . . . . . . . . . . . . . . . 40 4.4 Prediction Error Method (PEM) . . . . . . . . . . . . . . . . . . 42 4.5 Informativity of Closed-loop Data: Direct Method . . . . . . . . 42 4.5.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.6 Joint Input-Output Method . . . . . . . . . . . . . . . . . . . . . 47 4.6.1 Two-stage Method . . . . . . . . . . . . . . . . . . . . . 48 4.6.2 First Stage: Projection Method . . . . . . . . . . . . . . 49 4.6.3 Example . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.6.4 Second Stage: Long Range Predictive Identification . . . 52 4.6.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 viii Table of Contents 5 Control and Optimization of Two-stage Thermo Mechanical Pulping (TMP) Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.2.1 TMP Process Modeling . . . . . . . . . . . . . . . . . . 60 5.3 NMPC for Setpoint Tracking . . . . . . . . . . . . . . . . . . . . 65 5.3.1 Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.3.2 Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.4 NMPC for Economical Operation . . . . . . . . . . . . . . . . . 67 5.4.1 Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.4.2 Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6 Control and Optimization of TMP Processes with Multi-stage Low Consistency Refining . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.3 Process Description and Modeling . . . . . . . . . . . . . . . . . 82 6.3.1 Latency Chest Modeling . . . . . . . . . . . . . . . . . . 82 6.3.2 LC Refiner Modeling . . . . . . . . . . . . . . . . . . . 83 6.4 Scheme A: Conventional TMP Process with a Third-stage LC Re- finer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.4.1 NMPC for Setpoint Tracking and Economical Operation . 86 6.4.2 Closed-loop Simulation . . . . . . . . . . . . . . . . . . 87 6.5 Scheme B: Improved TMP Process with Single-stage of HC and Multi-stage LC Refiners . . . . . . . . . . . . . . . . . . . . . . 88 6.5.1 NMPC for Setpoint Tracking and Economical Operation . 89 6.5.2 Closed-loop Simulation . . . . . . . . . . . . . . . . . . 90 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 ix List of Tables 4.1 Three different case studies and respective control laws for Exam- ple 4.5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2 The complexity criterion F for the cases listed in Table 4.1 of Example 4.5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3 Model estimates for the cases listed in Table 4.1 of Example 4.5.1 47 5.1 Time constants and time delays of subprocesses of the TMP process 64 6.1 The chip bulk density and chip solid content changes during the closed-loop simulation of the TMP process in Scheme A . . . . . 87 x List of Figures 1.1 Basic schematic of a two-stage TMP process . . . . . . . . . . . . 2 1.2 Current control strategy of TMP refiner processes: the conven- tional two-layer architecture . . . . . . . . . . . . . . . . . . . . 5 2.1 The basic of MPC . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.1 The Bode plots of the true process G, and its estimated process models ¯GA, ¯GB, and ¯GC for the case studies A, B, C, respectively in Example 4.5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 The Bode plots of the true noise process H , and its estimated pro- cess models ¯HA, ¯HB, and ¯HC for the case studies A, B, C, respec- tively in Example 4.5.1 . . . . . . . . . . . . . . . . . . . . . . . 49 4.3 The Bode plots of the true process G0, the initial model Ginit , the estimated transfer functions G1 and GLRPI in Example 4.6.3 and 4.6.5, respectively . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.1 Controlled variables of the TMP process: Performance compari- son of asNMPC and ideal-NMPC techniques for setpoint tracking problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Manipulated inputs of the TMP process: Performance compari- son of asNMPC and ideal-NMPC techniques for setpoint tracking problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.3 The production rate, and motor load responses of the TMP process to the change of the chip bulk density . . . . . . . . . . . . . . . 72 xi List of Figures 5.4 The specific energies of the both primary and secondary TMP pro- cess: performance comparison of the sptNMPC with econNMPC techniques to the change of the chip bulk density . . . . . . . . . 73 5.5 The freeness, long fiber content and shive content responses of the TMP process to the change of the chip bulk density. . . . . . . . . 74 5.6 The production rate, consistencies and motor loads responses of the TMP process to the change of the chip bulk density and chip solid content, with random disturbances on the state variables . . . 77 5.7 The specific energies of the both primary and secondary the TMP process: performance comparison of the sptNMPC and econN- MPC techniques to the change of the chip bulk density and chip solid content, with random disturbances on the state variables . . . 78 5.8 The freeness, long fiber content and shive content responses of the TMP process to the change of the chip bulk density and chip solid content, with random disturbances on the state variables. . . . . . 79 6.1 Scheme A: basic schematic of a conventional TMP process with a third stage LC refiner . . . . . . . . . . . . . . . . . . . . . . . . 84 6.2 The manipulated input variables of the TMP process in Scheme A. 92 6.3 The controlled state variables of the TMP process in Scheme A. . 93 6.4 The specific energy consumptions of the primary, secondary, and tertiary refiners of TMP process in Scheme A. . . . . . . . . . . . 94 6.5 The pulp quality variables after the tertiary refining of the TMP process in Scheme A. . . . . . . . . . . . . . . . . . . . . . . . . 95 6.6 Scheme B: basic schematic of an improved TMP process with a two stages of LC refining . . . . . . . . . . . . . . . . . . . . . . 96 6.7 The manipulated input variables of the TMP process in Scheme B. 97 6.8 The controlled state variables of the TMP process in Scheme B. . 98 6.9 The specific energy consumptions of the primary refiner, tertiary refiner 1, and tertiary refiner 2 of TMP process in Scheme B. . . . 99 6.10 The pulp quality variables after the tertiary refining of the TMP process in Scheme B. . . . . . . . . . . . . . . . . . . . . . . . . 100 xii List of Figures 6.11 Specific energies of both Schemes A and B with sptNMPC and econNMPC techniques. . . . . . . . . . . . . . . . . . . . . . . . 101 xiii Symbols ∈ belongs to ⊆ subset of ∑ summation ‖x‖ norm of a vector x QT (xT ) transpose of a matrix Q (a vector x) ‖x‖2Q weighted norm xT Qx R > 0 a positive definite matrix R min minimum ∀ for all 6= not equal R set of real numbers R≥0 non-negative real numbers R n n-dimensional Euclidean space f : Rn×Rm → Rn a function f mapping a space Rn×Rm into a space Rn /0 empty set ln natural logarithm ∇ f gradient of a function f ∇2L Hessian of a function L diag(w) a diagonal matrix with w vector in the diagonal lim limit ∞ infinity → tends to deg degree xiv Abbreviations TMP Thermo Mechanical Pulping HC High Consistency RTO Real-Time Optimization LC Low Consistency NMPC Nonlinear Model Predictive Control econNMPC economically oriented NMPC LTI Linear Time Invariant NLP Nonlinear Programming IPOPT Interior Point OPTimizer asNMPC advanced step NMPC MPC Model Predictive Control MPHC Model Predictive Heuristic Control DMC Dynamic Matrix Control GPC Generalized Predictive Control CARIMA Control Auto-Regressive Integrated Moving-Average ISS Input to State Stability RPI Robustly Positively Invariant KKT Karush Kuhn Tucker LICQ Linear Independent Constraint Qualification SSOC Sufficient Second Order Condition SC Strict Complementary PEM Prediction Error Method LRPI Long Range Predictive Identification QP Quadratic Programming FIR Finite Impulse Response ARX Auto-Regressive Exogenous xv Abbreviations RLS Recursive Least Square MIMO Multi-Input Multi-Output SISO Single-Input Single-Output CSF Canadian Standard Freeness LFC Long Fiber Content SC Shive Content sptNMPC setpoint tracking NMPC AMPL A Modeling Language for Mathematical Programming MATLAB MATrix Laboratory MHE Moving Horizon Estimation EKF Extended Kalman Filter BJ Box-Jenkins xvi Acknowledgements First, I would like to offer my sincere gratitude to my advisor, Prof. Guy A. Dumont at the Department of Electrical Engineering, the Faculty of Applied Science, for giving me an opportunity to pursue my Ph.D. studies. I express my sincere grat- itude to him for providing an excellent research environment, guidance, valuable discussions, encouragement, and friendly attitude during the course of study. I also like to thank Prof. James A. Olson, the Director of Pulp and Paper centre (PPC), UBC, for providing me the financial support, guidance, and valuable discussions throughout the course. It was an honor to work with many research collaborators who have contributed significantly to this research work, and their help and guidance are greatly ac- knowledged. Many thanks to Prof. Lorenz T. (Larry) Biegler in the Department of Chemical Engineering at Carnegie Mellon University (CMU), and Dr. Rui Huang. I also like to thank Dr. Rodrigo López-Negrete and all the other members in the Centre for Advanced Process Decision-making (CAPD) at CMU, who helped me during my research visit. I also want to thank Prof. Frank Allgöwer in the Institute for Systems Theory and Automatic Control (IST) at University of Stuttgart, and Prof. Christian Ebenbauer, Dr. Shuyou Yu, Christoph Maier, and all other mem- bers in IST, who helped me during my research visit. Many thanks to Dr. Jun Yan who is now with Hirain Technologies, China, and XinMei Shi for having fruitful discussions. I gratefully acknowledge the supports given by all other members in the PPC. I gratefully acknowledge the financial contributions and supports from the Michael Smith Foreign Study Supplements Scholarship, NSERC, British Columbia Innovation Council, Advanced Fiber Technologies, Andritz, Arkema, BC Hydro, Canfor, Catalyst Papers, CEATI/ Ontario Power Authority, Honeywell, Howe Sound Pulp and Paper, Paprican, Westcan, and West Fraser Quesnel River Pulp. xvii Dedication To My Parents, Lisie and Peter To My Sisters To My Loving Wife Ruvini xviii Chapter 1 Introduction Thermo-Mechanical Pulping (TMP) is a dominant process in the pulping industry. In the TMP process, wood chip refiners, which are required to run uniformly, are used to convert wood chips to pulp. The TMP refiners consume large amounts of electric energy - in British Columbia, Canada, around 10% of the total electrical energy production is consumed by seven mechanical pulping plants [1]. Conven- tional TMP processes have two stages of High Consistency (HC) refining. The consistency can simply be defined as mass ratio of dry-fiber to the mixture of dry- fiber and water. The HC refiners play the main role in achieving the paper-making properties of pulp, yet are high energy consuming; nearly 80% of the total energy consumption in a TMP process is due to the HC refiners [1]. Although the TMP process is a mainstay of the pulping industry [2], TMP mills face a number of challenges due to the increase in electricity costs, and initiatives to reduce green house gas emissions. TMP mills need to reduce electrical energy consumption, thereby reducing production cost and environmental impact [3] [4]. In addition, pulp mills have to produce improved pulp with less variability in pulp quality in order to maintain competitiveness and retain customers [5]. At present, TMP mills are seeking advanced process control techniques to minimize energy consumption and improve pulp quality in order to address the above challenges. 1.1 Two-stage TMP Process Fig. 1.1 shows the basic operational units of a typical TMP process. There are mainly three operational areas in the TMP process [6]. In the wood chips pre- treatment, wood chips are first screened to remove under and over sized particles, and are then steamed at atmospheric pressure around 100oC for preheating. Fol- lowing this, the wood chips are fed to the washing system, which operates around 1 1.1. Two-stage TMP Process Washing, Preheating Screening, Steaming 1st stage HC refining Wood chips Pulp Wood chips pretreatment Refining 2st stage HC refining Screening, Reject refinging Latency removal Dewatering Pulp processing Figure 1.1: Basic schematic of a two-stage TMP process 70− 85oC to maintain chip moisture, to remove contaminants like sand, saw-dust and stones from the wood chips [2]. Before feeding to the primary refining stage, the wood chips are preheated to soften lignin binding wood fibers in order to facil- itate their further separation, maximizing long fiber content and minimizing shive content [6]. The smooth operation of wood chip refiners is critical, as they play a major role in achieving paper-making properties of the pulp. There are two stages of HC refining: primary and secondary. The wood chips are fed to the primary refiner using a screw feeder which is manipulated to control the production rate. The primary refiner is used to convert wood chips to individual wood fibers; the refining plates of the refiners are made such that wood chips are gradually broken into fibers in the refining zone. The second stage of refining is used to develop the pulp qualities. The primary and secondary refiners are operated at HC, and dilution water is fed to the refiners to control the consistencies. The pulp processing, done at Low Consistency (LC), is performed to achieve the required properties of the pulp. After the secondary stage of HC refining, the pulp is fed to the latency chest for latency removal, i.e. relieving of internal stresses and removal of curl, resulting in improved pulp properties such as strength. The screening and reject refining are then performed to remove contaminants and de- crease shive content of the pulp. In the dewatering process, optical properties of 2 1.1. Two-stage TMP Process the pulp are improved. The refining process is, inherently, a Multi-Input Multi-Output (MIMO) pro- cess with complex dynamics and severe interactions among process variables. To understand the refining process with a view to controlling it, the process variables are divided into four categories: manipulated variables, operating variables, pulp- quality variables, and disturbance variables [6, 7]. Manipulated Variables (MVs) The Manipulated Variables (MVs) are the main input variables for any refiner, and must be controlled in order to achieve the desired operational conditions. In the case of primary and secondary refining, the MVs include chip transfer screw speed (considered at the primary refiner only), plate gap, and dilution water flow [8, 9]. The temperature of the wood chips entering the primary refiner can also be considered as an input variable and assumed to be the same as the pre-heater temperature which softens the wood chips. However, the temperature of the wood chips entering the primary refiner is not considered as a MV [2, 7]. Operating Variables (OVs) The Operating Variables (OVs) are the key variables, and are highly correlated with pulp quality. By controlling the MVs, the OVs are kept at optimum setpoints (OSPs); the OSPs are determined by economical operation of the refining process. The main OVs in all stages of refining are: motor load, consistency, and production rate. Among them, the motor load is believed to be an important OV [7, 10, 11]. In the literature, specific energy and refining intensity are also considered as OVs; however, they are correlated with the above OVs for a given refiner. Therefore, the effect of specific energy and refining intensity on pulp quality can be controlled through motor load, consistency and production rate [12, 13]. Pulp Quality Variables (PQVs) The Pulp Quality Variables (PQVs) define the paper-making properties of the pulp, and the choice of different properties depends on the end users. The main PQVs that characterize the refining process are: freeness, fiber size distribution, and shive 3 1.2. Research Motivation content [11, 12]. Online measurement of PQVs is possible with a Pulp-Quality Monitor (PQM), and is important to determine the optimum economical operation of the refiner units. Disturbance Variables (DVs) The variables that cannot generally be manipulated with and are independent of the input variables are called Disturbance Variables (DVs). These variables appear in the refining process mainly through wood chips and the refiner unit itself. The DVs influence the operating conditions of the refiner process, thereby changing the pulp quality. The important DVs in wood chips are wood-chip species, chip size distribution, moisture content, and chip-bulk density. The back-flow of steam from the refiner and refiner plate wear are the main disturbances in the refiner unit [6, 7]. A comprehensive understanding of these process variables is necessary to de- velop refiner process control strategies. 1.2 Research Motivation Chips refiners are the most energy intensive units in the TMP process; approxi- mately 80% of the total energy is consumed by HC refiners [1]. Therefore, there is considerable literature on implementation of advanced control systems for them [5], [14], [8]. However, it has been reported that advanced control systems have failed to perform well over for a long period: over 90% of all advanced control applications have been turned off within a year [5]. Major reasons for poor per- formance of advanced control systems are 1) less attention being paid to low-level regulatory control and 2) the linearity assumption of the TMP process. In addition, operator training and servicing after installation need to be improved for reliable and continuing operation of advanced control systems. Multi-layer architecture is commonly used in the decision making hierarchy of industrial processes, and economical operation of a process is usually performed using an advanced process control system with a two-layer structure. Fig. 1.2 shows the two-layer control strategy for a TMP process. The two-layer structure consists of a higher layer, implementing Real-Time Optimization (RTO), followed 4 1.2. Research Motivation Real-Time Optimization (RTO) Layer Objective: Minimize specific energy consumption Strategy: Steady state optimization Setpoint Tracking Layer Objective: Track the setpoint with minimum transient time Strategy: Model Predictive Control (MPC) Refiner Process Primary, secondary, and tertiary refiner processes Figure 1.2: Current control strategy of TMP refiner processes: the conventional two-layer architecture 5 1.2. Research Motivation by a lower layer for setpoint tracking or regularization control. The RTO layer computes optimal steady-state setpoints using a steady-state process model, and provides these setpoints to the setpoint tracking layer, in order to maximize an economical objective while keeping product quality constraints within safe values. The performance of the advanced control systems installed in TMP mills di- rectly rely on the basic low-level regulatory control system; only the basic low- level control system has direct access to the TMP refiner process and can change the control input variables [5, 15], such as the refiner plate gap. Therefore, the proper monitoring of the regulatory layer is essential for economical operation of the TMP process. Most existing advanced control systems installed in TMP pro- cesses are based on a linear model of the process. Periodic update of this linear pro- cess model would be required for satisfactory performance of the controller, as the performance of the closed-loop system may deteriorate due to model/process mis- match caused by physical changes of the process over time, such as wear and tear of the refiner plates. However, setpoint tracking control would be difficult with a controller based on a linear model since TMP refiner processes are inherently com- plex processes, consisting of nonlinear dynamics. Moreover, the inconsistency of the two models used in the RTO and setpoint layers could result inconsistent and unreachable setpoint [16, 17]. Furthermore, optimum setpoints are not always at steady-state due to dynamically changing disturbances, like chip bulk density, chip moisture content, and model-plant mismatch. Therefore, it is sub-optimal to mini- mize transition time between setpoints. These problems have been noticed in one of our industrial partner’s TMP control system. On the other hand, LC refining technology is emerging. More attention has been paid towards research and devel- opment of LC refiners, mainly because the LC refiners are more energy efficient than HC refiners [1, 18]. The inclusion of the third stage LC refiner in the mainline has led to the reconsideration of advanced control system to optimally coordinate the HC and LC refiners. Therefore, time is ripe to develop an advanced control technique for the existing TMP processes with two-stage HC refiners, and new TMP processes which are to be developed with HC and LC refiners. 6 1.3. Research Objective 1.3 Research Objective The main objective of this thesis is to propose strategies to improve the perfor- mance of the existing TMP processes and the new TMP processes. For that pur- pose, closed-loop identification techniques, and nonlinear model based control and optimization strategies are presented for TMP processes. 1.4 Summary of Contributions The following are the summary of contributions in this research. 1. Development of two closed-loop identification schemes. The first scheme is a novel attempt while the second scheme is an extension of an existing method with a meaningful identification criteria. Both the schemes can be applied for Linear Time-Invariant (LTI) systems controlled by constrained MPC laws. 2. Design of a novel control and optimization strategy based on Nonlinear Model Predictive Control (NMPC) technique for a two-stage TMP process. 3. Design of economically oriented control strategy for a novel TMP process which consists of HC and LC refiners using NMPC framework. This re- search work is one of the first contributions that investigates dynamic opti- mization of TMP processes with multi-stages of LC refiners. 4. Analysis of nominal and robust stabilities of an infinite horizon NMPC with a general economic objective function. Since the stability analysis of this class of problem is still in early stage of its development, this research work will be a significant contribution to the field. 1.5 Organization of the Thesis The rest of the thesis is organized as follows. Chapter 2 presents the basic princi- ple of NMPC techniques and introduces some definitions and assumptions which are used throughout the thesis. The material presented in the Chapter 2 serves as 7 1.5. Organization of the Thesis the basis for the subsequent chapters. Chapter 3 discusses economically oriented NMPC (econNMC) framework for the optimal operation of a process. In partic- ular, infinite horizon econNMPC with a general economic objective function is discussed. The nominal and inherent robustness properties are discussed for econ- NMPC in Chapter 3. Chapter 4 presents closed-loop identification methods which can be applied to an LTI system controlled by a constrained MPC law. In Chapter 5, economical operation of two-stage TMP process is studied. The setpoint track- ing and econNMPC techniques are discussed, and closed-loop simulation studies also presented for two-stage TMP processes. Chapter 6 studies two types of TMP processes for optimal operation: a conventional two-stage TMP process with a LC refiner, an improved TMP process with two-stages of LC refining. Both pro- cesses are dynamically optimized with econNMPC techniques. Simulation studies are presented for both processes to realize potential economical benefits through reduction in total specific energies. 8 Chapter 2 Model Predictive Control (MPC) 2.1 Summary In this chapter, the basic principles of Nonlinear Model Predictive Control (NMPC) are presented for discrete time nonlinear systems. The nominal and robust stability of setpoint tracking NMPC problems are reviewed. Solution strategies for the re- sulting Nonlinear Programing (NLP) problem of NMPC framework are discussed; specifically, the basic formulation of the Interior Point OPTimizer (IPOPT) solver is presented. Furthermore, the advanced step NMPC (asNMPC) concept, which is a fast approximate solution technique of the resulting NLP to reduce on-line com- putational time, is discussed. The review of NMPC presented in this chapter is not exhaustive, but rather provides the basis for the following chapters. 2.2 Receding Horizon Control In MPC, at each sample time k, first the future dynamic behavior of the process is predicted using a model of the process over a prediction horizon N. Next the optimal control moves over a control horizon M are determined by minimizing an objective functional subject to the process dynamics, and input and state con- straints. Although M control moves are calculated, only the first control move is applied to the process until the next sample time k+1. At sample time k+1 , the process evolves to new states which, in practice, will not be equal to the model predicted states at sample time k+ 1, due to model-process mismatch and distur- bances. Hence, a new constrained optimization problem is formulated over the prediction horizon N with the state information at sample time k+1 and is solved for M optimal control moves. This strategy is also known as receding horizon or moving horizon control, and the MPC law is implemented in this receding horizon 9 2.3. Basic NMPC Problem Set-up k k+1 k+2 k+N Prediction horizon Time (a) At sample time k k k+1 k+2 k+N+1 Prediction horizon Time Predicted trajectory Predicted control input Applied control input Actual trajectory Past Future k+N (b) At sample time k+1 Figure 2.1: The basic of MPC framework. The basic concept of MPC is illustrated in Fig. 2.1 and summarized in the following algorithm. Algorithm 1 The basic steps of MPC framework at sample time k 1: Predict the future dynamic behavior of the process over the prediction horizon using the state information at sample time k. 2: Determine the optimal control inputs over the control horizon by solving the constrained optimization problem. 3: Apply only the first control move to the process until next sample time k+1. 4: k → k+1, go the step 1 2.3 Basic NMPC Problem Set-up Consider a discrete time nonlinear system described by xk+1 = f (xk,uk), k ≥ 0, (2.1a) 10 2.3. Basic NMPC Problem Set-up subject to the control and state constraints, uk ∈U and xk ∈X, (2.1b) where xk ∈ Rn and uk ∈ Rm are the state and control input variable vectors of the system, respectively. U ⊆ Rm is the control constraint set, and X ⊆ Rn is the state constraint set. f : Rn ×Rm → Rn is the nominal nonlinear system function. In the sequel, z ∈ Rn and v ∈ Rm are used to denote the model predicted state and computed control input vectors in the controller, respectively. Using the current process state xk as initial conditions, the dynamic optimiza- tion problem for NMPC is formulated at sample time k as follows. min v J(z(·),v(·)) (2.2a) s. t. : zi+1 = f (zi,vi), z0 = xk i = 0, . . . ,N−1, (2.2b) zi ∈ X i = 0, . . . ,N, and (2.2c) vi ∈ U i = 0, . . . ,N−1 (2.2d) where N is, for simplicity, the length of the horizons in sample times, without differentiating the control and the prediction horizons. xk is the current, measured or estimated, state vector of the process. v = {v0,v1, . . . ,vN−1} ∈ UN(xk) where UN(xk) is the set of feasible control moves which satisfy the constraints (2.2b)- (2.2d). J : Rn×Rm → R is a cost functional, sometimes called performance index or objective functional, is defined over the horizon N, J(z(·),v(·)) = N−1 ∑ i=0 l(zi,vi) (2.3) where l : Rn×Rm → R is the stage cost. There are mainly three ingredients in MPC: 1) the internal model, 2) the ob- jective functional, and 3) the constraints. A basic categorization of different MPC strategies can be obtained via the selection of the internal model and objective function in their formulations. An MPC formulation is, in general, referred to as linear or nonlinear if the model used is linear or nonlinear. With the type of stage 11 2.4. Nominal Stability of NMPC cost, for example setpoint tracking or economic, MPC techniques used in the linear and nonlinear MPC frameworks are different. Linear MPC techniques are widely popular and have been successfully applied in industrial processes [19]. The receding horizon concept, one of the key ideas of MPC, was first conceptu- alized in [20] to apply linear programming techniques in optimization-based con- trol problems. The description of perhaps the earliest application of MPC, Model Predictive Heuristic Control (MPHC), was presented in [21]. In the MPHC for- mulation, an impulse response model was used with a quadratic stage cost, and input and output constraints. Within the family of linear MPC strategies, the Dy- namic Matrix Control (DMC) [22] and the Generalized Predictive Control (GPC) [23, 24] were key formulations in the early stage of MPC development in both academia and industry, and have been successfully implemented in many indus- trial processes. The DMC formulation uses a linear step-response model while the GPC formulation uses a Controlled Auto-Regressive Integrated Moving-Average (CARIMA) model for setpoint tracking. Thorough reviews of MPC techniques are provided in [19, 25]. The theory of linear MPC is quite mature at present [26]. In this chapter, the setpoint tracking NMPC formulation is reviewed for its nominal and robust stability. 2.4 Nominal Stability of NMPC For the setpoint tracking NMPC, a quadratic stage cost is generally used in the cost functional, i.e. l(z,v) = ‖z− xs‖2Q +‖v−us‖2R (2.4) where xs ∈ Rn and us ∈ Rm are state and control setpoint vectors, respectively. Although the use of the reference control input vector us is not necessary if there is no degree of freedom in the controller i.e. n = m, this formulation might be useful when changing the operating conditions especially for nonlinear processes [27]. Q ∈ Rn×n > 0 and R ∈ Rm×m > 0 are positive definite weighting matrices for state and control vectors, respectively. To simplify the formulation of expressions in this chapter, the setpoint vectors are considered, without loss of generality, as xs = 0 12 2.4. Nominal Stability of NMPC and us = 0. The following definitions and lemmas have been used throughout the chapter. Definition 2.4.1 [28, 29] Let UN(x) be the set of feasible control moves which satisfy the constraints (2.2b)-(2.2d) with the initial state x. Then the admissible set ZN can be defined as, ZN = {(x,v)|x ∈ X and v ∈ UN(x)}, and the set of admissible states XN ⊆ X can be defined as, XN = {x|UN(x) 6= /0}. (2.5) Definition 2.4.2 [30] Consider the system xk+1 = f (xk,uk) with the control law κ :Rn →Rm defined by u= κ(x), and f (0,0) = 0. Then for the closed-loop system, xk+1 = f (xk,κ(xk)) = h(xk) (2.6) the set Xpi is said to be a positively invariant set if h(xk) ∈ Xpi, ∀ xk ∈ Xpi, ∀ k ≥ 0. Further, if Xpi ⊆ XN , then Xpi is called output admissible set which is a domain of attraction of the origin [31, 32]. Definition 2.4.3 [30] A continuous function α : [0,a)→ [0,∞), where a ∈ R≥0, is a class K function if it is strictly increasing and α(0) = 0. It is a class K∞ function if α(a)→ ∞ as a→ ∞. Definition 2.4.4 [30] A continuous function β : [0,a)× [0,∞)→ [0,∞) is a class K L function if, for each fixed s∈R≥0, β (·,s) is a class K function, and for each r ∈ R≥0, β (r,s)→ 0 as s→ ∞. 13 2.4. Nominal Stability of NMPC Definition 2.4.5 [30] A function V (·) is called a Lyapunov function for the system (2.6) if there exists a set M and K∞ functions αv1, αv2, and αv3 such that αv1(x) ≤V (x) ≤ αv2(x), V (h(x))−V (x) ≤−αv3(x) for all x ∈M. (2.7) Lemma 2.4.1 [32, 33] Let V (·) be a Lyapunov function for the system (2.6), and V (·) is defined in an admissible positive invariant set Xpi containing the origin in its interior. Then the origin is asymptotically stable. Further, if αv1 = av1‖x‖bv , αv2 = av2‖x‖bv , αv3 = av3‖x‖bv , for real positive av1, av2, av3, and bv then the origin is an exponentially stable equilibrium point. The constrained optimization problem is formulated for the setpoint tracking NMPC problem at sample time k to achieve the nominal stability of the closed-loop system as, VN(x) = min v N−1 ∑ i=0 l(zi,vi)+F(zN) (2.8a) s. t. : zi+1 = f (zi,vi), z0 = xk, (2.8b) zi ∈ X, vi ∈ U, (2.8c) zN ∈ XF , i = 0, . . . ,N−1, (2.8d) where F(·) is a terminal cost (or called terminal penalty term). XF ⊆ X is a termi- nal set (an output admissible set). The above constrained optimization problem is denoted as PN(k), and the solution of PN(k) is denoted by v∗ = {v∗0,v∗1, . . . ,v∗N−1}. Since the controller is implemented according to the receding horizon fashion, the MPC control law hmpc : Rn → Rm can be defined by hmpc(xk) = v∗0. (2.9) Then the closed-loop system can be described as xk+1 = f (xk,hmpc(xk)). (2.10) 14 2.4. Nominal Stability of NMPC In the following, some fundamental properties of the nonlinear system function f , stage cost l, terminal cost F , constraint sets U and X, and terminal set XF are assumed to show the nominal stability of the closed-loop system (2.10). Assumption 2.4.1 f is continuous, and f (0,0) = 0. Thus the origin is an equilib- rium point of the system (2.1a). Assumption 2.4.2 l is Lipschitz continuous with Lipschitz constant Ll . l(0,0) = 0, and αl(‖x‖) ≤ l(x,u) ≤ βl(‖x‖) where αl and βl are K functions. Assumption 2.4.3 U is compact, and X is closed. The origin (0,0) is an interior point of X×U. Assumption 2.4.4 XF ⊆ XN is compact and contains the origin in its interior, i.e. 0 ∈ XF . Assumption 2.4.5 F is Lipschitz in XF with Lipschitz constant LF . F(0) = 0, and αF(‖x‖) ≤ F(x) ≤ βF(‖x‖) where αF and βF are K functions. Assumption 2.4.6 For a given auxiliary (or called local) continuous control law hF : Rn → Rm, the following expressions hold: 1. hF is Lipschitz continuous in XF with Lipschitz constant Lh f , 2. hF(XF) ∈ U , 3. f (x,hF (x)) ∈ XF , and 4. F( f (x,hF (x)))−F(x)≤−l(x,hF (x)) for all x ∈ XF . The Assumption 2.4.4, and the conditions 2 and 3 in Assumption 2.4.6 imply XF is an output admissible set for the closed-loop f (x,hF (x)) system under the local control law hF . The condition 4 in Assumption 2.4.6 implies that F is de- creasing along the trajectory of closed-loop f (x,hF (x)), which is required to prove the closed-loop stability. The following theorem expresses the nominal stability of the closed-loop sys- tem (2.10). 15 2.4. Nominal Stability of NMPC Theorem 2.4.1 [28, 32] Let Assumptions 2.4.1-2.4.6 hold. Then the origin is an asymptotically stable equilibrium point for the closed-loop system, xk+1 = f (xk,hmpc(xk)), with a region of attraction XN . Further, if αl = al‖x‖bl , βF = aF‖x‖bl , for real positive al , aF , and bl , then the origin is an exponentially stable equilibrium point in XN . The proof of the theorem can be found in [28, 32] where VN(x), the value function, is shown as an Lyapunov function for the closed-loop (2.10) by showing that αl(‖x‖) ≤VN(x)≤ βF(‖x‖) (2.11a) VN( f (x,hmpc(x)))−VN(x)≤−αl(‖x‖) ∀ x ∈ XN . (2.11b) In the view of Lemma (2.4.1), asymptotic and exponential stability can be proven. For guaranteeing closed-loop stability, there are three main ingredients in the NMPC formulation (2.8) : 1) the terminal cost F , 2) the terminal set XF , and 3) the local control law hF . Depending on the different selection of these three F , XF , and hF , many techniques have been proposed to ensure the stability [34–40]. The formulation in [34] was first to propose the use of an equality constraint on the terminal state (zN = 0, i.e. XF = {0} in (2.8d)) to guarantee the stability of the closed-loop system, yet it faces some difficulties on implementation due to the requirement of exact satisfaction of nonlinear constraint, which is difficult and computationally expensive. The exact equality terminal constraint was therefore relaxed in a later development called the dual-mode MPC in [36]. In [36], the terminal state zN is forced to stay inside some neighborhood of the origin XF, i.e. zN ∈ XF. The basic idea of the dual-mode MPC is to apply the local stabilizing linear controller hF inside the terminal set XF and the MPC law hmpc outside the terminal set XF, to guarantee the stability of the closed-loop system. In another de- velopment [37], the terminal cost F was proposed to add to the objective functional (2.8a) to ensure the closed-loop stability without the need of terminal constraint (2.8d). In this formulation, the finite horizon objective functional with the terminal cost F approximates the infinite horizon objective functional. Later, different ideas 16 2.5. Robust Stability of NMPC were proposed in the literature with the combination of all three components: the terminal cost F , the terminal set XF , and the local control law hF . A key develop- ment called quasi-infinite horizon MPC was proposed in [39]. In the quasi-infinite horizon MPC technique, both the terminal constraint (2.8d) and the terminal cost F were employed to approximate the infinite horizon cost and to ensure the stability, respectively, yet the local controller hF was not applied inside the terminal set XF . Over the years, the classical Lyapunov stability theory has been applied to prove the stability of the closed-loop system with the MPC control law. The stan- dard approach is to consider the value function as a candidate Lyapunov function for the closed-loop system [34], and to prove the stability, assuming that both the nonlinear system and value functions are continuous (see Assumptions 2.4.1 and 2.4.2 ). The requirements of the nonlinear system and value function for being continuous were questioned in [41]. In [41], stability was proven without requiring that the MPC control law for being continuous, but by assuming the continuity of the nonlinear system function on a neighborhood of the origin. It was shown [42] that the MPC technique sometimes need to produce a discontinuous control law in order to ensure the closed-loop stability. 2.5 Robust Stability of NMPC In this section, the inherent robustness of the NMPC controller presented in sec- tion 2.4 is analyzed using the Input-to-State Stability (ISS) framework. The ISS definition [32, 43, 44] can be extended to analyze the stability of the system. 2.5.1 Inherent Robustness of NMPC Consider the nominal system f (2.1a) with uncertainty as xk+1 = f (xk,uk)+g(xk,uk,wk), := fcl(xk,uk,wk), (2.12a) xk ∈ X, uk ∈ U, wk ∈W, k ≥ 0, (2.12b) 17 2.5. Robust Stability of NMPC where wk ∈ Rd is the disturbance vector, and W ⊆ Rd is the compact uncertain region. g : Rn×Rm×Rd → Rn is the uncertain part of the system f and assumed to be Lipschitz continuous as follows. Assumption 2.5.1 The uncertain term g(·, ·, ·) is Lipschitz continuous with respect to all its arguments with Lipschitz constant Lg≥ 0 such that for all (z1,v1),(z2,v2)∈ Z, and w1,w2 ∈W, i.e. ‖g(z1,v1,w1)−g(z2,v2,w2)‖ ≤ Lg‖(z1,v1,w1)− (z2,v2,w2)‖. (2.13) This system (2.12) is considered to be controlled by the NMPC control law (2.9). Then the closed loop can be written as xk+1 = fcl(xk,hmpc(xk),wk). (2.14) For robust stability analysis, the following definitions, assumptions, and lemmas, which are common in the NMPC stability analysis, have been used in this section. Definition 2.5.1 [43, 44] A set Xrpi ∈XN is an admissible robust positive invariant (RPI) set for the system (2.12) if fcl(xk,κ(x),wk) ∈ Xrpi, ∀ xk ∈ Xrpi, ∀ wk ∈W, and k ≥ 0. (2.15) with the control law κ : Rn → Rm defined by u = κ(x). Correspondingly, ZRP = {(x,v)|x ∈Xrpi, and w ∈W.} Definition 2.5.2 [43, 44] The closed-loop system fcl is said to be ISS if there exists a K L function γ(·, ·), and a K function δ (·) such that, ‖xk‖ ≤ γ(x0,k)+δ (‖w‖), ∀k ≥ 0, ∀w ∈W, ∀x0 ∈Xrpi (2.16) Furthermore, an ISS-Lyapunov function is defined for the closed-loop system as follows. Definition 2.5.3 [44] A function V (·) is called an ISS-Lyapunov for the system (2.12a) if there exists a set N containing the origin in its interior, K∞ functions 18 2.5. Robust Stability of NMPC αv1(·),αv2(·),αv3(·), and a K function δv1(·) such that for all x∈N and all w∈W, V (x) ≥ αv1(‖x‖) V (x) ≤ αv2(‖x‖) V ( fcl(xk,κ(xk),wk))−V (x,w) ≤−αv3(‖x‖)+δv1(‖w‖). (2.17) Lemma 2.5.1 (see [32, 43]) Let Xrpi be an admissible RPI set for the system (2.12), containing the origin, and V (·) be an ISS-Lyapunov function defined in Xrpi for the system (2.12), then the resulting closed-loop system is ISS in Xrpi. The following theorem expresses the ISS property of the closed-loop system (2.14). Theorem 2.5.1 Let Assumptions 2.4.1 - 2.5.1 hold, then the closed-loop system (2.14) is ISS in an admissible RPI set Xrpi for any perturbation g(·, ·, ·) such that ‖g(x,uk,0)‖ ≤ r Lv αl(‖x‖) (2.18) where r ∈ R and 0 < r < 1 is an arbitrary real number. Lv is Lipschitz constant of V (·)N in (2.11). The proof of the theorem can be found in [32] where it is shown that VN( fcl(xk,hmpc(xk),wk))−VN(x,w) ≤−(1− r)αl(‖xk‖)+LvLg‖wk‖ (2.19) Since the closed-loop system has nominal stability, the condition (2.11b) still holds. Therefore, according to Definition 2.5.3, VN(·) is an ISS Lyapunov function for the closed-loop system. Then it can be concluded from Lemma 2.5.1 that the closed- loop system (2.14) is ISS at the origin. 2.5.2 Robust Design of NMPC In addition to the inherent robustness of the closed-loop system, provided by the nominal NMPC framework, the designing of the NMPC techniques to ensure ro- bust closed-loop stability of uncertain systems has gained much attention among 19 2.6. Solution of NMPC Optimization Problem researchers. These robust NMPC controller designing techniques should ensure 1) recursive feasibility, 2) robust constraint satisfaction, and 3) some closed-loop performance [44, 45]. In min-max approach, a H∞ objective function is first max- imized with respect to disturbances, and then minimized with respect to control moves over the horizon [46, 47]. Although it seems that the min-max approach is the general way to handle uncertain systems, they are computationally demand- ing. To overcome the computational complexity of min-max approaches, tight- ening constraints approaches are alternative methods where the robust stability is guaranteed by tightening the constraint of the nominal optimization problem for any admissible disturbances [48]. However, this approach can create to small con- straint sets or even lead to infeasible problems [32] due to the exponential increase of the uncertainty margin with the increase of the prediction horizon. The tube based NMPC approach, tube-NMPC, is also a constraint tightening method. Com- pared to the min-max approaches, this method provides stronger robust stability and is computationally attractive [49, 50]. In tube-NMPC, a feedback auxiliary controller is utilized with the nominal controller in order to force the trajectory of the closed-loop uncertain system to be inside a tube centered with nominal tra- jectory. The discussion of robust NMPC designing techniques can be found in [51, 52]. 2.6 Solution of NMPC Optimization Problem One drawback of the NMPC technique is that it needs to repeatedly solve a non- linear optimization problem at each sample time. Depending on the size of the problem, this could take a considerable amount of time without guarantee of a so- lution within a specified sample time. Therefore, the application of NMPC was restricted to slow processes at early stages of its development. NLP solvers can exploit a special structure of the NMPC problem to improve the computational performance. The NMPC problem has a special structure due to two main fea- tures: the sparsity of the problem and the availability of good initial solutions [53], [54]. One way of exploiting the block structure is to consider the system equations (output or/and state) as equality constraints in the NLP. Then the Hessian of the ob- jective function and the Jacobian of constraints can be made sparse. In the absence 20 2.6. Solution of NMPC Optimization Problem of disturbance and setpoint changes, a feasible solution can be found using the so- lution of the problem at previous time samples, and this feasible solution can be considered as an initial solution to the current problem [54], [55]. In recent years, researchers have developed fast MPC algorithms such as explicit NMPC [56] and Newton-type controllers [57]. While the explicit NMPC algorithms are suitable for systems with few states, the on-line NMPC algorithms such as Newton-type controllers are more efficient alternatives for large systems [58]. In this work, The IPOPT solver presented in [59] is used to solve the resulting NLP problem of the NMPC formulation (2.8), and the computational burden is further improved with the asNMPC concept presented in [58]. 2.6.1 IPOPT IPOPT uses a primal-dual interior-point algorithm, with a filter line-search strategy, to solve the general form of NLP problem as min w f (w) (2.20a) s. t. : c(w) = 0, (2.20b) w≥ 0, (2.20c) where w ∈ Rn is variable vector. f : Rn → R and c : Rnc → Rmc , with mc ≤ nc, are the objective function and the equality constraints, which are assumed to be twice continuously differentiable. In the context of NMPC formulation, the resulting NLP problem PN(k) (2.8) can be recast in the form of min w f (w,θ) (2.21a) s. t. : c(w,θ) = 0, (2.21b) w≥ 0, (2.21c) where θ is a parameter vector, and θ ≡ xk where xk is the initial state variable vector at kth sample time. Note that the sample time dependence is omitted in the equations. IPOPT computes a series of approximate solutions to a sequence of 21 2.6. Solution of NMPC Optimization Problem barrier problems min w f (w,θ)−δ nc∑ j=1 ln(w j)) (2.22a) s. t. : c(w,θ) = 0, (2.22b) with a decreasing sequence of barrier parameters, δ → 0. The barrier terms, ln(w j) where w j denotes the jth component of the vector w, are added to the objective function in order to handle the inequality constraints (2.21c). With a fixed value of δ , the IPOPT first computes an approximate solution to the barrier problem (2.22) by solving the Karush-Kuhn-Tucker (KKT) conditions: ∇w f (w,θ)+∇wc(w,θ)λ −η = 0 (2.23a) c(w,θ) = 0 (2.23b) ENe−δe = 0 (2.23c) where λ ∈ Rmn and η ∈ Rnc are Lagrange multipliers for the equality constraints and the bound constraints, respectively. Here E = diag(w) and N = diag(η) are diagonal matrices, and e = [1,1, · · · ,1]T is vector of ones. Then by decreasing the barrier parameter, a new barrier problem is solved with the approximate solution of the previous problem. IPOPT solves each barrier problem using Newton’s method for nonlinear equations applied to the KKT conditions (2.23). Let us consider a fixed value of δk, and an iteration i, (wi(θ),λi(θ),ηi(θ)), of the Newton step. Then the search directions (∆wi,∆λi,∆ηi) are obtained by solving following system of linear equations.   Hi Ai −I Ati 0 0 Ni 0 Ei     ∆wi ∆λi ∆ηi  −   ∇w f (wi,θ)+Aiλi c(wi,θ) EiNie−δke  = 0, (2.24) where the matrix Hi = ∇2wwL(wi,λi,ηi) is the Hessian of the Lagrange function L = f (w,θ)+ c(w,θ)T λ −ηT w, and the matrix Ai = ∇c(wi,θ) is the Jacobian of the constraint (2.21c). I and 0 denote identity and zero matrices with appropriate 22 2.6. Solution of NMPC Optimization Problem dimensions. Let us denote s∗k = [w ∗(θ)T ,λ ∗(θ)T η∗(θ)T ], (2.25) the solution for the linear system when δ = 0 at a fixed parameter value θ . The above description is a basic formulation of the IPOPT solver, and more details can be found in [59]. 2.6.2 asNMPC The asNMPC algorithm consists of two main steps: 1) solving a nominal NLP problem in advance, and 2) updating the nominal solution with new information. Suppose at sample time k the control action uk is known. Then the next state zk+1 of the process can be predicted through a forward simulation of the process model (2.8b), and solve the predicted nominal problem PN(k+ 1|k) with θ = θ0 = zk+1. With this strategy, the control policy for k+ 1 sample time, v∗(k+ 1|k), is readily available when time reaches k+ 1, provided that the predicted nominal problem P(k + 1|k) can be solved within a sample time. In the presence of disturbances there is a mismatch between the predicted state zk+1 and actual state xk+1. The asNMPC algorithm then updates the nominal solution s∗k+1|k = [w ∗(θ0)T ,λ ∗(θ0)T η∗(θ0)T ], (2.26) by solving the perturbed KKT system with θ = xk+1 K∗(θ0)∆s(θ) =−ϕ(s∗k+1|k,θ). (2.27) The above equation is closely related to (2.24) where K∗(θ0) is the KKT matrix evaluated at θ0 and s∗k+1|k. The vector ϕ(s∗k+1|k,θ) has the same form as in (2.24) but is evaluated at s∗k+1|k and θ = xk+1. With the assumption that the active-set at the nominal solution is as for the perturbed system, the factorization of the KKT matrix is readily available for the above system from (2.24). As a result, a fast approximate solution can be computed to the problem with a single back-solve. This section presents only a basic description of the method. Detailed discussion 23 2.7. Conclusion about optimality, stability and robustness of asNMPC formulation can be found in [58, 60]. 2.7 Conclusion This chapter provides the basics of NMPC formulation, establishing a foundation for processes that will be discussed in the following chapters. First, the basic idea of receding horizon control is depicted. Next, a mathematical representation of NMPC formulation is presented, along with their basic assumptions. Finally, nom- inal and robust stability of NMPC concepts are introduced. In particular, inherent robustness of the nominal NMPC technique and designing robust NMPC methods are discussed for uncertain nonlinear systems. The basic formulation of the IPOPT solver is presented to solve the resulting NLP problem of NMPC. Using the as- NMPC concept, a fast approximation of the NLP problem is sought with aims to reduce the on-line computational time. 24 Chapter 3 Economically Oriented NMPC 3.1 Summary In this chapter, an economically oriented process control technique is presented. The NMPC framework discussed in Chapter 2 is utilized to explicitly consider an economical objective function in the formulation. The nominal and robust stability properties are discussed for a nonlinear system with the infinite horizon econN- MPC technique. 3.2 Introduction The economical operation of a process is usually conducted using an advanced control system with a two-layer structure. The conventional two-layer structure consists of a higher layer to implement RTO, followed by a lower layer to track the setpoint. The RTO layer computes optimal steady-state setpoints using a steady- state process model, and provides them to the setpoint tracking layer in order to maximize an economical objective while keeping product quality constraints within safe values. The setpoint tracking layer uses advanced control strategies, most commonly NMPC techniques. The conventional two layer structure has some drawbacks mainly due to 1) low frequency of operation of the RTO layer, and 2) assumption of slow variation of the disturbances on the process [15]. Due to the use of nonlinear steady-state modeling of the process in the RTO layer, the re- sulting constrained optimization problem can take considerable amount of time. Therefore, it might not be able to operate simultaneously with the NMPC dynamic optimization layer. The slow variation of disturbances during setpoint update is also not a valid assumption in many industrial processes. Literature on this topic generally discusses few alternatives to this two lay- 25 3.2. Introduction ers approach. The basic idea is to recalculate or modify the setpoints using less- detailed steady-state process models, most commonly linear models within the lower-level [61]. This modification can be performed with a separate optimization within the setpoint tracking layer or can be integrated with the dynamic optimiza- tion problem in NMPC formulation [62]. The main issues with the application of NMPC with detailed nonlinear models are due to computational burden and lack of theoretical support for stability [63, 64]. However, with the increase of compu- tational power and development of efficient Nonlinear Programing (NLP) solvers, the economical objective and setpoint tracking objective tasks could be simulta- neously implemented with a single layer, with detailed nonlinear models in the advanced control system [65–68]. With an application of economical objective function of specific energy, we have recently demonstrated a potential reduction of about 13% in total specific energy for a two-stage TMP process in a simula- tion study [67]. In another simulation study [69], a cost reduction of about 30% has been realized for an air separation unit by taking into account the periodically changing electricity price. Practical economical benefits are discussed in [70–72] for optimal operation of processes. Although the stability theories of these eco- nomically oriented NMPC techniques are still in early stages, there are interesting developments in the literature [29, 69, 73–76]. Nominal stability properties of two economically oriented NMPC schemes for cyclic processes are presented in [69], and robust stability of economically oriented infinite horizon NMPC is presented in [74]. A discount factor, which projects future profit to the present value, is introduced in [77] for NMPC formulation to maximize the net present value of the economical objective function. In our recent paper [69], nominal stability of infinite horizon NMPC with a discount factor for cyclic processes was proved by introducing a transformed system. In this chapter, a general infinite horizon NMPC technique is discussed with an economic objective, and more work on this topic can be found in [69, 74]. 26 3.3. Infinite Horizon economically oriented NMPC (econNMPC) 3.3 Infinite Horizon economically oriented NMPC (econNMPC) Consider a discrete time nonlinear system described by xk+1 = f (xk,uk), k ≥ 0, (3.1a) subject to the control and state constraints, uk ∈U and xk ∈X, (3.1b) where x ∈ Rn and u ∈ Rm are the state and control variable vectors of the system, respectively. U⊆Rm is the control constraint set, and X⊆Rn is the state constraint set. f :Rn×Rm→Rn is the nominal nonlinear system function, and k is the sample time. In a sequel, z and v are denoted as the predicted state and control vectors, respectively, by the model f . In the econNMPC framework, the following performance index J : Rn×Rm → R 1 is considered with a general economic stage cost functional l : Rn×Rm → R1, for the system (3.1a) - (3.1b). J(x(·),u(·)) = limN→∞ 1 N N ∑ k=0 l(xk,uk). (3.2) The objective of the infinite horizon optimal control problem at the kth sample time is to find a control law uk = h(xk) where h : Rn → Rm, by minimizing the performance index J subject to the nonlinear dynamic system (3.1a), and the state and control constraints (3.1b), such that the closed-loop system (3.3) is feasible. xk+1 = f (xk,h(xk)). (3.3) This infinite horizon optimal control problem is referred as P∞, assuming that there exists a unique optimum steady-state (x∗,u∗) such that l(x∗,u∗) < l(xk,uk), (xk,uk) 6= (x ∗,u∗), for all uk ∈ U, xk ∈ X. Further, the following assumptions are made in this chapter, as discussed in [29, 69]. 27 3.3. Infinite Horizon economically oriented NMPC (econNMPC) Assumption 3.3.1 f (·, ·) and l(·, ·) are Lipschitz continuous on the admissible set ZN (see Definition 2.4.1. Here, ZN is defined accordingly for the problem P∞), with Lipschitz constants L f and Ll ≥ 0 such that for all (z1,v1), (z2,v2) ∈ ZN , i.e. ‖ f (z1,v1)− f (z2,v2)‖ ≤ L f‖(z1,v1)− (z2,v2)‖ ‖l(z1,v1)− l(z2,v2)‖ ≤ Ll‖(z1,v1)− (z2,v2)‖ (3.4) Assumption 3.3.2 There exists a K∞ function γ(·) such that for every state x∈XN , where XN is the admissible state set (see Definition 2.4.1. Here, XN is defined accordingly for the problem P∞), there exists u such that (x,u) ∈ ZN and N−1 ∑ k=0 ‖uk−u ∗‖ ≤ γ(‖x− x∗‖), (3.5) and (xk,uk) = (x∗,u∗) for all k≥ N. The Assumption 3.3.2 means that some degree of system controllability is needed [29, 34]; starting from an initial state x, the system can be steered to the steady-state in N steps with a bounded input sequence, while satisfying the input and state constraints. In the following section, an econNMPC is presented with a modified stage cost. The econNMPC with a Regularization Term In this section, an econNMPC problem is formulated with a regularization term in the economic stage as min (zi,vi) ∞ ∑ i=0 l(zi,vi)+‖zi− zi+1‖2Q +‖vi− vi+1‖2R (3.6a) s.t. : zi+1 = f (zi,vi) (3.6b) z0 = xk (3.6c) zi ∈ X, and vi ∈ U, i = 1, . . . ,∞. (3.6d) The regularization terms ‖zi − zi+1‖2Q + ‖vi− vi+1‖2R, where Q and R are positive definite weighting matrices, force the system to converge to a steady-state. And also, it ensures the optimization problem (3.6) is well-posed and the solution is 28 3.4. Nominal Stability of the econNMPC locally unique. The discussion on strong convexity requirement for the stage cost can be found in our recent paper [74]. The following assumption states the well- posedness of the problem (3.6) [78]. Assumption 3.3.3 The constrained optimization problem (3.6) satisfies the Linear Independent Constraint Qualification (LICQ), Sufficient Second order Condition (SSOC), and Strict Complementarity (SC) at the solution. 3.4 Nominal Stability of the econNMPC The nominal stability of the econNMPC is analyzed in this section. The notion of transformed system [69] is considered and nominal stability is proved using Lyapunov stability theory. The transformed system is defined by subtracting the optimal steady-state from the original system (3.1a)-(3.1b) as z̄i+1 = f (z̄i + x∗, v̄i +u∗)− x∗ := ¯f (z̄i, v̄i), (3.7) v̄i ∈ ¯U and z̄i ∈ ¯X, (3.8) where v̄i = vi−u∗ and z̄i = zi− x∗. ¯f : Rn×Rm → Rn is the nominal transformed nonlinear system function. ¯U⊆Rm and ¯X⊆Rn are corresponding constrained sets for the transformed variables z̄ and v̄, respectively. Here it can be seen that when (z̄i, v̄i)→ (0,0), then (zi,vi)→ (x∗,u∗), i.e. the original system will converge to the optimal steady-state (x∗,u∗) if the transformed system is steered from an initial state to 0. The nominal stability of the original system is proved using Lemma (3.4.1). Lemma 3.4.1 (see [69]) The stability of the transformed system (3.7) at (0,0) is equivalent to the stability of the original closed-loop system (3.1a) at the steady- state (x∗,u∗). 29 3.4. Nominal Stability of the econNMPC Further, the original objective function (3.6a) is modified to have the trans- formed objective as V∞(x̄k) = ∞ ∑ i=0 Li(z̄i, v̄i) := ∞ ∑ i=0 l(z̄i + x∗, v̄i +u∗)− l(x∗,u∗)+ ‖(zi− x ∗)− (zi+1− x ∗)‖2Q +‖(vi−u ∗)− (vi+1−u ∗)‖2R = ∞ ∑ i=0 ¯l(z̄i, v̄i)+‖(zi− x∗)− (zi+1− x∗)‖2Q +‖(vi−u∗)− (vi+1−u∗)‖2R, (3.9) where x̄k = xk − x∗ at sample time k, and ¯l = l(z̄i + x∗, v̄i +u∗)− l(x∗,u∗). Similar to the original system, the following assumptions are made on Lipschitz continuity and weak controllability for the transformed system. Assumption 3.4.1 ¯f (·, ·) and Li(·, ·) are Lipschitz continuous on the admissible set of the transformed system, ¯ZN , with Lipschitz constants L ¯f and LL ≥ 0 such that for all (z1,v1), (z2,v2) ∈ ¯ZN , i.e. ‖ ¯f (z1,v1)− ¯f (z2,v2)‖ ≤ L ¯f‖(z1,v1)− (z2,v2)‖ ‖Li(z1,v1)−Li(z2,v2)‖ ≤ LL‖(z1,v1)− (z2,v2)‖. (3.10) Assumption 3.4.2 There exists a K∞ function γ∞(·) such that for every state x̄ ∈ ¯XN , there exists ū such that (x̄, ū) ∈ ¯ZN and N−1 ∑ k=0 ‖ūk‖ ≤ γ∞(‖x̄‖), (3.11) where ū is a feasible input sequence, and ¯XN is the admissible state set for the transformed system, respectively. The validity of Assumptions 3.4.1 and 3.4.2 can be easily seen in the view of Assumptions 3.3.1 and 3.3.2, respectively. The Assumption 3.4.2 means that starting from an initial state x̄, the transformed system can be steered to (0,0) 30 3.4. Nominal Stability of the econNMPC in N steps with a bounded input sequence, while satisfying the input and state constraints of the transformed system. For this particular feasible solution (ẑi, v̂i), the expression (3.12) holds. Li(ẑi, v̂i) = 0 i≥ N (3.12) Furthermore, Assumption 3.4.3 is made on the transformed stage cost Li(·, ·). Assumption 3.4.3 There exits a K∞ function β∞(·) such that the stage cost Li(·, ·) defined in (3.9), satisfies Li(z̄, v̄)≥ β∞(‖z̄−0‖). (3.13) Note that Assumption 3.4.3 may not hold for general economic objective functions. An alternative approach, a modification of Li(z̄, v̄), is proposed using rotated stage cost in [29] while a generalization of that is developed in our recent paper [74]. The constrained optimization problem in the transformed econNMPC formu- lation is expressed as min (z̄i,v̄i) ∞ ∑ i=0 Li(z̄, v̄) (3.14a) s.t : z̄i+1 = ¯f (zi,vi) (3.14b) z̄0 = xk− x ∗ (3.14c) z̄i ∈ ¯X, and v̄i ∈ ¯U, i = 1, . . . ,∞. (3.14d) The following theorem expresses asymptotic stability of the nominal system at the optimal steady-state (x∗,u∗) in the view of asymptotic stability of the trans- formed system at origin (0,0). Theorem 3.4.1 Let Assumptions 3.3.1 - 3.4.3 hold, then at time step k, the value function V∞(x̄k) defined by (3.9) is a Lyapunov function, and the transformed sys- tem (3.7) is asymptotically stable at (0,0) with the controller based on (3.14). Consequently, the original closed-loop system controlled by the econNMPC (3.6) is asymptotically stable at the optimal steady-state (x∗,u∗). 31 3.5. Robust Stability of the econNMPC The proof of the Theorem 3.4.1 is given in [69]. That V∞(x̄k) is a Lyapunov function for the transformed system is proven in [69], by showing β∞(‖x̄k‖)≤V∞(x̄k)≤ α∞(‖x̄k‖), (3.15) V∞( ¯f (x̄k, ūk))−V∞(x̄k)≤−β∞(‖x̄k‖), (3.16) where α∞(·) is a K∞ function. In the view of Lemma (3.4.1), the proof is com- pleted for the original system. Further, V∞(·) is Lipschitz continuous with Lipschitz constant Lv ≥ 0. 3.5 Robust Stability of the econNMPC In this section, the inherent robustness of the econNMPC controller presented in the previous section is analyzed using the ISS framework. Consider the nominal system f with uncertainty as xk+1 = f (xk,uk)+g(xk,uk,wk) := h(xk,uk,wk), xk ∈X, uk ∈U, wk ∈W, (3.17) where wk ∈ Rd is the disturbance vector, and W ⊆ Rd is the compact uncertain region. g : Rn ×Rm ×Rd → Rn is the uncertain part of the system f , and the Lipschitz continuity of g is assumed with respect to all its arguments, as follows. Assumption 3.5.1 The uncertain term g(·, ·, ·) is Lipschitz continuous with respect to all its arguments with Lipschitz constant Lg ≥ 0 such that ∀ (z1,v1),(z2,v2) ∈ Z, and w1,w2 ∈W, i.e. ‖g(z1,v1,w1)−g(z2,v2,w2)‖ ≤ Lg‖(z1,v1,w1)− (z2,v2,w2)‖. (3.18) 32 3.5. Robust Stability of the econNMPC The transformed system can be defined for the system (3.17), similar to the nominal system in (3.7), x̄k+1 = f (x̄k + x∗, ūk +u∗)+g(x̄k + x∗, ūk +u∗,wk) = ¯f (x̄k, ūk)+ ḡ(x̄k, ūk,wk) := ¯h(x̄k, ūk,wk), (3.19) where ḡ : Rn×Rm×Rd → Rn is the uncertainty part of the transformed system ¯f , and ḡ(x̄k, ūk,wk) = g(x̄k + x∗, ūk +u∗,wk). The transformed nominal system ¯f (·, ·) is Lipschitz continuous with the Lipschitz constant L ¯f ≥ 0 (3.4.1). In the view of Assumptions (3.3.1) and (3.5.1), it can be seen that ḡ(·, ·) is also Lipschitz contin- uous with the Lipschitz constant Lḡ for all (z1,v1), (z2,v2) ∈ ZN , and w1,w2 ∈W, such that ‖ḡ(z̄1, v̄1,w1)− ḡ(z̄2, v̄2,w2)‖ ≤ Lḡ‖(z̄1, v̄1,w1)− (z̄2, v̄2,w2)‖. (3.20) For the robust stability analysis of the transformed system (3.19), the following definitions and lemmas which are adopted for the transformed system are used. Definition 3.5.1 (Robust positively invariant set: Transformed system) [32, 43, 44] A set ¯XN is a RPI set for the transformed system (3.19) if ¯h(x̄k, ūk,wk) ∈ ˆXN , for all x̄k ∈ ˆXN , wk ∈W. The ISS definition [32, 43, 44] can be extended to analyze the stability of the transformed system. Definition 3.5.2 (Input-to-state stability: Transformed system) [74] The trans- formed closed-loop system, the system (3.19) controlled by the controller (3.14), is said to be ISS if there exits a K L function γs(·, ·), and a K function δ (·) such that, ‖x̄k‖ ≤ γs(x̄0,k)+δ (‖w‖),∀k ≥ 0, ∀x̄0 ∈ ˆXN . (3.21) 33 3.5. Robust Stability of the econNMPC Further, ISS-Lyapunov function is defined for the transformed system as fol- lows. Definition 3.5.3 (ISS-Lyapunov function: Transformed system) [74] A function V∞(·) is called an ISS-Lyapunov for the transformed system (3.19) if there exists a RPI set ˆXN containing the origin in its interior, K∞ functions α1(·),α2(·),α3(·) and δ1(·) such that ∀x̄ ∈ ˆXN and ∀w ∈W , V∞(x̄)≥ α1(‖x̄‖) V∞(x̄)≤ α2(‖x̄‖) ∆V∞(x̄,w) =V∞(¯h(x̄, ū,w))−V∞(x̄)≤−α3(‖x̄‖)+δ1(‖w‖). (3.22) Lemma 3.5.1 (see [32, 43]) Let ˆXN be a RPI set for the transformed system (3.19), containing the origin, and let V∞(·) be an ISS-Lyapunov function for the trans- formed system (3.19), then the resulting closed-loop system is ISS in ˆXN . Similar to Lemma (3.4.1), the ISS of the original closed-loop system, the sys- tem (3.1a) controlled by the controller (3.6), can be investigated using Lemma (3.5.2). Lemma 3.5.2 (see [74]) The ISS stability of the transformed closed-loop system, the system (3.19) controlled by the controller (3.14), at 0 is equivalent to the ISS stability of the original closed-loop system at the optimal steady-state (x∗,u∗) when x̄ ∈ ˆX, x ∈ X, and w ∈W. The following theorem expresses the robust stability property of the original closed-loop system. Theorem 3.5.1 Let Assumptions 3.3.1 - 3.5.1 hold, then the transformed closed- loop system is ISS in ˆXN for any perturbation ḡ(·, ·, ·) such that ‖ḡ(x̄, ūk),0‖ ≤ r Lv β∞(‖x̄−0‖) (3.23) where r ∈ R and 0 < r < 1 is an arbitrary real number. Then from Lemma 3.5.2, the original closed-loop system is ISS at the optimal steady-state (x∗,u∗). 34 3.6. Conclusion Proof 3.5.2 Since the original system has the nominal stability, the inequalities (3.15) still holds. With the Lipschitz continuity of V∞(·) for the transformed system, it can be seen that V∞(¯h(x̄k, ūk,wk))−V∞(x̄k) =V∞(¯h(x̄k, ūk,wk))−V∞( ¯f (x̄k, ūk))+V∞( ¯f (x̄k, ūk))−V (x̄k) ≤ Lv‖ḡ(x̄k, ūk,wk)‖−β∞(‖x̄k −0‖), from (3.20) and (3.23), V∞(¯h(x̄k, ūk,wk))−V (x̄k)≤ Lv‖ḡ(x̄k, ūk,0)‖+Lvlḡ‖wk‖−β∞(‖x̄k −0‖) ≤ rβ∞(‖x̄k−0‖)+Lvlḡ‖wk‖−β∞(‖x̄k−0‖) ≤−(1− r)β∞(‖x̄k −0‖)+Lvlḡ‖wk‖ (3.24) Hence, according to Definition 3.5.3, V∞(·) is an ISS Lyapunov function for the transformed system. Then the transformed system is ISS at (0,0). In the view of Lemma 3.5.2, the original system is ISS at the optimal steady-state (x∗,u∗) 3.6 Conclusion An econNMPC formulation with a general cost function is presented in this chap- ter. The economical objective function is modified by adding regularization terms. The regularization terms enforce strong convexity of the stage cost and also force the system to a steady-state. The nominal and robust stability of the econNMPC are proven by considering the value function as a Lyapunov function. More general re- sults of the econNMPC can be found in [69, 74], where the nonlinear systems with cyclic behavior have been considered. The analysis in this Chapter will provide the theoretical basis for the following Chapters 5 and 6 where optimal operation of thermo mechanical pulping processes is performed by using the economically oriented NMPC framework. 35 Chapter 4 Closed-loop Identification with Constrained MPC 4.1 Summary In this chapter, two closed-loop identification schemes, which can be applied to a LTI system controlled by a constrained MPC law, are developed. In the first scheme, the informativity of the normal operational closed-loop data is analyzed, and a criterion is presented on the closed-loop data set for applying the direct identification technique with the Prediction Error Method (PEM). In the second scheme, a joint input-output, the two-stage method for closed-loop identification, is discussed. In the first stage of two-stage identification method, the proper com- ponent in the control signal is extracted by using the projection method while in the second stage, a model is identified with the use of the PEM. Furthermore, the Long Range Predictive Identification (LRPI) criterion is utilized in the second stage of identification to obtain a better model estimate. 4.2 Introduction Performance of a closed-loop system with a model based controller mainly de- pends on the internal model used in the controller; when the accuracy of the pro- cess model depreciates, the performance deteriorates greatly [79]. In the case of MPC, the performance mainly depends on the process model to predict the future outputs/states of the process. In practice, the process model is first identified us- ing open-loop identification methods in the commissioning stage, then a stabilizing and well-performing MPC is setup for the process. With time, the performance of 36 4.2. Introduction the closed-loop gracefully degrades mainly due to a different dynamic behavior of the process, caused by wear and tear or other physical changes in the process [80]. Therefore, to keep the closed-loop system performance at the desired level, it may be necessary to carry out the re-commissioning on a periodic basis [81]. How- ever, these re-commissioning stages can be lengthy and costly; the cost of process modeling is around 75% of the total costs associated with advanced process con- trol projects [79, 82]. Moreover, the open-loop test can not be performed during the operation of a process due to the production commitment. To avoid opening the closed-loop and performing the re-commissioning, closed-loop identification should be applied to the system to identify the process model during the normal operation. In general, closed-loop identification techniques can be categorized into three main groups, based on assumptions on the feedback mechanism and the type of the controller [83–85]: 1. Direct identification uses PEM to identify the open-loop system with the input-output data ignoring any effect of feedback, as in the open-loop case. This approach, when it works, is the simplest and the most attractive one, as it does not need any special identification algorithm. 2. Indirect identification consists of two steps: first, identification of either the closed-loop system with the use of measured external independent excitation input (or reference) and output data [86], or sensitivity function with the use of measured input-output data. The open-loop model is then determined using the knowledge of the linear feedback. 3. Joint input-output identification assumes that the input is generated using a certain form of the controller, but the exact knowledge of the controller is not necessary. This approach regards both the control and output signals to- gether as the joint outputs of the independent joint inputs consists of process excitation or reference signal, and the process noise. The joint input-output identification method first extracts the excitation or reference signal compo- nent in the control signal, then uses this extracted signal and the output data to fit a process model. 37 4.3. Output Feedback MPC The indirect identification method heavily depends on linearity of the con- troller though it can be applied to known nonlinear controller with the considerable amount of work [84]. The constrained MPC is usually nonlinear [54]. In this chap- ter, first an informativity criterion for normal operation data, which is generated by MPC, is analyzed for direct closed-loop identification method. Next, the joint input-output methods, particularly the projection method in which the controller can be nonlinear [87], are studied for closed-loop identification. 4.3 Output Feedback MPC In this chapter, output feedback MPC for an LTI discrete-time Single-Input Single- Output (SISO) process (4.1) is considered. yk = G(z)uk +nk, (4.1a) nk = H(z)ek, k ≥ 0, (4.1b) where G(z) and H(z) are the process and noise transfer functions, respectively. uk ∈ R and yk ∈ R are the control input and the process output, respectively. ek ∈ R ∼ N(0,σ 2e ) is a Gaussian white noise signal with the variance σ 2e . Both G(z) and H(z) are rational and proper functions. Further, H(z) is inversely stable. The internal model, a model of the process (4.1), which is used in the MPC formulation is given by ŷk = ˆG(z)ûk + n̂k, (4.2a) n̂k = ˆH(z)êk, (4.2b) where ˆG(z) = z−d B(z)A(z) is the model transfer function. H(z) = C(z) D(z) is the inversely stable model noise filter. The degrees of the polynomials A, B, C, and D are deg(A) = deg(B) = n and deg(C) = deg(D) = m. ŷ ∈ R, û ∈ R and ê ∈ R are used to denote the model predicted output, calculated control input, and model predicted white noise, respectively. The model delay d ≥ 1 in ˆG(z) is enforced to avoid circulating dependence between ŷk and ûk. The dynamic optimization 38 4.3. Output Feedback MPC problem is formulated for the MPC at sample time k as follows. min û J(yk, û) J(yk, û) = N−1 ∑ i=0 (ŷk+i+1− rk+i+1)2 +λ û2k+i, (4.3a) s. t. : Aû≤ b, (4.3b) where N is, for simplicity, the length of the horizons in sample times, without dif- ferentiating the control and the prediction horizons. û = {ûk, ûk+1, . . . , ûk+N−1} ∈ UN(k) where UN(k) is the set of feasible control moves over the horizon N which satisfy the constraints (4.3b). rk ∈R, the reference signal, is the setpoint of the pro- cess output yk. λ ∈ R+ is a positive real number. A ∈ Rnc×N is the constraint ma- trix. nc is the number of linear constraints, and b ∈ Rnc . The solution of the above constrained optimization problem (4.3) is denoted by û∗ = {û∗k , û∗k+1, . . . , û∗k+N−1}. The control input applied to the process is uk = û∗k according to the receding hori- zon technique. In the following, the output prediction is described for the MPC formulation. 4.3.1 Output Prediction The minimum variance predictor for ŷ can be written as follows [54, 88]. ŷk+i = Ei(z)uk+i−d + { Fdi (z) C(z) yk + z−1Fi(z)C(z)− z−dFdi (z)B(z) A(z)C(z) uk } (4.4a) ˆdk+i = Fdi (z) D(z) n̂k, (4.4b) where n̂k = ˆH−1(z)(yk − ˆG(z)uk) n̂k+i = 0, i≥ 1. 39 4.3. Output Feedback MPC The polynomials Ei, Fi, Edi , and Fdi in (4.4) are determined by solving the Dio- phantine equations (4.5). C(z) D(z) = Edi (z)+ z −i Fdi (z) D(z) , (4.5a) B(z) A(z) = Ei(z)+ z−(i−d+1) Fi(z) A(z) (4.5b) The polynomials Ei, Fi, Edi , and Fdi can be calculated via long division such that degrees of polynomials satisfy deg(Ei) ≤ (i− d), deg(Fi) = (n− 1), deg(Edi ) ≤ (i−1), and deg(Fdi ) = (m−1). When the delay d ≥ 1 in (4.4), the predictor ŷk+i can be decomposed into two parts: 1) the first term in (4.4a) is the future process behavior which depends on current and future control actions determined by the the optimization problem (4.3) 2) the second term in (4.4a) is the past process behavior which is a combination of past control input and process output. 4.3.2 Explicit Solution of MPC Similar to the MPC formulation with a state space model [56], the constrained optimization problem in the MPC formulation (4.3) with the input output model (4.2a) can be written as a Quadratic Programming (QP) problem [88]. min û ûT Mû+2fTk Lû, s. t. : Aû≤ b, (4.6) 40 4.3. Output Feedback MPC where M = LL +λ I, fTk = (fyk + fuk − r)T L. With the predictor (4.4), L, fyk, fuk , and r are given in [88] as follows. L =   e0 0 . . . 0 0 . . . 0 e0 e1 . . . 0 0 . . . 0 . . . . . . . . . . . . . . . . . . . . . e0 e1 . . . eN−d 0 . . . 0   , (4.7a) Ei(z) = e0 + e1z+ . . .+ ei−dz−(i−d) (4.7b) fyk = [ Fd1 (z) C(z) , Fd2 (z) C(z) , . . . , FdN (z) C(z) ]T yk, (4.7c) fuk =   zF1(z)C(z)−z−dFd1 (z)B(z) A(z)C(z) . . . zFN(z)C(z)−z−dFdN (z)B(z) A(z)C(z)  uk, (4.7d) r = [rk+1, . . . ,rk+N ] T . (4.7e) The above constrained QP problem (4.6) is explicitly solved using KKT conditions for each set of active constraints. The final explicit control law can be written as in [88], uk =Cp +Fyp (z)yk +Frp(z)rk, (4.8) where Cp is a constant, Fyp is a linear causal filter, and Frp is a linear non-causal transfer function. p = {p1, p2, . . . , pq} ⊆ {1, . . . ,nc} is the index set of active con- straints q. The regional constraint for p is being active constraint set depends on the free response (fyk + f u k) of the transfer function model (4.2a) while it depends on the current state in the state space formulation [56]. Here, it can be noticed that the constrained MPC law is nonlinear, and it can be represented by a piecewise affine controller with each piece given by (4.8). 41 4.4. Prediction Error Method (PEM) 4.4 Prediction Error Method (PEM) Here, a brief description of PEM is given for the parameter estimation problem, and the complete description about PEM including its convergence property can be found in [83, 84]. Consider a particular model in the following model structure M. M(θ) : yk = ¯G(z,θ)uk + ¯H(z,θ)ēk , (4.9) where ¯G and ¯H are the process and noise model, respectively. θ ∈ D ⊂ Rp is the parameter vector with the dimension p, and D is the parameter set. ē is the innovation. The one-step-ahead predictor [84] for (4.9) is given by ȳk|k−1 = ¯H(z,θ)−1 ¯Guk + ( I− ¯H−1(z,θ) ) yk. (4.10) Then the prediction error for the system (4.9) is given by εk(θ) = yk− ȳk|k−1 = ¯H(z,θ)−1 ( yk− ¯G(z,θ)uk ) (4.11) With the quadratic error criterion and a data set ZNd = {y1,u1,y2,u2, . . . ,yNd ,uNd}, the prediction error estimate can be formulated as ¯θ = arg min θ∈D VNd (θ ,ZNd ), (4.12a) VNd (θ ,ZNd ) = 1 Nd Nd∑ k=1 εk(θ)2. (4.12b) 4.5 Informativity of Closed-loop Data: Direct Method In this section, direct closed-loop identification with PEM [83–85] is discussed for applying it on a closed-loop data set generated with the constrained MPC. The direct closed-loop identification uses PEM to identify the process using the closed- loop data set as in the open-loop case ignoring any effect of feedback. A consistent estimate of the true process can be obtained using PEM if 1) the data set ZNd is informative enough with respect to a given model set, and 2) the model set con- 42 4.5. Informativity of Closed-loop Data: Direct Method tains the true process, i.e. ¯G(z,θ0) = G(z) and ¯H(z,θ0) = H(z) for some θ0 ∈D for the true process (4.1) and the model (4.9) [84]. The second condition, a globally identifiable model structure, is generally assumed in order to achieve the desired global solution, and a model structure has to be carefully selected for unbiased es- timation; the disturbance model should be independently parameterized and given enough complexity [84, 87]. The first condition guarantees the uniqueness of the global minimum of PEM with respect to the model structure, and experimental conditions should be applied to ensure that the data set is informative [89]. Since the data set is generated by the true process, these experimental conditions depend on both the true process and the feedback controller for the closed-loop system [89]. A condition, on the closed-loop data set (4.14) generated by the process (4.1) with MPC laws (4.8), for informativity is given in [88]. Here, it is summarized as follows. Consider a closed-loop data set ZNd generated by the process (4.1) ZNd = {y1,u1,y2,u2, . . . ,yNd ,uNd}, (4.13) with s pieces of active MPC laws (4.8) as u j k =C j +F y j (z)y j k +F r j (z)rk, j = 1, . . . ,s, and (uk,yk) ∈ ZNd . (4.14) Then the prediction error (4.11) corresponding to the data set (4.14) for the model (4.9) can be written as ε jk = ˆH −1(z)(y jk− ˆG(z)u j k), j = 1, . . . ,s. (4.15a) Substituting (4.14) in the process (4.1), the prediction error ε jk is given by ε jk = ˆH −1(z) [( 1−G(z)Fyj (z) )−1 G(z)− ˆG(z) ( I +Fyj (z) ( I−G(z)Fyj (z) )−1 G(z) )] [( C j +Frj (z)rk ) + ˆH−1(z) ( 1− ˆG(z)F yj (z) )−1 H(z)nk ] . (4.15b) 43 4.5. Informativity of Closed-loop Data: Direct Method It is shown in [88] that the global asymptotic minimum of (4.12b) is achieved when ε jk = ek, and if the matrix F = ( I . . . I 0 . . . 0 −Fy1 . . . −F y s F r1 . . . Frs ) (4.16a) has full row rank, then ¯G = G, and ¯H = H. (4.16b) is unique solution of ε jk = ek. With the criterion (4.16a), the following algorithm can be expressed to identify the models ¯G and ¯H for the process (4.1). Algorithm 2 The basic steps of the direct closed-loop identification with con- strained MPC 1: Select the model structure M (4.9). 2: Check the condition (4.16a) for the closed-loop data set Znd (4.14). 3: if condition (4.16a) is passed, then 4: estimate the models ¯G and ¯H using PEM (4.12). 5: else 6: Perform on-line experiments to generate an informative data set Znd . 7: end if The design and performance of on-line experiments in order to pass the condi- tion (4.16a) are discussed in [88]. These experiments are merely done to activate some pieces of the affine MPC law through the regional constraints. More details can be found in [90] on the idea of the Algorithm 2 under the topic of least costly identification; first, identify models based on the operation closed-loop data, and if necessary, design and perform the least disturbing experiment in the closed-loop. 44 4.5. Informativity of Closed-loop Data: Direct Method 4.5.1 Example The objective of this example is to study several cases to show the effect of com- plexity of MPC for consistent estimates. Consider the following process, yk = G(z)uk +H(z)ek, (4.17a) G(z) = 0.6z −1 1−1.6z−1 +0.7z−2 , H(z) = 1+0.8z−1 1−1.2z−1 +0.2z−2 , (4.17b) where e ∼ N(0,0.01). For the MPC formulation, the model described in 4.18 is used as its internal model. ŷk = ˆG(z)uk + ˆH(z)êk, (4.18a) ˆG(z) = 0.5z −1 1−0.9z−1 , ˆH(z) = 1+0.5z−1 1− z−1 . (4.18b) Then the constrained optimization problem in the MPC can be written as min û J(yk, û) = N−1 ∑ i=0 (ŷk+i+1− rk+i+1)2 +λ û2k+i, (4.19a) s. t : ûk+i ≤ 1, ∀ i = 0, . . . ,N−1. (4.19b) Then for the prediction horizon N = 2 and the controller weighing λ = 0.01, MPC laws are given as follows. For no being active constraints, u(k) = Fy1 (z)yk +F r 1 (z)rk, (4.20a) Fy1 (z) = −0.754+0.6786z−1 1−0.4054z−1−0.2642z−2 , (4.20b) Fr1 (z) = −3.232z2 +5.027z−0.0395−1.68z−1 1−0.4054z−1−0.2642z−2 . (4.20c) 45 4.5. Informativity of Closed-loop Data: Direct Method For second constraint being active i.e. uk+1 = 1, u(k) = Fy2 (z)yk +F r 2 (z)rk +1.0465, (4.21a) Fy2 (z) = −3.081+2.773z−1 1−1.06z−1−0.00973z−2 , (4.21b) Fr2 (z) = 0.973z2 +0.6919z−0.8703−0.4865z−1 1−1.06z−1−0.00973z−2 . (4.21c) For first constraint being active, uk = 1. (4.22) The closed-loop is simulated in MATLAB, and 2.4× 104 samples were col- lected for experiments. Box-Jenkins (BJ) model structure was selected, and identi- fication task was performed via MATLAB command bj. Here, for simplicity, order of numerator polynomials for both the estimates ˆG and ˆH are assumed the same as the true process 4.17. Following cases are studied. Case r Constraint Control law A r = 0 Unconstrained uk = Fy1 (z)yk B r ∼ N(0,0.01) Unconstrained uk = Fy1 (z)yk +Fr1 (z)rk C r = 0 Constrained uk =   Fy1 (z)yk or Fy2 (z)yk +1.0465 or 1 Table 4.1: Three different case studies and respective control laws for Example 4.5.1 The complexity criterion (4.16a) discussed for direct closed-loop identification are given in Table 4.2. Note that consistent estimates are expected when F has full raw rank 2. For the case A, the rank of F is 1. Therefore, the estimates are quite biased as shown in Fig. 4.1 and Fig. 4.2. For the each case B and C, the rank of F is 2, and the estimates are consistent with the true process as shown in Fig. 4.1 and Fig. 4.2. 46 4.6. Joint Input-Output Method Case Complexity criterion F Rank of F A FA = (1 Fy1 )T 1 B FB = ( 1 0 Fy1 F r 1 ) 2 C FC = ( 1 1 Fy1 F y 2 ) 2 Table 4.2: The complexity criterion F for the cases listed in Table 4.1 of Example 4.5.1 Case ˆG ˆH A ¯GA = −2.369z −1+2.014z−2 1−0.1165z−1−0.802z−2 ¯HA = 1−0.6437z −1 −0.4329z−1−0.5408z−2 B ¯GB = 0.6023z −1+0.001607z−2 1−1.601z−1+0.6986z−2 ¯HB = 1+0.7978z−11−1.196z−1+0.1965z−2 C ¯GC = 0.6001z −1−0.0003913z−2 1−1.602z−1+0.7014z−2 ¯HC = 1+0.8014z −1 1−1.198z−1+0.1978z−2 Table 4.3: Model estimates for the cases listed in Table 4.1 of Example 4.5.1 4.6 Joint Input-Output Method In this section, a joint input-output closed-loop identification method is described for the closed-loop data set ZNd (4.14). The identification task is to estimate the models ¯G and ¯H (4.9) for the process (4.1). The basic idea of joint input-output closed-loop identification [83] is to consider the control input signal u and the process output y together as the joint outputs (u,y) of the joint inputs (r,e), which consists of the process reference signal r and underlying process disturbance e. The two input signals (r,e) are typically independent. Therefore, the joint outputs (u,y) can be separated as two independent parts: one due to r and one due to e. uk = u r k +u e k, (4.23a) yk = yrk + y e k = G(z)urk +H∗ ◦ ek, H∗ ◦ ek = {G(z)uek +H(z)ek}, (4.23b) where ur and ue are input signal components due to reference signal r and process disturbance e, respectively. yr and ye are corresponding components of process output y. (yr,ur)⊥(ye,ue), and H∗ ◦ e represents the noise portion of y. Then if 47 4.6. Joint Input-Output Method Bode Diagram Frequency (rad/sec) M ag n it u d e , d B P h as e , d eg G GA GB GC -20 -10 0 10 20 30 -180 -90 0 90 180 10-3 10-2 10-1 10 0 101 Figure 4.1: The Bode plots of the true process G, and its estimated process models ¯GA, ¯GB, and ¯GC for the case studies A, B, C, respectively in Example 4.5.1 one can separate (yr,ur) from (y,u) based on measured (u,y,r), a model ¯G can be fit between ur and y. In this work, the two-stage method is considered, i.e. first separate ur from u, then fit a model ¯G between ur and y. 4.6.1 Two-stage Method Here, the two-stage method [91] is briefly described to identify the model ¯G from closed-loop data. If the feedback is a linear causal relation, for instance uk = −κ(z−1)yk, then H∗ ◦ ek in (4.23b) simply becomes H∗ek, and H∗ is the sensitivity not the disturbance model H . In this case yk = G(z)Sur(z)rk +Sye(z)ek, (4.24a) uk = Sur(z)rk +Sue(z)ek, (4.24b) 48 4.6. Joint Input-Output Method M ag n it u d e, d B P h as e, d eg Bode Diagram Frequency (rad/sec) -50 0 50 100 150 -180 -90 0 90 H H A H B HC 10-5 10-4 10-3 10 -2 10-1 100 101 Figure 4.2: The Bode plots of the true noise process H , and its estimated process models ¯HA, ¯HB, and ¯HC for the case studies A, B, C, respectively in Example 4.5.1 and the optimal separation of ur from u is to know the exact Sur , i.e. ur = Surr. Then two-stage method can be summarized as in the Algorithm 3. In the first stage, causal high-order models, for example Finite Impulse Re- sponse (FIR) or Auto-Regressive Exogenous (ARX), are used for ¯Sur[84, 87, 91] while in the second stage, PEM can be used to identify the model ¯G. 4.6.2 First Stage: Projection Method A more general method of joint input-output identification is the projection method [87], in which the controller is possibly nonlinear while the process (G,H) is still linear. In this case, the separation of ur from u is done via fitting a non-causal model between u and r. If infinite impulse model were allowed, then the optimal 49 4.6. Joint Input-Output Method Algorithm 3 The basic steps of two-stage method for closed-loop identification. 1: First stage: Estimate ¯Sur for Sur in (4.24b), based on the data (u,r). Then, calculate ur via ūr = ¯Surr. 2: Second stage: Using the data (ūr,y), identify a model ¯G of G in (4.24a) with the following structure. yk = ¯G(z,θ)ūrk + ¯H∗(z)ek (4.25) solution of min sk ‖ uk− ∞ ∑ k=−∞ skrk ‖ 2, (4.26) s∗k would give ūrk = ∑∞k=−∞ s∗krk, and the residue ūek = uk− ūrk is orthogonal to r. For practical reasons, a high order non-causal FIR is used, then ūrk = ∑M2k=−M1 s∗krk [92]. Once ūr is generated, PEM is often used to fit the model between ûr and y. Then the prediction error for (4.23b) is εk(θ) = ¯H−1∗ (z) ( yk− ¯G(z,θ)ūrk ) . (4.27) The PEM minimizes the cost function VNd (θ) (4.28) with respect to θ . VNd (θ) = 1 Nd Nd∑ k=1 ε2k (4.28) Then w.p.1 as Nd → ∞ and using Parseval’s theorem, the estimation problem be- comes [87] argmin θ∈D ∫ pi −pi | G− ¯G− ¯B |2 Φū r | ¯H∗ |2 dω , (4.29) where Φūr is the spectrum of ūr, and ¯B= GΦūeūrΦūr is the bias term. Thus the consistent estimation of G can be archived when Φūeūr = 0, i.e. optimum projection of u onto r. Hence this estimation problem fits into the class of open-loop identification for the model described by (4.23b) with fixed ¯H∗. The above only presents an overview of the method. Detailed discussion about 50 4.6. Joint Input-Output Method joint input-output identification and projection method can be found in [83, 87, 92]. The important features of the projection method are: 1) consistent estimate of G can be achieved even though controller is non-linear, provided that ur and ue are uncorrelated, and 2) it allows fitting the model to the data with arbitrary frequency weighting in PEM. 4.6.3 Example In this section, the joint input-output identification method is applied to closed- loop data generated by a constrained MPC, where the first stage uses the projection method to obtain excitation and noise contents separation followed by a standard PEM in the second stage. With MPC in the loop, H∗ ◦ ek is not a linearly filtered process. As the unbiasedness of ¯G does not rely on the knowledge of the noise portion but the de-correlation between ur and ue, a linear approximation ¯H∗ek for the noise portion can be used. In that case, ¯H∗ serves as a data filter that can be chosen arbitrarily. Consider the following process with input constraint, yk = G(z)uk +H(z)ek, (4.30a) G(z) = 0.6z −1 +0.45z−2 1−2.5z−1 +2.14z−2−0.63z−3 , H(z) = 1+0.5z−1 1− z−1 , uk ≤ 1, ∀k ≥ 0, (4.30b) and the process noise is ek ∼ N(0,0.1). The constrained MPC is formulated for (4.30) with the following initial model. ˆGinit(z) = −0.1011z−1 +0.7689z−2 1−1.702z−1 +0.8593z−2 , ˆHinit(z) = 1+0.5z−1 1− z−1 . (4.31) Then for the prediction horizon N = 5 and the controller weighting λ = 0.03, the MPC problem (4.3) can written as min û J(yk, û) = N−1 ∑ i=0 (ŷk+i+1− rk+i+1)2 +λ û2k+i (4.32a) s. t : ûk+i ≤ 1, ∀i = 0, . . . ,N−1. (4.32b) 51 4.6. Joint Input-Output Method The objective of the MPC is to improve the loop performance, and the achieved performance is measured by the cost function defined as Vp(y,u,N) = 1 N N ∑ k=1 [y2k +u 2 k]. (4.33) For the above set-up, the cost (4.33) was 1.6978 from 16384 samples simulation, and the major active controller modes were 55.10% from unconstrained case and 38.28% from u = 1. The objective of this simulation is to identify an improved, but still reduced-order (second order) model ¯G. The closed-loop experiment data was generated from a 8192 sample experiment with extra excitation rk ∼ N(0,0.36). The non-causal model used in the first stage is given by ūr = 11 ∑ j=−29 s jrk+ j. Then in the second stage, the process model ¯G1 ¯G1(z) = 0.4472z−1 +0.3871z−2 1−1.821z−1 +0.8859z−2 (4.34) was determined via MATLAB routine b j() with the freedom to choose H∗ among 4th order transfer functions. Next, the closed-loop was simulated with the updated models ( ¯G1, ˆH) in the MPC with the same prediction horizon N and weighting λ . Then the achieved cost was reduced to 0.2650 where the unconstrained controller mode accounted for 96.73% samples while the saturated mode accounted for only 2.48%. 4.6.4 Second Stage: Long Range Predictive Identification In the identification for control context [79], an identification criterion for model estimate bears control performance goals, especially when modeling error is un- avoidable such as the under modeling case, and often can be transformed into stan- dard PEM with special data filters/noise models. Since MPC techniques determine the control action based on predicted future process outputs over a horizon and the performance of an MPC technique largely depends on the model’s prediction 52 4.6. Joint Input-Output Method ability, one may question the use of models which only minimize one-step ahead prediction errors (4.10) in MPC techniques. In [93], it is shown that the combina- tion of the pure Recursive Least Square (RLS) and GPC does not produce a good results. Alternatively, the LRPI method which is used to minimize the averaged prediction error throughout the MPC horizon is introduced in [93]. Good discus- sion about MPC relevant identification methods and LRPI can be found in [94, 95]. In this section, the LRPI is adopted as a MPC relevant criterion in the second stage of identification in the joint input-output identification method. The LRPI criterion for the second stage identification is given by JLRPI = limNd→∞ 1 Nd Nd−N∑ k=1 1 N N ∑ i=1 (yk+i− ȳk+i|k)2, (4.35) where ȳi+1|k is the predicted output from data {y0,...,k, ūr0,...,k+N} with the model (4.40). It is shown in [93] that the criterion (4.35) is equivalent to PEM with a data filter that depends on the noise model parameters. Though the original formula was derived for an ARX model, it can be extended for general models ( ¯G, ¯H∗) as follows. For simplicity, suppose that ¯H∗ = ¯H0 = ¯C(z) ¯D(z) (4.36) is known. Then the output prediction, similar to (4.4), can be written as ȳk+i|k = ¯G(z)urk+i + z−i ¯Fdi (z) ¯C(z) (yk+i− ¯G(z)urk+i), (4.37a) ¯C(z) ¯D(z) = ¯Edi (z)+ z −i ¯Fdi (z) ¯D(z) , (4.37b) where the degrees of the polynomials are deg(Edi )≤ (i−1) and deg(Fdi )= deg(C)− 53 4.6. Joint Input-Output Method 1. Thus the prediction error is yk+i− ȳk+i|k = ( 1− z−i ¯Fdi (z) ¯C(z) )( yk+i− ¯G(z)urk+i ) , = ¯Edi ¯H −1 0 ( yk+i− ¯G(z)urk+i ) . (4.38) Therefore, by Parseval’s theorem, the LRPI criterion (4.35) becomes JLRPI = limNd→∞ 1 Nd Nd∑ k=1 (ε fk (θ)) 2, (4.39) ε fk (θ) = L(z) [ ¯H−10 (yk− ¯Gu r k) ] , L(z)L∗(z−1) = 1 N N ∑ i=1 | ¯Edi (z) | 2 . Since the noise model ¯H0 is assumed to be known, this is a simple filtered PEM problem. Otherwise, in the iterative approaches of seeking better models and con- trol, the current noise model ¯H can be used to simplify matters. The two-stage method with the projection method and the LRPI criterion can be summarized as in Algorithm 4. Algorithm 4 The basic steps of the two-stage method with the projection method and the LRPI criterion for closed-loop identification. 1: First stage: Estimate ¯Sur for Sur in (4.24b) based on the data (u,r), using the following non-causal FIR model structure. ¯Sur : ūrk = M2∑ k=−M1 s∗krk Then, calculate ur via ūr = ¯Surr. 2: Second stage: Using the PEM with LRPI (4.39) as the identification criterion, identify a model ¯G of G in (4.24a) based on the data (ūr,y), with the following structure. yk = ¯G(z,θ)ūrk + ¯H∗(z)ek (4.40) 54 4.7. Conclusion 4.6.5 Example Here, Example 4.6.3 is continued to illustrate the Algorithm 4. In the second stage, the following first order LRPI filter is used. L = 1.605+0.3879z−1 1−0.6083z−1 (4.41) Then a model ¯GLRPI is identified using the PEM with the LRPI filer. ¯GLRPI = 0.4323z−1 +0.3671z−2 1−1.827z−1 +0.9658z−2 . (4.42) For the same initial conditions and noise process used in Example 4.6.3, the achieved cost (4.33) of the new MPC based on ¯GLRPI is 0.2504, slightly better than ¯G1, where the unconstrained controller mode accounted for 96.82% samples and the saturated mode accounted for only 1.93%. The Bode plots of the true, the initial model, and estimated process models are shown in Fig. 4.3. 4.7 Conclusion In this chapter, two closed-loop identification schemes are discussed to apply on the input output data set generated with a constrained MPC in the loop. In the first scheme, the direct closed-loop identification method is investigated, and the infor- mativity criterion, which ensure that the constrained MPC law is complex enough, is presented on the normal operational data set. One would first check this crite- rion, and if it passes, then estimate models using PEM directly with the closed-loop data. Otherwise, on-line experiments have to be conducted to achieve informative data in order to have consistent model estimates. Although the informativity crite- rion is originally presented to the SISO input output process, it can be extended to Multi-Input Multi-Output (MIMO) state space processes. In the second scheme, the joint input-output method is discussed. The projec- tion method is applied in the first stage to fit a non-causal FIR model. A successful first stage then gives the freedom of choosing frequency weighting arbitrarily in the identification criterion during the second stage. In particular, in the case of un- 55 4.7. Conclusion M ag n it u d e, d B P h as e, d eg Bode Diagram Frequency (rad/sec) G0 Ginit G1 GLPRI -40 -20 0 20 40 60 -360 -180 0 180 360 10-3 10-2 10-1 100 101 Figure 4.3: The Bode plots of the true process G0, the initial model Ginit , the estimated transfer functions G1 and GLRPI in Example 4.6.3 and 4.6.5, respectively 56 4.7. Conclusion der modeling, meaningful MPC-relevant identification criterion, such as the LRPI, become viable so that the resultant MPC may give satisfying performance with low controller complexity. Moreover, the projection method might require a large num- ber of parameters in the non-causal FIR model to achieve a better separation of the input signal. One could utilize the explicit piece-wise MPC law in the constrained MPC to revise the first stage that requires less parameters. 57 Chapter 5 Control and Optimization of Two-stage Thermo Mechanical Pulping (TMP) Processes 5.1 Summary In this chapter, NMPC techniques are presented to control and optimize for TMP processes. The TMP process considered in this chapter has two-stages of refining with a primary and a secondary refiners. In the first part of this work, a setpoint tracking problem is formulated with the NMPC technique for the TMP process. Then the performances of ideal-NMPC described in Section 2.3 and asNMPC de- scribed in Section 2.6.2 are compared in a simulation study. In the second part, an econNMPC is formulated for the TMP process. In econNMPC, a dynamic opti- mization problem is formulated to simultaneously regulate and optimize the TMP process in the presence of disturbances. In a simulation study, potential economical benefits of the proposed econNMPC method are demonstrated, through a reduction in total specific energy. The computational burden of the resulting NLP problems in NMPC formulations are handled with the IPOPT solver described in Section 2.6.1. 5.2 Introduction MPC has been applied to control and optimize chemical processes [68], [65]. One of the main advantages of MPC is its ability to cope with hard constraints on con- trol and state variables. Therefore, MPC has attracted considerable research effort, 58 5.2. Introduction and has been widely applied in process industries for low-level regulatory control under constraints [96], [45]. The main issues with the applications of NMPC with detailed nonlinear models were due to computational burden and lack of theoretical support for stability [63, 64]. However, with the improvement of today’s computa- tional power and development of efficient NLP solvers [59] and tailored algorithms [58], there is growing tendency to use NMPC in the optimization layer of the deci- sion making hierarchy, with detailed nonlinear dynamic models [65, 67, 68]. Computer control of wood-chip refiners has been studied in literature and in- dustry for a long time [10], and there have been a few attempts to control and optimize TMP processes using MPC strategy. Ruscio [97] [7] demonstrated a lin- ear MPC control technique for a TMP process. The objective of that work was to control specific energy and consistency of the TMP process using screw speed, di- lution flow and plate gap. First, a subspace identification technique was presented, and then a linear state-space model was identified with time series data to use for MPC. Du et al. [98] presented a GPC algorithm based on a nonlinear Laguerre model to control motor load and consistency of a TMP process. In addition, there are few applications reported to control pulp quality and reject refining process using MPC [3], [99]. Recently, Sidhu et al. [8] applied MPC strategy to an indus- trial two-stage TMP process. They presented a two-level control strategy in order to control motor loads and consistencies, and to improve pulp quality. The strat- egy was to first decouple the control variables according to their bandwidths, and then to formulate linear stabilizing and quality controllers. However, the separation was merely due to the placement of the pulp quality sensor after the latency tank, but not due to the dynamics of the refiners. The coordination between the quality controller and the stabilizing controller was performed in steady-state. In contrast to the above developments, a general NMPC framework to control and optimize a two-stage TMP process is considered in this chapter. First, NMPC for setpoint tracking problem is designed, and then an economical objective is in- tegrated in the NMPC formulation. The economical benefit is realized through a reduction in specific energy consumption of the TMP process. The resulting NLP problems are solved using the IPOPT solver described in Section 2.6.1. in order to handle the computational burden. Further, the computational issue is improved with asNMPC controller concept described in Section 2.6.2. 59 5.2. Introduction 5.2.1 TMP Process Modeling TMP processes are, inherently, MIMO processes with complex dynamics and se- vere interactions among process variables. Modeling of TMP is challenging due to complex conditions inside the pulp refiners. It is not yet clearly understood what the desired mechanism for efficient fiber development is and how energy trans- forms into the pulp inside the refiners [100]. In this chapter, the following models are used for the process variables to develop a discrete-time nonlinear model for the TMP process. Production Rate Production rate or throughput of the TMP process is proportional to chip-transfer screw speed and is simply calculated by multiplying the chip-transfer screw speed by a constant in most industrial processes [12]. Changes in chip bulk density and chip solid content are considered as more dominant disturbances in TMP processes. In order to consider these variations, here the following model is used for the pro- duction rate, P (tonnes/day) [6]. P = 1.44kpscdcR, (5.1) where kp (m3/rev) is proportional constant, sc (%) is chip solid content, dc (kg/m3) is chip bulk density, and R (rpm) is chip-transfer screw speed, respectively. Motor Loads The plate gap, which determines refining zone surface area to volume ratio, is a key variable for modeling the motor load [6], [11]. The number of wood fibers per unit volume entering the refiner, which is a function of production rate and dilution water flow rate, is another variable which determines the motor load. In this work, the following model1 is used to describe motor load as a function of the production rate, plate gap, and dilution flow rate. 1This model has been used in the case study of MPC Supervisory Control of a Two Stage Thermo- Mechanical Pulping Process, Model Predictive Control Toolbox DEMOS, MATLAB 7.8.0 (R2009a). 60 5.2. Introduction Mi = kmiP Di (1− exp(−10Gi))(ai−biGi), (5.2) where Mi (MW ) is the motor load, Di (l/min) is the dilution water flow rate, Gi (mm) is the gap, and kmi, ai and bi are parameters of the refiner. The subscript index i = p,s represent the primary and secondary refiners, respectively. Consistencies Refining consistency is defined as the ratio of dry fiber flow rate to fiber and water flow rate. The refining consistency of each refining process changes with produc- tion rate, motor load, and dilution flow rate and can be derived using mass balance [7], [6]. In this work, it is assumed that water and energy losses in the refiners are negligible, and back-flow of steam is a function of motor load for each refiner. Then outlet consistencies of the primary and secondary refiners can be expressed as Cp = 100P P+1.44Dp− kepMp , (5.3a) Cs = 100P (P/(0.01Cp)+1.44Ds− kesMs) , (5.3b) where kep and kes are parameters for the primary and secondary refiners, respec- tively. Pulp Quality Variables In this work, pulp quality variables such as Canadian Standard Freeness (CSF), Long Fiber Content (LFC), and Shive Content (SC) are considered after the sec- ondary refining process, as the target pulp quality variables. For the pulp quality variables, the nonlinear model developed in [101] in terms of specific energy and specific refining power is used. The specific energy and specific refining power are implicit functions of the production rate, motor loads and consistencies of the both primary and secondary refiners. Therefore, for simplicity, the relationship between 61 5.2. Introduction pulp quality variables and the production rate, motor loads and consistencies is expressed as a general nonlinear implicit function gtmp as follows. 0 = gtmp(CSF,LFC,SC,P,MP,Cp,Ms,Cs) (5.4) Here, the aim is to form a nonlinear discrete-time model for the TMP refining process; both primary and secondary refiners are considered together as a single process model. The production rate, primary motor load, primary consistency, secondary motor load, and secondary consistency are considered as the discretized differential state variables while the pulp quality variables: CSF, LFC, and SC, are considered as the discretized algebraic state variables of the TMP refining process model. Manipulated variables are chip-transfer screw speed, plate gap and dilution water flow rates of the both primary and secondary refiners. By introducing linear dynamics for the discretized differential state variables and superimposing it on the steady-state relationships (5.1)-(5.3), a Wiener type discrete-time nonlinear model for this TMP refining process can be formulated at sample time k as follows. zk+1 = Atmpzk +htmp(zk,vk) (5.5a) 0 = gtmp(zk,yk) (5.5b) where zk ∈Rnz , yk ∈Rny , and vk ∈Rnv denote the vectors of discretized differential state, algebraic state, and manipulated variables of the TMP refining model. nz, ny, and nv are dimensions of the variables z, y, and v, respectively. Then 62 5.2. Introduction z =   Production rate Primary refiner motor load Primary refiner consistency Secondary refiner motor load Secondary refiner consistency   , y =   CSF after primary refining LFC after primary refining SC after primary refining CSF after secondary refining LFC after secondary refining SC after secondary refining   , and v =   Chip-transfer screw speed Primary refiner plate gap Primary refiner dilution water flow rate Secondary refiner plate gap Secondary refiner dilution water flow rate   . (5.6) Atmp ∈ Rnz×nz is the dynamic matrix which can be identified for the TMP pro- cess by using any linear system identification method. On the other hand, for simplicity, one can use the time constant and time delay information of each sub- processes to form Atmp as follows. ¯Atmp(z) =   gps(z) 0 0 0 0 0 gpmpg(z) 0 0 0 0 0 gpcpd(z) 0 0 0 0 0 gsmsg(z) 0 0 0 0 0 gscsd(z)   (5.7) where ¯Atmp(z) is the dynamic transfer function matrix of the TMP process. gps(z) is the transfer function between production rate and chip-transfer screw speed. gpmpg(z) is the transfer function between primary refiner motor load and primary refiner gap. gpcpc(z) is the transfer function between primary refiner consistency and primary refiner dilution water flow rate. gsmsg(z) is the transfer function be- tween secondary refiner motor load and secondary refiner plate gap. gscsd(z) is the transfer function between secondary refiner consistency and secondary refiner di- lution water flow rate. The transfer functions of the subprocesses have following 63 5.2. Introduction form. gps = b1z−d1 z−a1 (5.8) gpmpg = b2z−d2 z−a2 (5.9) gpcpd = b3z−d3 z−a3 (5.10) gsmsg = b4z−d4 z−a4 (5.11) gscsd = b5z−d5 z−a5 (5.12) where ai and di for i = 1 . . .5 are the poles and time delays of the subprocess, respectively. bi for i = 1 . . .5 are parameters, and can be found as bi = 1− ai, for unity dynamic gains of the each subprocesses. Then the dynamic matrix can be expressed as Atmp = diag([a1,a2,a3,a4,a5]T ). (5.13) The nonlinear state function, htmp :Rnz+nv 7→Rnz , and htmp = (I−Atmp)Htmp where Htmp maps steady-state control input variables to steady-state differential state vari- ables of the TMP refining model. Note that the time delay information is repre- sented in the the state function htmp as htmp(zk,vk−d) where d = [d1,d2,d3,d4,d5]T is the time delay vector in sample times. In this work, we assigned time constants and time delays of subprocesses as given in Table 5.1, and then ai and bi are calcu- lated. Subprocess Time constant (s) Time delay, d a b gps 2 2 0.6065 0.3935 gpmpg 2 1 0.6065 0.3935 gpcpd 2 0 0.6065 0.3935 gsmsg 2 1 0.6065 0.3935 gscsd 2 0 0.6065 0.3935 Table 5.1: Time constants and time delays of subprocesses of the TMP process 64 5.3. NMPC for Setpoint Tracking For simplicity, a general form of state equation for (5.5a) is considered in fol- lowing NMPC formulations. zk+1 = ftmp(zk,vk) (5.14) Note that in the remainder of this chapter, the TMP process is simulated up to the secondary refiner, assuming the pulp quality measurements are available at that point. The following section describes the NMPC strategy for setpoint tracking problem for the TMP process. 5.3 NMPC for Setpoint Tracking Using the current process state xk, the dynamic optimization problem is formulated for the NMPC at sample time k as follows. min v Jspt , Jspt := N ∑ i=1 ‖ξi−ξ sk+i‖2Qξ +‖vi−usk+i−1‖2Qv (5.15a) s. t. : z0 = xk, (5.15b) zi = ftmp(zi−1,vi−1), (5.15c) 0 = gtmp(zi−1,yi−1), (5.15d) zmin ≤ zi ≤ zmax, (5.15e) ymin ≤ yi ≤ ymax, (5.15f) umin ≤ vi−1 ≤ umax, i = 1, . . . ,N (5.15g) where N is the length of the prediction horizon in sample times. v= {v0,v1, . . . ,vN−1} a feasible control moves which satisfy the constraints (5.15b)-(5.15g). zmin and zmax, ymin and ymax, and umin and umax are lower and upper bound vectors of differ- ential and algebraic state, and manipulated variables of the TMP process, respec- tively. ξ ∈Rnξ ⊂Rnz+ny is a vector of controlled variables. With the TMP process, the discretized differential variables - production rate, motor loads and consisten- cies of both primary and secondary refiners - are the controlled variables, i.e. ξ ≡ z, while the discretized algebraic state variables: CSF, LFC, SC are controlled within a desired range given in (5.15f). The diagonal matrices Qξ and Qv are weights for 65 5.3. NMPC for Setpoint Tracking controlled state variables and input variables, respectively. ξ s is the optimum set- point vector, obtained normally solving a steady-state optimization, for controlled variables, and us is the corresponding reference input vector for the control vari- ables. Although the use of reference input vector us is not necessary if there is no degree of freedom in the controller i.e. nξ = nv, this information might be useful when changing the operating conditions especially for nonlinear processes [27]. 5.3.1 Example 1 In this study, the closed-loop is simulated for the TMP process for the setpoint tracking with the NMPC technique. The simulation was set-up so that there was no process/model mismatch and disturbances. In the beginning of the simulation, the production rate, primary motor load and consistency, and secondary motor load and consistency are set to 208 tonnes/day, 10.5 MW and 50 %, and 4.75 MW and 43%, respectively. At time 50s, new target setpoints are introduced as 250 tonnes/day, 8.8 MW, 44 %, 4.82 MW and 37.2%, respectively, along with the new reference input vector, for the setpoint tracking NMPC formulation. The sampling time is 1 s, and the control and prediction horizons are set to 10 s. The controlled and manipulated weighting matrices are kept as Qξ = diag([0.01 1 1 1 1]T ) and Qv = diag([0.8 8 0.8 8 0.8]T ×10−3). First, the resulting NLP problem was solved using fmincon solver in MATLAB with interior-point option. In this case, it took 5−15 s CPU seconds to find the optimum solution. In contrast, the IPOPT solver took 0.02−0.04 s to find the optimum solution. The simulation studies were performed with Intel Core 2 CPU with 2.4 GHz PC. The computational delay for NMPC optimization introduces feedback delay in the system. This feedback delay may cause poor performance and deterioration of the stability of the system [102], [103]. The computational burden is one of main drawbacks of NMPC techniques when applied to real systems. In [103, 104], it showed how this computational delay can be explicitly included in the NMPC algo- rithm to ensure minimal degradation of the performance, and moreover, discussed the conditions under which the stability of the closed-loop can be guaranteed. In this chapter, it is assumed that states can be directly measured for the TMP process. However in practice, all states of the TMP process, for example outlet consisten- 66 5.4. NMPC for Economical Operation cies of the refiners, are not measured, but need to be estimated. Incorporation of a state estimation technique to NMPC formulation would increase the computational burden of the overall problem. In the following example, the closed-loop simulation was performed to demon- strate the performance of the asNMPC technique described in Section 2.6.2 as a way to reduce on-line computational time, along with the ideal-NMPC technique described in Section 2.3 for the setpoint tracking problem of the TMP process. 5.3.2 Example 2 This simulation study considers the presence of random disturbances on the state variables of the TMP process to illustrate the performance of asNMPC technique. The asNMPC technique computes approximate solution to the resulting NLP prob- lem of NMPC quickly. The ideal NMPC and asNMPC have identical performance under nominal conditions. The difference between the solution of ideal-NMPC and its approximation through asNMPC is identified by the second order effects of model/process mismatch and disturbances [58], which in turn leads to predic- tion error of the states of the process. As in the first simulation study 5.3.1, the asNMPC and ideal-NMPC techniques are given new setpoint targets, along with the corresponding reference manipulated inputs vector, at time 50s. At 100s, the targets are set back to their original values. The tuning parameters, N, Qξ , and Qv, are the same as in Example 5.3.1. The IPOPT solver took 0.02 s to solve the NLP problem in ideal-NMPC while the on-line computational time of asNMPC was 0.004s. By providing the values of the variables of the previous solution, the warm-start strategy was used in the both cases. The simulation results are shown in Fig. 5.1 and Fig. 5.2. The simulation results illustrate that the asNMPC tech- nique was able to track the setpoints, similar to the ideal-NMPC technique, in the presence of disturbances. 5.4 NMPC for Economical Operation As discussed in Section 3.2, the conventional two-layer architecture for advanced control systems has several drawbacks, mainly due to 1) low frequency of opera- 67 5.4. NMPC for Economical Operation 0 50 100 150 150 200 250 300 P ro d u ct io n r at e to n n e/ d ay 0 50 100 150 8 10 12 P ri m ar y m o to r lo ad M W 0 50 100 150 40 45 50 55 P ri m ar y c o n si st en cy % 0 50 100 150 4 4.5 5 5.5 S ec o n d ar y m o to r lo ad M W 0 50 100 150 30 35 40 45 S ec o n d ar y c o n si st en cy % Time, s Time, s Time, s Time, s Time, s asNMPC ideal-NMPC set-point Figure 5.1: Controlled variables of the TMP process: Performance comparison of asNMPC and ideal-NMPC techniques for setpoint tracking problem. 68 5.4. NMPC for Economical Operation 0 50 100 150 10 15 20 25 30 S cr ew s p ee d R P M 0 50 100 150 0 0.5 1 1.5 2 P ri m ar y g ap m m 0 50 100 150 100 150 200 250 300 P ri m ar y d il u ti o n f lo w l /m in 0 50 100 150 0 1 2 3 S ec o n d ar y g ap m m 0 50 100 150 0 50 100 150 S ec o n d ar y d il u ti o n f lo w l /m in Time, s Time, s Time, s Time, s asNMPC ideal-NMPC set-point Figure 5.2: Manipulated inputs of the TMP process: Performance comparison of asNMPC and ideal-NMPC techniques for setpoint tracking problem. 69 5.4. NMPC for Economical Operation tion of the RTO layer, and 2) assumption of slow variation of the disturbances on the process. In particular, disturbances in wood chip species, chip size distribution, the chip solid content, and the chip bulk density are rapidly changing in the TMP process. The consequences of the slow operation of the RTO layer and the distur- bances are that the operating setpoints will be sub-optimal and also the economical optimum may not even be at steady-state. In this section, the economical objective is integrated with the setpoint tracking objective in the formulation of NMPC technique, as discussed in 3.3 and compare the economical benefits with the setpoint tracking NMPC technique when distur- bances are present in the TMP process. The simulation studies are performed with the asNMPC control strategy along with IPOPT solver. The specific energy, de- fined as in Section (5.16), is considered as the economical indicator for the TMP process. The overall economical objective is to minimize the total specific energy while keeping the pulp quality variables within a specific range. Specific Energy = Total motor load Production rate MW/(tonnes/day) = Primary motor load + Secondary motor load Production rate . (5.16) Thus the objective function for economically-oriented NMPC formulation (econ- NMPC) is described by Jecon = αtmp N ∑ i=1 Specific Energy+βtmp N ∑ i=0 ‖ξi−ξ sk+i‖2Qξ +‖vi−1−usk+i−1‖2Qv . (5.17) The αtmp and βtmp are positive weights for the specific energy and the tracking objectives, respectively. Throughout the rest of this chapter, sptNMPC and econ- NMPC denote for the NMPC formulations (5.15) with setpoint tracking objective function (5.15a), and with the economically-oriented objective function (5.17), re- spectively. With the objective function (5.17), the econNMPC technique tries to minimize the specific energy consumption and simultaneously tries to minimize 70 5.4. NMPC for Economical Operation transition time to reach the target setpoint from the current state. Note that first, it was to try to use only the economic objective to formulate the NMPC, i.e. βtmp = 0 in (5.17) because the overall objective is to minimize the total specific energy con- sumption while respecting the process constraints. However, as noted by several researchers [65] [105], the difficulty in solving the resulting NLP problem was experienced. This is simply because with βtmp = 0 it cannot be guaranteed the uniqueness of the solution of the resulting NLP problem, the second order opti- mality conditions might not hold, and it may lead to ill-posed search directions in the IPOPT solver which makes the NMPC problems difficult to solve [65]. By adding the setpoint tracking objective to the econNMPC formulation in (5.17), the objective function is regularized. In this case, it is necessary to set αtmp and βtmp as 86 and 1, respectively. The following closed-loop simulation studies demon- strate the potential economical benefits of the econNMPC formulation against the sptNMPC formulation when the disturbances are present in the TMP process. 5.4.1 Example 1 In this study, changes in the chip bulk density are considered as a disturbance to the TMP process; variation in the chip bulk density is considered as the main cause of disturbances in this process [6]. Changes in the chip bulk density dc in (5.1) in turn cause some variations in the production rate (5.1), motor loads (5.2) and consistencies (5.3) instantaneously and hence in specific energy (5.16). This is mainly because the TMP refiners operate based on mass throughput [8]. First, the chip bulk density is set to 90% of its nominal value until 60s of the simulation. Then, the chip bulk density is increased up to 110% of it nominal value for the period of 1 min. Then it is set back to 90% of the nominal value at 120s. Random disturbances on state variables were not considered in order to demonstrate the main effects of the chip bulk density changes. The sampling time is 1s, and the prediction and control horizons are set to 10s for the both sptNMPC and econNMPC techniques. Since the overall economical objective is to minimize the specific energy of the TMP process, the controllers are tuned so that there is some freedom to optimize the specific energy, rather than being constrained by tightly controlled the production rate, motor loads [8]. 71 5.4. NMPC for Economical Operation 0 20 40 60 80 100 120 140 160 180 220 240 260 280 0 20 40 60 80 100 120 140 160 180 6 7 8 9 0 20 40 60 80 100 120 140 160 180 4 4.5 5 P ro d u ct io n r at e to n n e/ d ay P ri m ar y m o to r lo ad M W S ec o n d ar y m o to r la o d M W econNMPC sptNMPC set-point Time, s Time, s Time, s Figure 5.3: The production rate, and motor load responses of the TMP process to the change of the chip bulk density 72 5.4. NMPC for Economical Operation 0 20 40 60 80 100 120 140 160 180 2.2 2.4 2.6 2.8 3 0 20 40 60 80 100 120 140 160 180 1.4 1.5 1.6 1.7 1.8 1.9 0 20 40 60 80 100 120 140 160 180 3.5 4 4.5 5 Time, s Time, s Time, s econNMPC sptNMPC P ri m ar y s p ec if ic e n er g y G J/ to n n e S ec o n d ar y s p ec if ic e n er g y G J/ to n n e T o ta l sp ec if ic e n er g y G J/ to n n e Figure 5.4: The specific energies of the both primary and secondary TMP pro- cess: performance comparison of the sptNMPC with econNMPC techniques to the change of the chip bulk density 73 5.4. NMPC for Economical Operation 0 20 40 60 80 100 120 140 160 180 100 110 120 130 140 150 62 63 64 65 66 61 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 0.75 0.8 0.85 0.9 0.95 S ec o n d ar y sh iv e co n te n t, % S ec o n d ar y l o n g fi b er co n te n t, % S ec o n d ar y fr ee n es s, C S F Time, s Time, s Time, s econNMPC sptNMPC Figure 5.5: The freeness, long fiber content and shive content responses of the TMP process to the change of the chip bulk density. 74 5.4. NMPC for Economical Operation The closed-loop responses are shown in Fig. 5.3 - 5.5. Fig. 5.3 shows the production rate, primary motor load and secondary motor load responses while Fig. 5.4 shows the specific energy changes during the simulation. As shown in Fig. 5.5, the both asNMPC and econNMPC techniques were able to keep the pulp quality variables: freeness, long fiber content, and shive content, within the operating limits of (0− 175) CSF, (50− 70%), and (0− 1.2%), respectively. As it is clear from the simulation results, there is a potential reduction, about 13%, in total specific energy in the system with the econNMPC technique, in contrast to the system with the sptNMPC technique. The changes in the chip bulk density and the chip solid content are considered as disturbances in Example 5.4.2, where the closed-loops are simulated with some random disturbances on the state variables. 5.4.2 Example 2 The objective of this simulation study is to show the economic benefits of the econ- NMPC technique with the presence of state disturbances. The simulation set-up is same as Example 5.4.1. Here, changes in the chip solid content is considered as an additional disturbance from 180s to 300s of the simulation. The variation in the chip solid content sc in (5.1) or chip moisture content is also a key disturbance, which causes production rate (5.1), motor loads (5.2), and consistencies (5.3) to change in the TMP process. Up to 180 s the chip bulk density changes are the same as in Example 5.4.1. After that, from 180s to 240s the chip solid content is set to increase to 110% and decrease to 90% of its nominal value at 240s, while keeping the chip bulk density as 90% of its nominal value. The closed-loop sim- ulation results are shown in Fig. 5.6 - 5.8. Fig. 5.6 shows the production rate, motor loads and consistencies for the both primary and secondary refiners. As it is clearly seen from Fig. 5.7, there is still reduction in specific energy, of about 12%, in the TMP process, even with some random disturbances on the state vari- ables. The both asNMPC and econNMPC techniques were able to keep the pulp quality variables: freeness, long fiber content, and shive content, within the oper- ating limits of (0−175) CSF, (50−70%), and (0−1.2%), respectively. Here the pulp quality variables are considered being measured once a second just after the secondary refiner. 75 5.5. Conclusion 5.5 Conclusion In this chapter, an NMPC strategy for the setpoint tracking problem, described in Section 2.3, for the TMP refiner process is first presented. The burden of compu- tational delay for NMPC was handled by using the asNMPC concept, described in Section 2.3, with the IPOPT solver described in Section 2.6.2. The performances were the same as the system with ideal NMPC strategy 2.3. Next, an economi- cal objective was integrated to the NMPC formulation to dynamically optimize the TMP process with the presence of disturbances. A potential economic benefit was realized through reductions, of about 13%, in the total specific energy of the TMP process when chip bulk density changes with time, compared to the closed-loop system with the setpoint tracking NMPC technique. In the final simulation study, it was realized approximately 12% reduction in the total specific energy when the TMP process undergoes the both chip bulk density and chip solid content changes, and some random disturbances on the state variables. 76 5.5. Conclusion 0 50 100 150 200 250 300 220 240 260 280 P ro d u ct io n r at e t o n n e/ d ay 0 50 100 150 200 250 300 6 7 8 9 Time, s P ri m ar y m o to r lo ad , M W 0 50 100 150 200 250 300 35 40 45 50 Time, s P ri m ar y co n si st en cy , % 50 100 150 200 250 300 3 4 5 6 Time, s S ec o n d ar y m o to r lo ad , M W 0 50 100 150 200 250 300 30 35 40 Time, s S ec o n d ar y co n si st en cy , % Time, s econNMPC sptNMPC set-point 0 Figure 5.6: The production rate, consistencies and motor loads responses of the TMP process to the change of the chip bulk density and chip solid content, with random disturbances on the state variables 77 5.5. Conclusion 0 50 100 150 200 250 300 2 2.5 3 3.5 P ri m ar y r ef in er sp ec if ic e n er g y , G J/ to n n e 0 50 100 150 200 250 300 1.4 1.6 1.8 2 S ec o n d ar y r ef in er s p ec if ic e n er g y , G J/ to n n e 0 50 100 150 200 250 300 3 3.5 4 4.5 5 T o ta l sp ec if ic e n er g y G J/ to n n e 1.2 Time, s Time, s Time, s econNMPC sptNMPC Figure 5.7: The specific energies of the both primary and secondary the TMP pro- cess: performance comparison of the sptNMPC and econNMPC techniques to the change of the chip bulk density and chip solid content, with random disturbances on the state variables 78 5.5. Conclusion 0 50 100 150 200 250 300 80 100 120 140 160 Time, s S ec o n d ar y fr ee n es s, C S F 0 50 100 150 200 250 300 58 60 62 64 66 Time, s 0 50 100 150 200 250 300 0.7 0.75 0.8 0.85 0.9 0.95 Time, s S ec o n d ar y lo n g f ib er c o n te n t, % S ec o n d ar y sh iv e co n te n t, % econNMPC sptNMPC Figure 5.8: The freeness, long fiber content and shive content responses of the TMP process to the change of the chip bulk density and chip solid content, with random disturbances on the state variables. 79 Chapter 6 Control and Optimization of TMP Processes with Multi-stage Low Consistency Refining 6.1 Summary In this chapter, nonlinear model based dynamic optimization approaches for opti- mal control for TMP processes are presented. Two types of TMP refining processes are considered. The first TMP process is a conventional refining process with two- stages of HC refining with a primary and a secondary refiners followed by a latency chest and a third-stage LC refiner, while the second TMP process has only a HC primary refiner, and the secondary stage HC refiner is replaced by multiple-stages of LC refiners after the latency test. Both TMP plants are dynamically optimized against disturbances to produce the same target final pulp qualities while minimiz- ing specific energy consumptions of the processes. In this work, a NMPC technique with an economic objective function is considered. The simulations result shows the potential economical benefits of an improved TMP process, the TMP process with multi-stage LC refining, through a reduction in total specific energy, com- pared to a conventional TMP process. The NMPC controller is implemented with the asNMPC concept described in Section 2.6.2, and by using the IPOPT described in Section 2.6.1. 80 6.2. Introduction 6.2 Introduction The advanced control of TMP processes, which aims at optimal operation of the process, has become of prime important [5, 8, 14, 67]. The primary objective of these advanced control systems is to minimize energy consumption while reducing variability in pulp quality against inevitable disturbances. Most of the industrial implementations of advanced TMP control systems has been considered for main line primary and secondary refiners and, in some instances, for main line reject refiners and screening processes [8, 106]. The mainline refiners are, however, high energy-consuming; approximately 60% of total energy in the TMP process is con- sumed by these refiners as they operate at HC [1]. Therefore, in recent years, researchers have leaned toward developing LC refining technology, which would be more energy efficient than conventional HC refining [1, 18]. At the same time, the pulping industry tends to install LC refiners in the main line, most commonly following the latency chest and before the screening process [18, 107, 108]. The inclusion of the third stage LC refiner in the mainline has led to the reconsideration of advanced control system to optimally coordinate the HC and LC refiners. In this chapter, the econNMPC framework is investigated for two TMP pro- cesses to dynamically optimize the processes when the inevitable disturbances present. In Scheme A, a conventional TMP refining process is considered, where it has two-stages of HC refining with a primary and a secondary refiners followed by a latency chest and a third-stage LC refiner. In Scheme B, the improved TMP process [1] only has one HC primary refiner, and the secondary stage HC refiner is replaced by multiple-stages of LC refiners after the latency test. Both the schemes are compared with setpoint tracking NMPC and econNMPC techniques. In this work, the resulting NLP problem is solved using the IPOPT solver, described in Section 2.6.1, in order to handle the computational burden. The feedback compu- tational is further improved with the asNMPC concept, described in Section 2.6.2. To the best of author’s knowledge to date, this work is one of the first contribu- tions that investigates dynamic optimization of TMP processes with multi-stage LC refining using an economically oriented NMPC framework. 81 6.3. Process Description and Modeling 6.3 Process Description and Modeling In this work, the HC refiners are modeled as described in Section 5.2.1. In the latency chest, latency removal is characterized by freeness drop, and LFC and SC are assumed remain unchanged. The model for latency chest freeness presented in [109] is used in this work. For LC refiners, recently developed steady state models in [1] are used. 6.3.1 Latency Chest Modeling Consistency in the latency chest can be modeled using mass balance. Assuming water and energy losses are negligible, the following model can be expressed, sim- ilar to the model described in (5.2.1). Cl = 100P (P/(0.01Cs)+1.44Dl) , (6.1a) where Cl(%) is latency chest consistency, P(tonne/day) production rate, Cs(%) secondary refiner consistency, and Dl(l/min) latency chest dilution water flow rate. The model presented in [109] for CSF in the latency chest is adapted for steady state modeling as follows. CSFoutl =CSF inl  1− 1( 1+al exp −( bl Tl ) )   , (6.1b) where CSFoutl and CSF inl are freeness of the pulp after and before the latency chest, respectively. Tl(◦C) is the temperature of the latency chest. al and bl are parame- ters. 82 6.4. Scheme A: Conventional TMP Process with a Third-stage LC Refiner 6.3.2 LC Refiner Modeling In a most recent work [1], the following models are presented for motor load, LFC, and CSF for LC refiners, respectively. Mt = kmt ( Dt Gt )at (FLin Dt )bt , (6.2a) where Mt(MW ) is the LC refiner motor load, Dt(m) is the diameter of the LC refiner, Gt(mm) is the gap, FLin(mm) is the length of the incoming fiber to the LC refiner. kmt , at , and bt are parameters for a given LC refiner. CSFoutt =CSF int exp ( CSFint D3t )−ct (DtGt ) dt NSEt ω2D2t (6.2b) where CSFoutt (CSF), and CSF int (CSF) are the freeness of the pulp before and after the LC refining, respectively. NSEt(kWh/t) is the net specific energy, and ω(1/s) is the LC refiner rotational speed. ct and dt are parameters for a given LC refiner. LFCoutt = { al f cLFCint if Gt ≥ Gcrit ;( hl f c ln ( Gt Gcrit ) +al f c ) LFCint if Gt < Gcrit , (6.2c) where LFCoutt (%) and LFCin(%) are LFC after and before LC refining, respec- tively. Gcrit(mm) is the critical gap. al f c and hl f c are parameters. 6.4 Scheme A: Conventional TMP Process with a Third-stage LC Refiner Two-stage TMP processes are most common in pulp and paper industry, and have extensively been studied for optimal operation [8, 67]. However, with the intro- duction of more energy efficient LC refining technology, a third stage LC refiner is installed just after the latency chest as shown in Fig. 6.1. The precise physical mechanism for the increased efficiency of LC refining is not yet clearly understood [1]. Despite the lack of a theoretical understanding of LC refining, LC refiners are being applied in the pulp and paper industry, and significant energy savings 83 6.4. Scheme A: Conventional TMP Process with a Third-stage LC Refiner Screening, Steaming Washing, Preheating 1st stage HC refining Wood chips Pulp Wood chips pretreatment Refining 2st stage HC refining 3rd stage LC refining Screening, Reject refining Latency removal Dewatering Pulp processing Figure 6.1: Scheme A: basic schematic of a conventional TMP process with a third stage LC refiner are reported [108]. In this work, a Wiener type discrete-time nonlinear model for this TMP refining process is used in the NMPC formulation, similar to the model described in Section 5.5. zk+1 = AAtmpzk +hAtmp(zk,vk) (6.3a) 0 = gAtmp(zk,yk), (6.3b) where (·)A denotes Scheme A, zk ∈ Rnz , yk ∈ Rny , and vk ∈ Rnv denote the vectors of discretized differential state, algebraic state, and manipulated variables of the TMP refining model. nz, ny, and nv are dimensions of the variables z, y, and v, respectively. Then 84 6.4. Scheme A: Conventional TMP Process with a Third-stage LC Refiner z =   Production rate Primary refiner motor load Primary refiner consistency Secondary refiner motor load Secondary refiner consistency Consistency after latency chest CSF after latency chest Tertiary refiner motor load CSF after tertiary refiner   , y =   CSF after primary refining LFC after primary refining SC after primary refining CSF after secondary refining LFC after secondary refining SC after secondary refining LFC after latency chest SC after latency chest LFC after tertiary refiner SC after tertiary refiner   , and v =   Chip-transfer screw speed Primary refiner plate gap Primary refiner dilution water flow rate Secondary refiner plate gap Secondary refiner dilution water flow rate Latency chest temperature Latency chest dilution water flow rate Tertiary refiner plate gap   . (6.4) AAtmp ∈Rnz×nz is the dynamic matrix, and the nonlinear state function, hAtmp :Rnz+nv 7→ R nz , maps steady-state control input variables to steady-state differential state vari- ables of the TMP refining model. gAtmp is the nonlinear implicit function of variables z and y. For simplicity, a general form of state equation for (6.3a) is considered in following NMPC formulations. zk+1 = f Atmp(zk,vk) (6.5) Note that in Scheme A, the TMP process is simulated up to the tertiary LC refiner, assuming the pulp quality measurements are available at that point. The following section describes the NMPC strategy for setpoint tracking problem for the TMP process 6.1. 85 6.4. Scheme A: Conventional TMP Process with a Third-stage LC Refiner 6.4.1 NMPC for Setpoint Tracking and Economical Operation Similar to Section 5.3, using current process state xk, the dynamic optimization problem is formulated for the NMPC at sample time k as follows. min v JAspt , JAspt := NA ∑ i=1 ‖ξ Ai −ξ Ask+i‖2Qξ A +‖vi−1−uAsk+i−1‖2QvA (6.6a) s. t. : z0 = xk, (6.6b) zi = f Atmp(zi−1,vi−1), (6.6c) 0 = gAtmp(zi−1,yi−1), (6.6d) zmin ≤ zi ≤ zmax, (6.6e) ymin ≤ yi ≤ ymax, (6.6f) umin ≤ vi−1 ≤ umax, i = 1, . . . ,NA (6.6g) where NA is the length of the prediction horizon in sample times. v= {v0,v1, . . . ,vNA−1} a feasible control moves which satisfy the constraints (6.6b)-(6.6g). zmin and zmax, ymin and ymax, and umin and umax are lower and upper bound vectors of differential and algebraic state, and manipulated variables of the TMP process (6.1), respec- tively. ξ A ∈ Rnξ ⊂ Rnz+ny is a vector of controlled variables. For Scheme A, ξ A ≡ z. Since the final freeness, after the tertiary LC refiner, is controlled to a desired value, it is considered as a controlled variable in this formulation. The diagonal matrices Qξ A and QvA are weights for the controlled state variables and input variables, respectively. ξ As is the optimum setpoint vector, obtained solving a steady-state optimization, for controlled variables, and uAs is the corresponding reference input vector for the control variables. In the econNMPC formulation, the following objective function (6.7) is minimized whereas in the NMPC optimization (6.6), JAspt is minimized. JAecon = αAtmp NA ∑ i=1 Specific Energy+β Atmp NA ∑ i=1 ‖ξi−ξi+1‖2Qξ A +‖vi−1− vi‖2QvA , (6.7) 86 6.4. Scheme A: Conventional TMP Process with a Third-stage LC Refiner where Specific Energy is defined as Specific Energy = Total motor load Production rate MW/(tonnes/day) = Primary motor load + Secondary motor load + Tertiary motor load Production rate . (6.8) αAtmp and β Atmp are positive weights for the specific energy and the tracking objec- tives, respectively. Similar to the description in Section 5.3, with this objective function (6.7), the econNMPC technique tries to minimize the specific energy con- sumption and simultaneously tries to minimize transition time to reach a steady state from the current state. Throughout this chapter, sptNMPC and econNMPC denote for the NMPC formulations (6.6), with setpoint tracking objective function (6.6a), and with the economically-oriented objective function (6.7), respectively. 6.4.2 Closed-loop Simulation The objective of this simulation study is to show the economical benefits of the econNMPC technique with the presence of disturbances. Changes in the chip bulk density and chip solid content are considered as the main disturbances, and ad- ditionally, random disturbances on the state variables are also considered in the the TMP process (6.1). As discussed in Example 5.4.1, changes in the chip bulk density dc and chip solid content sc in (5.1) in turn cause some variations in the production rate, motor loads, consistencies and specific energy (6.8). Table 6.1 shows the values of the chip bulk density and chip solid content from their nominal values during the closed-loop simulation. Time, s 0 -60 60 - 120 120 - 180 180 - 240 240 - 300 Chip bulk density 95% 110% 100% 110% 90% Chip solid content 100% 100 % 90 % 90% 90% Table 6.1: The chip bulk density and chip solid content changes during the closed- loop simulation of the TMP process in Scheme A The sptNMPC and econNMPC techniques are both implemented with the same set of disturbances. The sampling time is 1s, and the prediction and control hori- 87 6.5. Scheme B: Improved TMP Process with Single-stage of HC and Multi-stage LC Refiners zons are set as NA = 30 for both the sptNMPC and econNMPC techniques. The controlled and manipulated weighting matrices are kept as Qξ A = diag([0.01 0.5 1 0.5 1 1 1 0.5 1]T × 10−1) and QvA = diag([0.8 8 0.8 8 0.8 0.8 0.8 8]T × 10−1). αAtmp = 1000 and β Atmp = 1. The closed-loop simulation results are shown in Fig. 6.2 - 6.5. As evident from the result in Fig. 6.4, there is 188 kWh/tonne reduc- tion in average specific energy, i.e. about 13%, in the system with econNMPC technique when disturbances are present, compared to the system with sptNMPC technique. Both sptNMPC and econNMPC techniques were able to keep the pulp quality variables: long fiber content, and shive content, within the operating lim- its of (50−75%), and (0−2.5%), respectively, while keeping final freeness after the tertiary LC refining at the desired value of 106.5 (CSF). Here, the pulp quality variables are considered being measured once every second just after the tertiary LC refiner. 6.5 Scheme B: Improved TMP Process with Single-stage of HC and Multi-stage LC Refiners In this section, an improved TMP process presented in [1] is studied for optimal operation, where the secondary stage refining in the conventional TMP process, see Fig. 6.1, is replaced by multi-stages of LC refining after the latency chest as shown in Fig. 6.6. The process is dynamically optimized when the disturbances are present. As discussed in Section 6.3, the improved TMP process is modeled for the NMPC formulation. In Scheme B, the two-stages of LC refining are found to be optimal at steady state to replace the secondary stage HC refiner of the TMP process discussed in Scheme A [1]. Hence in this work, two LC refiners are placed after the latency chest as shown in Fig. 6.6. For this process, the following general nonlinear model can be expressed. zk+1 = f Btmp = ABzk +hB(zk,vk) (6.9a) 0 = gB(zk,yk) (6.9b) 88 6.5. Scheme B: Improved TMP Process with Single-stage of HC and Multi-stage LC Refiners where z, y, and v can be given by z =   Production rate Primary refiner motor load Primary refiner consistency Consistency after latency chest CSF after latency chest Tertiary refiner 1 motor load CSF after tertiary refiner 1 Tertiary refiner 2 motor load CSF after tertiary refiner 2   , y =   CSF after primary refining LFC after primary refining SC after primary refining LFC after latency chest SC after latency chest LFC after tertiary refiner 1 SC after tertiary refiner 1 LFC after tertiary refiner 2 SC after tertiary refiner 2   , and v =   Chip-transfer screw speed Primary refiner plate gap Primary refiner dilution water flow rate Latency chest temperature Latency chest dilution water flow rate Tertiary refiner 1 plate gap Tertiary refiner 2 plate gap   . (6.10) 6.5.1 NMPC for Setpoint Tracking and Economical Operation Similar to the optimization problem in Scheme A, the dynamic optimization prob- lem is formulated for the NMPC at sample time k (6.11). Here all the variables are corresponding to the TMP process in Scheme B, and have been defined similar to Scheme A. First, the setpoint tracking problem is solved. Then econNMPC is formulated with the economic objective function 6.12. For Scheme B, ξ B ≡ z. In oder to compare performances and specific energy savings of the TMP processes in Schemes A and B, the final freeness after the tertiary LC refiner 2, is controlled 89 6.5. Scheme B: Improved TMP Process with Single-stage of HC and Multi-stage LC Refiners to that of the same desired value of the TMP process in Scheme A. min v JBspt , JBspt := NB ∑ i=1 ‖ξ Bi −ξ Bsk+i‖2Qξ B +‖vi−1−uBsk+i−1‖2QvB (6.11a) s. t. : z0 = xk, (6.11b) zi = f Btmp(zi−1,vi−1), (6.11c) 0 = gBtmp(zi−1,yi−1), (6.11d) zmin ≤ zi ≤ zmax, (6.11e) ymin ≤ yi ≤ ymax, (6.11f) umin ≤ vi−1 ≤ umax, i = 1, . . . ,NB (6.11g) In the econNMPC formulation, the following objective function (6.12) is mini- mized. JBecon = αBtmp NB ∑ i=1 Specific Energy+β Atmp NB ∑ i=1 ‖ξi−ξi+1‖2Qξ A +‖vi−1− vi‖2QvA , (6.12) where Specific Energy is defined as Specific Energy = Total motor load Production rate MW/(tonnes/day) = Primary refiner motor load + Tertiary refiner 1 motor load + Tertiary refiner 2 motor load Production rate . (6.13) αBtmp and β Btmp are positive weights. 6.5.2 Closed-loop Simulation In this simulation study, the sptNMPC and econNMPC techniques are formulated for the TMP process in Scheme B discussed in Section 6.6. The simulation set-up is the same as that in Scheme A discussed in Section 6.4.2, and the disturbances presented in Table 6.1 are placed during the closed-loop simulation. The sam- pling time is 1s, and the prediction and control horizons are set as NB = 30 for 90 6.6. Conclusion both the sptNMPC and econNMPC techniques. The controlled and manipulated weighting matrices are kept as Qξ B = diag([0.01 0.5 1 1 1 0.5 1 0.5 1]T × 10−1) and QvB = diag([0.8 8 0.8 0.8 0.8 8 8]T ×10−1). αBtmp = 1000 and β Btmp = 1. The closed-loop simulation results are shown in Fig. 6.7 - 6.10. There is 88 kWh/tonne reduction in average specific energy, i.e. about 7.5%, in the system with econN- MPC technique when disturbances are present, compared to the system with spt- NMPC technique. Similar to Scheme A 6.4.2, both the sptNMPC and econNMPC techniques were able to maintain the pulp quality variables: long fiber content, and shive content, within the operating limits of (50− 75%), and (0− 2.5%), respec- tively, while keeping final freeness after the tertiary LC refining at the desired value of 106.5 (CSF). 6.6 Conclusion In this chapter, optimal operations of two TMP process schemes have been studied. The Scheme A considers a conventional TMP process while Scheme B considers an improved TMP process with multi-stages of LC refining. The dynamic opti- mizations are performed within the NMPC framework for both the processes. Sig- nificant savings in specific energy are noticed when the disturbances are present, compared to setpoint tracking NMPC formulation, with the econNMPC technique for both the Schemes, while respecting the pulp quality requirements. Fig. 6.11 shows the specific energy comparison for both the Schemes with sptNMPC and econNMPC techniques. The average reduction in specific energy from the system with sptNMPC in Scheme A to the system with sptNMPC in Scheme B is around 252 kWh/tonne, about 18%. The average reduction is around 340 kWh/tonne, or 24%, from the system with sptNMPC in Scheme A to the system with econNMPC in Scheme B. Note that in [1], it showed that around 20% of specific energy reduc- tion can be achieved with an improved TMP process at the steady state, compare to a convention TMP process in Scheme A. Finally, this simulation study concludes that multi-stage LC refining technology can save specific energy of approximately 18% with existing setpoint tracking control when disturbances are present. Fur- ther, approximately 6% reduction can be achieved with the proposed econNMPC technique. 91 6.6. Conclusion 0 50 100 150 200 250 300 20 22 24 26 S cr ew s p ee d R P M Time, s 0 50 100 150 200 250 300 1.5 2 P ri m ar y g ap m m 1 0 Time, s 0 50 100 150 200 250 300 249.5 250 250.5 P ri m ar y d il u ti o n f lo w , l/ m in Time, s S ec o n d ar y g ap m m 0.5 0 1.5 2 0 50 100 150 200 250 300 Time, s S ec o n d ar y d il u ti o n fl o w , l/ m in 0 50 100 150 200 250 300 Time, s 99.6 100 100.4 0.6 T er ti ar y g ap m m 0.3 0.4 0.5 0 50 100 150 200 250 300 Time, s econNMPC sptNMPC set-point Figure 6.2: The manipulated input variables of the TMP process in Scheme A. 92 6.6. Conclusion 39 41 43 45 P ri m ar y co n si st en cy , % 200 240 280 P ro d u ct io n r at e t o n n e/ d ay 0 50 100 150 200 250 300 Time, s 7.5 8.5 9.5 P ri m ar y m o to r lo ad , M W 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 4.5 5.5 6.5 S ec o n d ar y m o to r l o ad , M W 31 33 35 37 S ec o n d ar y co n si st en cy , % 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s T er ti ar y m o to r lo ad , M W 0 50 100 150 200 250 300 Time, s 1 1.3 1.6 4.15 4.25 4.35 4.45 T er ti ar y co n si st en cy , % 0 50 100 150 200 250 300 Time, s econNMPC sptNMPC set-point Figure 6.3: The controlled state variables of the TMP process in Scheme A. 93 6.6. Conclusion 500 600 700 800 P ri m ar y r ef in er s p ec if ic en er g y , k W h /t 400 500 600 700 S ec o n d ar y r ef in er sp ec if ic e n er g y , k W h /t 100 140 180 T er ti ar y r ef in er sp ec if ic e n er g y , k W h /t 1000 1200 1400 1600 T o ta l sp ec if ic e n er g y k W h /t econNMPC sptNMPC 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s Figure 6.4: The specific energy consumptions of the primary, secondary, and ter- tiary refiners of TMP process in Scheme A. 94 6.6. Conclusion 58 60 62 64 66 T er ti ar y L F C % 0.65 0.75 0.85 0.95 T te rt ia ry S C % econNMPC sptNMPC set-point T er ti ar y f re en es s C S F m l 101 103 105 107 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s Figure 6.5: The pulp quality variables after the tertiary refining of the TMP process in Scheme A. 95 6.6. Conclusion Screening, Steaming Washing, Preheating HC refining Wood chips Pulp 2rd stage LC refining Screening, Reject refining Latency removal Dewatering LC refining 1st stage LC refining Figure 6.6: Scheme B: basic schematic of an improved TMP process with a two stages of LC refining 96 6.6. Conclusion 20 24 28 S cr ew s p ee d R P M 0 0.5 1 1.5 P ri m ar y g ap m m 249.8 250.2 250.6 P ri m ar y d il u ti o n f lo w , l/ m in econNMPC sptNMPC set-point 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0.32 0.36 0.40 0.44 T er ti ar y r ef in er 1 g ap , m m 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0.32 0.36 0.4 0.44 0.48 T er ti ar y r ef in er 2 g ap , m m 0 50 100 150 200 250 300 Time, s Figure 6.7: The manipulated input variables of the TMP process in Scheme B. 97 6.6. Conclusion 200 240 280 P ro d u ct io n r at e t o n n e/ d ay 7 8 9 10 P ri m ar y m o to r lo ad , M W 39 41 43 45 P ri m ar y c o n si st en cy , % 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s T er ti ar y r ef in er 1 m o to r lo ad , M W 1.6 1.7 1.8 1.9 2.0 2.1 0 50 100 150 200 250 300 Time, s T er ti ar y r ef in er 2 m o to r lo ad , M W 0 50 100 150 200 250 300 Time, s 1.4 1.6 1.8 2.0 econNMPC sptNMPC set-point Figure 6.8: The controlled state variables of the TMP process in Scheme B. 98 6.6. Conclusion 600 800 1000 140 180 220 100 150 200 250 900 1100 1300 P ri m ar y r ef in er s p ec if ic en er g y , k W h /t T er ti ar y r ef in er 1 sp ec if ic e n er g y , k W h /t T er ti ar y r ef in er 2 sp ec if ic e n er g y , k W h /t T o ta l sp ec if ic e n er g y k W h /t 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s econNMPC sptNMPC Figure 6.9: The specific energy consumptions of the primary refiner, tertiary refiner 1, and tertiary refiner 2 of TMP process in Scheme B. 99 6.6. Conclusion 100 104 108 112 T er ti ar y r ef in er 2 L F C , % T te rt ia ry r ef in er 2 S C , % T er ti ar y r ef in er 2 fr ee n es s C S F , m l 0 50 100 150 200 250 300 Time, s 67.5 68.5 69.5 70.5 71.5 1.9 2 2.1 2.2 0 50 100 150 200 250 300 Time, s 0 50 100 150 200 250 300 Time, s econNMPC sptNMPC set-point Figure 6.10: The pulp quality variables after the tertiary refining of the TMP pro- cess in Scheme B. 100 6.6. Conclusion 1000 1100 1200 1300 1400 1500 1600 S p ec if ic e n er g y , k W h /t 900 sptNMPC, Scheme A econNMPC, Scheme A sptNMPC, Scheme B econNMPC, Scheme B 0 50 100 150 200 250 300 Time, s Figure 6.11: Specific energies of both Schemes A and B with sptNMPC and econ- NMPC techniques. 101 Chapter 7 Conclusions In this work, we first studied closed-loop system identification techniques which can be applied to a Linear Time-Invariant (LTI) system controlled by a constrained Model Predictive Control (MPC) law. We investigated two techniques: 1) direct identification technique using Prediction Error Method (PEM), and 2) joint input- output technique using the PEM with the Long Range Predictive Identification (LRPI) criteria. In the direct identification technique, we studied how to use nor- mal operational data for identification task without preforming any identification experiments. In the joint input-output identification technique, we studied LRPI criteria, a more meaningful identification criteria than the one-step ahead criteria in the context of identification for control using MPC. These techniques can be ap- plied in the advanced control systems of the existing Thermo-Mechanical Pulping (TMP) processes to re-identify process models for improved performances, in par- ticular, their low-level control systems where linear process model is being used in the MPC formulation. In the second part of the thesis, we studied economical operation of the two- stage TMP process. A dynamic optimization problem is formulated to simultane- ously regulate and optimize the TMP process when disturbances are present. The proposed economically oriented Nonlinear Model Predictive Control (econNMPC) minimizes the total specific energy consumption of both the primary and secondary refiners, respecting process constraints of such as motor loads, plate gaps, and pulp properties. In the simulation study, we first implemented setpoint tracking NMPC technique for the TMP process, where the controller task was to track the optimal steady state setpoint given by upper-level Real-Time Optimization (RTO) layer. Next, the econNMPC technique was implemented. Both the control schemes showed similar performances while the closed-loop system with econN- MPC technique showed around 12% reduction in total specific energy compared 102 Chapter 7. Conclusions to the closed-loop system with setpoint tracking NMPC technique. In the third part of the thesis, an improved TMP process is studied for its op- timal operation. The improved TMP process has a High Consistency (HC) refiner and two Low Consistency (LC) refiners. These refiners are dynamically coordi- nated using econNMPC technique. In a previous study [1], it is shown that around 20% specific energy reduction can be achieved with an improved TMP process at the steady state, compared to a convention TMP process consisting of two HC refiners and one LC refiner. In this study, we also showed that additional approxi- mately 6% of specific energy reduction can be realized when the process is dynam- ically optimized. Future Work We have assumed that the state variables are measured for all NMPC formulations. However in practice, all the state measurements of TMP processes are not readily available to integrate with real-time controllers. For example, outlet consistencies of the refiners are not measured, but need to be estimated instantaneously. On the other hand, pulp qualities are usually measured just after the latency chest. There- fore, integration of those measurements with the NMPC scheme would require an estimation technique. However, this would increase complexity of the problem and require more computational power, and also would need to address theoretical is- sues such as stability of the overall scheme. Estimation methods such as Moving Horizon Estimation (MHE) technique or Extend Kalman Filter (EKF) would be more suitable for the TMP processes. Incorporation of such an estimation tech- nique with NMPC algorithms for TMP processes is therefore our future research interest. At the same time, we are in the process of first preparing mill trials and then implementing the proposed techniques on a real HC and LC TMP plant. The implementation of the proposed optimal control strategy is based on the state space modeling of the TMP refiners. 103 Bibliography [1] A. Luukkonen. Development of a methodology to optimize low consis- tency refining of mechanical pulp. PhD thesis, The University of British Columbia, June 2011. [2] J. Sundholm. Mechanical Pulping. Papermaking Science and Technology. Fapet Oy, 1999. [3] L. Euhus, B. C. Strand, D. Shearin, G. Fralic, and M. Edvall. Optimiza- tion of LWC TMP using advanced quality control. In 2003 International Mechanical Pulping Conference, pages 375–378, Jun 2003. [4] O. Johansson, M. Jackson, and N.W. Wild. Three steps to improved TMP operating efficiency. In International Mechanical Pulping Conference 2007, TAPPI, volume 1, pages 356–369. TAPPI, 2007. [5] B. C. Strand and Greg Fralic. Economic benefits from advanced quality control of TMP mills. In Control Systems 2000 Conference, pages 11–15, May 2000. [6] H. Du. Multivariable Predictive Control of a TMP Plant. PhD thesis, Uni- versity of British Columbia, October 1998. [7] D. D. Ruscio. Topics in Model Based Control with Application to the Thermo Mechanical Pulping Process. PhD thesis, The Norwegian Institute of Technology, 1993. [8] M. S. Sidhu, R. V. Fleet, M. R. Dion, D. W. Anderson, and B. W. Weger. Modeling and advanced control of TMP refiner system. In Control Systems 2004 Conference, Jun 2004. 104 Bibliography [9] F. Rosenqvist, D. Berg, A. Karlstrom, K. Eriksson, and C. Breitholtz. In- ternal interconnections in tmp processes. In IEEE Conference on Control Applications, pages 1010–1015, Sep 2002. [10] G. A. Dumont and K. J. Åström. Wood chip refiner control. IEEE Control Systems Magazine, 8(2):38–43, 1988. [11] H. M. Schwartz, G. Chang, and Y. Liu. Method of modeling, predicting and controlling tmp pulp properties. In IEEE Conference on Control Applica- tions, pages 846–851, Sep 1996. [12] I. Lama, M. Perrier, and P. Stuart. Controllability analysis of a tmp- newsprint refining process. Pulp and Paper Canada, 107(10):44–48, 2006. [13] A. Karlstrom, K. Eriksson, D. Sikter, and M. Gustavsson. Refining models for control purposes. Nordic Pulp and Paper Research Journal, 23(1):129– 138, 2008. [14] M. Elsinga. TMP optimization using multivariable analysis. IEEE Transac- tions on Industry Applications, (3):893–898, May/June 2003. [15] A. M. Brdys and P. Jatjewski. Iterative Agorithms for Multilayer Optimizing Control. Imperial College Press, 2005. [16] J. B. Rawlings, D. Bonne, J. B. Jorgensen, A. N. Venkat, and S. B. Jor- gensen. Unreachable setpoints in model predictive control. Automatic Con- trol, IEEE Transactions on, 53(9):2209–2215, oct. 2008. [17] R. Huang. Nonlinear Model Predictive Control and Dynamic Real Time Optimization for Large-scale Processes. PhD thesis, Carnegie Mellon Uni- versity, December 2010. [18] M. Sabourin, J. Aichinger, and R. Eksteen. Minimizing tmp energy con- sumption using a combination of chip pre-treatment, rts and multiple stage low consistency refining. In International Mechanical Pulping Conference 2007, TAPPI, pages 839–893, May 2007. 105 Bibliography [19] S. J. Qin and T. A. Badgwell. A survey of industrial model predictive control technology. Control Engineering Practice, 11(7):733–764, July 2003. [20] A. I. Propoi. Use of linear programming methods for synthesizing sampled- data automatic systems. Automation and Remote Control, 24(7):837–844, 1963. [21] J. Richalet, A. Rault, J. L. Testud, and J. Papon. Model predictive heuristic control: Applications to industrial processes. Automatica, 14(5):413–428, 1978. [22] C. R. Cutler and B. L. Ramaker. Dynamic matrix control-a computer control algorithm. In AICHE national meeting, April 1979. [23] D. W. Clarke, C. Mohtadi, and P. S. Tuffs. Generalized predictive controlpart i. the basic algorithm. Automatica, 23(2):137–148, 1987. [24] D. W. Clarke, C. Mohtadi, and P. S. Tuffs. Generalized predictive controlpart i. the basic algorithm. Automatica, 23(2):149–160, 1987. [25] C. E. Garcı́a, D. M. Prett, and M. Morari. Model predictive control: Theory and practice-a survey. Automatica, 25(3):335–348, 1989. [26] R. Findeisen and F. Allgöwer. An introduction to nonlinear model predic- tive control. In C.W. Scherer and J.M. Schumacher, editors, The Impact of Optimization in Control, pages 3.1–3.45. Dutch Institute of Systems and Control, DISC, 2001. [27] P. Tatjewski. Advanced Control of Industrial Processes. Advances in Indus- trial Control. Springer London, 2007. [28] D. Q. Mayne. Nonlinear model predictive control: Challenges and opportu- nities. In Frank Allgöwer and Alex Zheng, editors, Nonlinear Model Pre- dictive Control, pages 23–44. Birkhäuser, 2000. [29] M. Diehl, R. Amrit, and J. B. Rawlings. A Lyapunov function for economic optimizing model predictive control. IEEE Transactions on Automatic Con- trol, 56(3):703–707, 2011. 106 Bibliography [30] H. K. Khalil. Nonlinear Systems. Prentice Hall, 2nd edition, 1996. [31] E. G. Gilbert and K.T. Tan. Linear systems with state and control constraints: the theory and application of maximal output admissible sets. Automatic Control, IEEE Transactions on, 36(9):1008–1020, sep 1991. [32] L. Magni and R. Scattolini. Robustness and robust design of mpc for nonlin- ear discrete-time systems. In Rolf Findeisen, Frank Allgöwer, and Lorenz Biegler, editors, Assessment and Future Directions of Nonlinear Model Pre- dictive Control, volume 358 of Lecture Notes in Control and Information Sciences, pages 239–254. Springer Berlin / Heidelberg, 2007. [33] M. Lazar, M. Heemels, A. Bemporad, and S. Weil. On the stability and ro- bustness of non-smooth nonlinear model predictive control. In in Workshop on Assessment and Future Directions of NMPC, Freudenstadt-Lauterbad, pages 327–334, 2005. [34] S. S. Keerthi and E. G. Gilbert. Optimal infinite-horizon feedback laws for a general class of constrained discrete-time systems: Stability and moving- horizon approximations. Journal of Optimization Theory and Applications, 57(2):265–293, 1988. [35] D. Q. Mayne and H. Michalska. Receding horizon control of nonlinear systems. IEEE Transactions on Automatic Control, 35(7):814–824, 1990. [36] H. Michalska and D. Q. Mayne. Robust receding horizon control of con- strained nonlinear systems. IEEE Transactions on Automatic Control, 38(11):1623–1633, 1993. [37] R. R. Bitmead, M. Gevers, and V. Wertz. Adaptive Optimal Control. Prentice Hall, 1990. [38] T. Parisini and R. Zoppoli. A receding-horizon regulator for nonlinear systems and a neural approximation. Automatica, 31:1443–1451, October 1995. 107 Bibliography [39] H. Chen and F. Allgöwer. Quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica, 34(10):1205–1217, 1998. [40] G. De Nicolao, L. Magni, and R. Scattolini. Stabilizing receding-horizon control of nonlinear time-varying systems. IEEE Transactions on Automatic Control, 43(7):1030–1036, 1998. [41] M. Alamir and G. Bornard. On the stability of receding horizon control of nonlinear discrete-time systems. Systems and Control Letters, 23(4):291– 296, 1994. [42] E. S. Meadows, M. A. Henson, J. W. Eaton, and J. B. Rawlings. Receding horizon control and discontinuous state feedback stabilization. International Journal of Control, 62(5):1217–1229, 1995. [43] Z. P. Jiang and Y Wang. Input-to-state stability for discrete-time nonlinear systems. Automatica, 37(6):857–869, June 2001. [44] D. Limon, T. Alamo, D. M. Raimondo, D. de la Peña, J. Bravo, A. Fer- ramosca, and E. F. Camacho. Input-to-state stability: A unifying framework for robust model predictive control. In Lalo Magni, Davide Raimondo, and Frank Allgöwer, editors, Nonlinear Model Predictive Control, volume 384 of Lecture Notes in Control and Information Sciences, pages 1–26. Springer Berlin / Heidelberg, 2009. [45] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert. Constrained model predictive control: Stability and optimality. Automatica, 36(6):789– 814, June 2000. [46] L. Magni, G. De Nicolao, R. Scattolini, and F. Allgöwer. Robust model predictive control for nonlinear discrete-time systems. International Journal of Robust and Nonlinear Control, 13(3-4):229–246, 2003. [47] D. M. Raimondo, D. Limon, M. Lazar, L. Magni, and E. F. Camacho. Min- max model predictive control of nonlinear systems: A unifying overview on stability. European Journal of Control, 15(1):5–21, 2009. 108 Bibliography [48] D. L. Marruedo, T. Alamo, and E. F. Camacho. Input-to-state stable mpc for constrained discrete-time nonlinear systems with bounded additive un- certainties. In Decision and Control, 2002, Proceedings of the 41st IEEE Conference on, volume 4, pages 4619–4624, dec. 2002. [49] D. Q. Mayne, M. M. Seron, and S. V. Raković. Robust model predictive con- trol of constrained linear systems with bounded disturbances. Automatica, 41(2):219–224, 2005. [50] S. V. Raković, A. R. Teel, D. Q. Mayne, and A. Astolfi. Simple robust control invariant tubes for some classes of nonlinear discrete time systems. pages 6397–6402, San Diego, CA, United states, 2006. [51] A. Bemporad and M. Morari. Robust model predictive control: A survey. In A. Garulli and A. Tesi, editors, Robustness in identification and control, volume 245 of Lecture Notes in Control and Information Sciences, pages 207–226. Springer Berlin / Heidelberg, 1999. [52] J. B. Rawlings and D. Q. Mayne. Model Predictive Control: Theory and Design. Nob Hill Publishing, 2009. [53] S. J. Wright. Applying new optimization algorithms to model predictive control. In Proceedings of the CPC-V, pages 147–155. CACHE Publica- tions, 1997. [54] J. M. Maciejowski. Predictive Control with Constraints. Prentice Hall, 2002. [55] E. F. Camacho and C. Bordons. Model Predictive Control. Advanced Text- books in Control and Signal Processing. Springer, 2 edition, 2004. [56] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos. The explicit linear quadratic regulator for constrained systems. Automatica, 38:3–20, January 2002. [57] M. Diehl, H. G. Bock, and J. P. Schlöder. A real-time iteration scheme for nonlinear optimization in optimal feedback control. SIAM Juournal on Cotrol and Optimization, 43(5):1714–1736, 2005. 109 Bibliography [58] V. M. Zavala and L. T. Biegler. The advanced step nmpc controller: Opti- mality, stability and robustness. Automatica, 45(1):86–93, January 2009. [59] A. Wächter and L. T. Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathe- matical Programming, 106(1):25–57, May 2006. [60] V. M. Zavala. Computational Strategies for the Optimal Operation of Large- Scale Chemical Processes. PhD thesis, Carnegie Mellon University, 2008. [61] D. E. Kassmann, T. A. Badgwell, and R. B. Hawkins. Robust steady-state target calculation for model predictive control. AIChE Journal, 46(5):1007– 1024, May 2000. [62] A. C. Zanin, M. T. de Gouvea, and D. Odloak. Integrating real-time op- timization into the model predictive controller of the fcc system. Control Engineering Practice, 10(8):819–31, 2002. [63] L. T. Biegler. Efficient Solution of Dynamic Optimization and NMPC Prob- lems. Progress in Systems and Control Theory. Birkhäuser, 2000. [64] D. Q. Mayne. Nonlinear model predictive control: An assessment. Number 316, pages 217–217, 1997. [65] V. M. Zavala and L. T. Biegler. Optimization-based strategies for the oper- ation of low-density polyethylene tubular reactors: nonlinear model predic- tive control. Computers & Chemical Engineering, 33(10):1735–46, October 2009. [66] V. Adetola and M. Guay. Integration of real-time optimization and model predictive control. J. Process Control, 20(2):125–133, February 2010. [67] E. Harinath, L. T. Biegler, and G. A. Dumont. Control and optimization strategies for thermo-mechanical pulping processes: Nonlinear model pre- dictive control. Journal of Process Control, 21(4):519–528, April 2011. 110 Bibliography [68] R. Huang, V. M. Zavala, and L. T. Biegler. Advanced step nonlinear model predictive control for air separation units. J. Process Control, 19(4):678– 685, April 2009. [69] R. Huang, E. Harinath, and L. T. Biegler. Lyapunov stability of econom- ically oriented NMPC for cyclic processes. Journal of Process Control, 21(4):501–509, April 2011. [70] E. M. B. Aske, S. Strand, and S. Skogestad. Coordinator mpc for maximiz- ing plant throughput. Computers and Chemical Engineering, 32(1-2):195– 204, 2008. [71] S. Engell. Feedback control for optimal process operation. Journal of Pro- cess Control, 17(3):203–219, 2007. [72] J. Kadam and W. Marquardt. Integration of economical optimization and control for intentionally transient process operation. In Rolf Findeisen, Frank Allgöwer, and Lorenz Biegler, editors, Assessment and Future Direc- tions of Nonlinear Model Predictive Control, volume 358 of Lecture Notes in Control and Information Sciences, pages 419–434. Springer Berlin / Hei- delberg, 2007. [73] J. B. Rawlings and R. Amrit. Optimizing process economic performance us- ing model predictive control. In Lalo Magni, Davide Raimondo, and Frank Allgöwer, editors, Nonlinear Model Predictive Control, volume 384 of Lec- ture Notes in Control and Information Sciences, pages 119–138. Springer Berlin / Heidelberg, 2009. [74] R. Huang, L. T. Biegler, and E. Harinath. Robust stability of economically oriented infinite horizon NMPC that include cyclic processes. Journal of Process Control, 22(1):51–59, January 2012. [75] M. Heidarinejad, J. Liu, and P. D. Christofides. Economic model predictive control of nonlinear process systems using lyapunov techniques. AIChE Journal, 2011. 111 Bibliography [76] L. Grüne. Optimal invariance via receding horizon control. Submitted in 50th IEEE Conference on Decision and Control and European Control Con- ference, 2011. [77] L. Würth, J. B. Rawlings, and W. Marquardt. Economic dynamic real-time optimization and nonlinear model predictive control on infinite hoirzons, in international symposium on advanced control of chemical process. volume 7 of 7th IFAC International Symposium on Advanced Control of Chemical Processes, 2009. [78] J. Nocedal and S. J. Wright. Numerical Optimization. Operations Research. Springer, 2nd edition, 2007. [79] H. Hjalmarsson. From experiment design to closed-loop control. Automat- ica, 41(3):393–438, March 2005. [80] J. S. Conner and D. E. Seborg. Assessing the need for process re- identification. Industrial & Engineering Chemistry Research, 44(8):2767– 2775, 2005. [81] O. A. Z. Sotomayor, D. Odloak, and L. F. L. Moro. Closed-loop model re-identification of processes under mpc with zone control. Control Engi- neering Practice, 17(5):551–563, 2009. [82] M. A. Hussain. Review of the applications of neural networks in chemical process control – simulation and online implementation. Artificial Intelli- gence in Engineering, 13(1):55–68, 1999. [83] U. Forssell and L. Ljung. Closed-loop identification revisited. Automatica, 35(7):1215–1241, Jul 1999. [84] L. Ljung. System Indentification: Theory for the user. Prentice-Hall, 2 edition, 1999. [85] T. Söderström and P. Stoica. System Identification. Prentice Hall Interna- tionl, 1989. 112 Bibliography [86] P. M. J. Van den Hof and R. A. de Callafon. Multivariable closed-loop identification: from indirect identification to dual-youla parametrization. In Decision and Control, 1996., Proceedings of the 35th IEEE, volume 2, pages 1397–1402 vol.2, dec 1996. [87] U. Forssell and L. Ljung. Projection method for closed-loop identification. IEEE Transactions on Automatic Control, 45(11):2101–2106, Nov 2000. [88] J. Yan, E. Harinath, and G. A. Dumont. Closed-loop identification for model predictive control: Direct method. In Proceedings of the IEEE Conference on Decision and Control, pages 2592–2597, December 2009. [89] M. Gevers, A. Bazanella, and L. Mikovic. Informative data: How to get just sufficiently rich? pages 1962–1967, 2008. [90] X. Bombois, G. Scorletti, M. Gevers, P. M. J. Van den Hof, and R. Hilde- brand. Least costly identification experiment for control. Automatica, 42(10):1651–1662, 2006. [91] M. J. Van Den Hof and R. J. P. Schrama. Indirect method for transfer func- tion estimation from closed loop data. Automatica, 29(6):1523–1527, Nov 1993. [92] U. Forssell. Closed-loop Identification: Methods, Theory, and Applications. Linköping studies in science and technology. thesis no 566, 1999. [93] D. S. Shook, C. Mohtadi, and S. L. Shah. Identification for long-range predictive control. IEE Proceedings D: Control Theory and Applications, 138(1):75–84, 1991. [94] R. B. Gopaluni, R. S. Patwardhan, and S. L. Shah. The nature of data pre- filters in MPC relevant identification - open- and closed-loop issues. Auto- matica, 39(9):1617–1626, 2003. Bias distribution;. [95] R. B. Gopaluni, R. S. Patwardhan, and S. L. Shah. MPC rele- vant identification–tuning the noise model. Journal of Process Control, 14(6):699–714, 2004. 113 Bibliography [96] M. Morari and J. H. Lee. Model predictive control: past, present and future. Computers and Chemical Engineering, 23(4-5):667–682, May 1999. [97] D. D. Ruscio. Model predictive control and identification: a linear state space model approach. In Proceedings of the 36th IEEE Conference on Decision and Control. IEEE, December 1997. [98] H. Du, G. A. Dumont, and Y. Fu. Nonlinear control of a wood chip refiner. In Proceedings of the 4th IEEE Conference on Control Applications, pages 1065–6. IEEE, 1995. [99] S. Kidd. Implementation of a tmp advanced quality control system at a newsprint manufacturing plant. Technical report, Department of Energy, Washington, DC, 2005. [100] M. Illikainen. Mechanisms of thermomechanical pulp refining. PhD thesis, University of Oulu, October 2008. [101] X. Qian and P. Tessier. A mechanistic model for predicting pulp properties from refiner operating conditions. Tappi journal, 78(78):215–222, April 1995. [102] L. O. Santos, P. A. F. N. A. Afonso, J. A. A. M Castro, N. M. C. Oliveira, and L. T. Biegler. Online implementation of nonlinear mpc: an experimental case study. Control Engineering Practice, 9(8):847–857, August 2001. [103] W. H Chen, D. J. Ballance, and J. O’Reilly. Model predictive control of nonlinear systems: Computational burden and stability. IEE Proc Control Theory Applications, 147(4):387–394, July 2000. [104] R. Findeisen and F. Allgöwer. Computational delay in nonlinear model predictive control. In Int. Symp. Adv. Control of Chemical Processes, AD- CHEM’03, pages 427–432, 2004. [105] S. Skogestad. Self-optimizing control: The missing link between steady- state optimization and control. Computers and Chemical Engineering, 24(2- 7):569–575, July 2000. 114 Bibliography [106] B. C. Strand, J. Straight, and D. Norris. Factors affecting energy reduction of newsprint from southern pine tmp. [107] G. Mitchell, R. Musselman, and M. Sabourin. Third stage low consistency refining of tmp. In 52nd Appita general pulping conference, pages 475–483, 1998. [108] R. Musselman, D. Letarte, and R. Sinard. Third stage low consistency refin- ing of tmp for energy saving and quality enhancement. In PIRA 4th Inter- national Refining Conference, 1997. [109] J. Gao, C. P. J. Bennington, D. M. Martinez, and J. A. Olson. Latency re- moval in mechanical pulping processes. In PAPTAC Annual Meeting, pages 1293–1325. CD–ROM, 2012. 115
Thesis/Dissertation
eng
Electrical and Computer Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Identification, control and optimization strategies for thermo-mechanical pulping (TMP) processes
Text