Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Stochastic averaging level control and its application to broke management in paper machines Ogawa, Shiro 2003

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

Item Metadata


831-ubc_2004-902439.pdf [ 6.1MB ]
JSON: 831-1.0065478.json
JSON-LD: 831-1.0065478-ld.json
RDF/XML (Pretty): 831-1.0065478-rdf.xml
RDF/JSON: 831-1.0065478-rdf.json
Turtle: 831-1.0065478-turtle.txt
N-Triples: 831-1.0065478-rdf-ntriples.txt
Original Record: 831-1.0065478-source.json
Full Text

Full Text

Stochastic Averaging Level Control and Its Application to Broke Management in Paper Machines by SHIRO OGAWA B. Eng., The University of Tokyo, 1967 M . Eng., The University of Tokyo, 1969 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of Electrical and Computer Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 12, 2003 ©Shiro Ogawa, 2003 Abstract Averaging level control refers to liquid level control of storage tanks, where the objective is to keep the outlet flow u(t) as smooth as possible against the fluctuating inlet flow, while at the same time keeping the tank level y{t) within high and low limits. The thesis treats the stochastic averaging level control problem, where the input disturbance is a stochastic process. The problem is formulated to minimize a weighted sum of Var[u(r)] and Var[«(/)] subject to the target Var[y(r)]. The state-space linear quadratic optimal control method is used, resulting in a linear state feedback controller. When the input disturbance is modelled as the output of a first-order low-pass filter driven by white noise, the optimal controller is a phase-lag network. Broke storage tank level control is important in stabilizing the paper machine wet end. It is treated as a special type of averaging level control, where the input disturbance Fb(t) is a two-state continuous-time Markov process. The spectrum of Fb(t) is obtained and the linear optimal con-troller is designed with the same methodology as for the general averaging level control problem. Taking this very specific nature of Fb(t), a new nonlinear control scheme called the minimum overflow probability controller (MOPC) is designed, and tested against data collected from a paper machine. The MOPC performs better than the optimal linear controller and manual control. A new theorem on the state probability distribution of a continuous-time Markov jump system is presented, which leads to new methods for evaluating the mean and the variance of the state of a linear jump system, and a new reliable numerical method to calculate the state distributions of jump systems. These results are utilized to evaluate the overflow probabilities of controllers. ii Contents Abstract ,,. ii List of Tables vii List of Figures ix Acknowledgements x 1 Introduction and Overview 1 1.1 Background and Motivations 1 1.2 Overview and Structure of the Thesis 2 1.2.1 Format and Notation 4 1.3 Contributions of the Thesis 4 2 Averaging Level Control 7 2.1 Introduction 7 2.2 Mathematical Models 8 2.2.1 Process Model 8 2.2.2 Disturbance Model . . ! 9 2.2.3 Flow Smoothness/Roughness Model 11 2.3 Historical Perspective and Literature Review 13 2.3.1 P/PI Detuning 13 2.3.2 Deterministic Constrained Optimization 16 2.3.3 Constrained Minimum Variance Controller 17 3 Linear Optimal Controllers 19 3.1 Introduction 19 3.2 State-Space Method 20 3.3 Evaluation of Variances 22 3.4 Stationary Input 23 iii 3.4.1 System Equation and Performance Index 23 3.4.2 State-Space Solution . 24 3.4.3 Noise-Free Observer 25 3.4.4 Controller Parameterization 26 3.4.5 Comparison to PI Controller 29 3.5 Minimum Variance of Outlet Flow 32 3.5.1 Comparison to P Controller 33 3.6 Random Walk Input 35 3.6.1 Effects of Damping Factor 37 3.7 Summary 38 3.8 Mathematical Details 39 3.8.1 Solution of Section 3.4.2 39 3.8.2 Rv(y) and Rv(u) in 3.4.4 39 3.8.3 Rv(y) and Rv(u) in 3.4.5 41 4 Downstream Processes 43 4.1 Introduction 43 4.2 Downstream Process Model 44 4.3 State-Space Solution 47 4.4 Performance Comparison 49 4.4.1 Effects of v(0 weight 50 4.4.2 Var[y2] reduction versus CJP 52 4.5 Summary 54 5 Markov Jump Systems 55 5.1 Introduction 55 5.2 Jump Markov Systems 56 5.3 Holding Time Density 57 5.4 State Distribution 59 5.4.1 Entry Probabilities 63 5.4.2 Examples 63 5.5 State Moments of Linear Jump Systems 71 5.5.1 Example 73 5.6 Linear Quadratic Optimal Control 75 5.7 Summary 77 5.8 Appendix - Mathematical Details 78 5.8.1 Covariance of x(f) 78 5.8.2 Linear Quadratic Control 80 i v 6 Broke Storage Tank Level Control 84 6.1 Introduction 84 6.1.1 Paper Machine 84 6.1.2 Benefits of Wet End Stabilization 85 6.1.3 Broke Handling 87 6.1.4 Literature Review 87 6.2 Characteristics of Broke Flow 88 6.2.1 Spectrum of Broke Flow ^ 89 6.3 Flow Smoothness Requirements' . 92 6.4 Linear Optimal Controller 92 6.4.1 Jump System Theory 94 6.4.2 Implementation Issues 95 6.5 Nonlinear Controllers . 95 6.5.1 Minimum Overflow Probability Controller 96 6.6 Optimal Change of u(i) 100 6.6.1 Asymptotic Analysis 101 6.6.2 Linear u(i) 103 6.7 Summary 103 6.8 Maximum Principle 104 7 Simulation With Industrial Data 110 7.1 Introduction HO 7.2 Data HO 7.2.1 Process HO 7.2.2 Data Set C 112 7.2.3 Incoming Flow Estimate 113 7.2.4 Statistical Analysis 115 7.2.5 Probability Plots 117 7.2.6 Chi-Square Test 118 7.2.7 Summary 119 7.3 Controller Design 120 7.3.1 QMOPC 120 7.3.2 Linear Controller 122 7.4 Simulation 124 7.4.1 Manual Control 125 7.4.2 Linear Controller 125 7.4.3 QMOPC -127 7.5 Summary 127 v 8 Conclusions 130 8.1 Summary 130 8.2 Future Research 132 Bibliography 134 vi List of Tables 7.1 Chi-square test of break and normal durations 119 7.2 Overflow probabilities for QMOPC and MOPC 120 7.3 Overflow probabilities for QMOPC 121 7.4 The controller parameter for <x[y] = 20 123 7.5 Statistics of simulations 124 vii List of Figures 1.1 Averaging level control problem 1 2.1 Storage tank material balance 8 3.1 Parameters of C L vs. K 27 3.2 Gain of C L for K = 0.1,1 and 10 28 3.3 Magnitude of S for C L and C P I 28 3.4 Comparison of C L and PI controller with n - 0.7 30 3.5 Comparison of C L and the PI controller with various n 31 3.6 The ratio of the variance of u(t) between CPD and Cp. x-axis is Rv(y) 34 3.7 Performance change by rj for PI controller 37 4.1 Downstream process and F u 44 4.2 Gain plot of the downstream process transfer function 46 4.3 System for tank level control loop and downstream process 47 4.4 Gains of Co and CL 51 4.5 Magnituedof5 forC D andC/, 51 4.6 Variance reduction of y>i(j) vs variance of v(0 52 4.7 Effects of o)p on variance of y2{t) 52 4.8 Effects of a>p on variance of ^ ( 0 with normalized downstream process cut-off frequency 53 5.1 Symmetrical binary signal d(t) is filtered by a low-pass filter. 64 5.2 Probability densities of filtered binary signal for B = 0.5,1, and 2 65 5.3 Probability density of filtered binary signal for fi = 10 66 viii 5.4 0i, 02 and <fo calculated with N = 300 69 5.5 Effects of state discretization 69 5.6 Comparison to simulation results 70 6.1 Simplified diagram of a paper machine wet end 85 6.2 Alternating Poisson process 88 6.3 Implementation of Ci 96 6.4 The minimum overflow probability strategy. 98 6.5 The phase plane of (u(t),y(t)) under MOPC 99 6.6 The logic of QMOPC : ; 99 6.7 Graphs of the optimum change strategy 102 6.8 Jc/JminvsaT 104 7.1 Simplified diagram of the broke management system I l l 7.2 Graphs of tank level y(t), outlet flow u(t), and break detector 112 7.3 Graphs of estimated incoming flow and level 115 7.4 Graphs of model level y m ( t ) 116 7.5 Histogram of the estimated incoming flow Fb{t) 116 7.6 Probability plots of break and normal duration 118 7.7 Sample distributions and theoretical distributions of break and normal durations . . 119 7.8 Probability distribution of y(t) under the QMOPC 122 7.9 Graphs of cr(u) and cr(u) vs damping factor 123 7.10 Probability distribution ofy(0 and u(t) under C L 124 7.11 Comparison of the QMOPC, Cz,, and manual control 125 7.12 Histograms ofy(r) (left) and u(t) (right) by manual control 126 7.13 Histograms of y{i) (left) and u(t) (right) by C L 126 7.14 Histograms of y(t) (left) and u(t) (right) by QMOPC 127 7.15 Comparison of u(t) by manual, C L and QMOPC 128 7.16 Comparison oiy{t) by manual, CL and QMOPC 129 ix .1.1 V, Acknowledgements During the research of this thesis, I have been fortunate to receive help from many people, and I would like to express my gratitude to them. First, to my supervisors, Prof. Guy Dumont, Prof. Michael Davies, and Dr. Bruce Allison. They warmly encouraged me to take up the challenge, and provided me with steady guidance throughout my study. Their patience with and critiques for my ever changing ideas, are very much appreciated. Dr. Allison introduced me to a fruitful area of averaging level control, and kick-started the research. It was one of highlights of my student life to exchange ideas with him on many subjects. Also thanks to Dr. Maryam Khanbaghi for introducing Markov jump systems and sharing her knowledge on broke tank level control with me. To the Network of Centres of Excellence, Mechanical Wood-Pulps, the Science Council of British Columbia, and the Pulp and Paper Institute of Canada, I thank for their financial assistance. I would like to thanks Messrs. Rick Harper and Leo Pelletier of NorskeCanada for supplying the industrial data. Without their help, this thesis could not have been completed. To my fellow graduate students, Dr. Junqiang Fan, Stevo Mijanovic, Manny Sidhu, and Dr. Michael Chong Ping, I would like to express my sincere gratitude for their friendship, and kindness to share knowledge with me. I appreciate that they treated me as their peer even though I am much older than them. I felt that I was young again. I thank Dr. Leonardo Kammer of PAPRICAN for his critique and advice on my work. The library services in the UBC and the Pulp and Paper Centre was essential to my research, and I thank the people their. Special thanks goes to Ms. Reta Penco for helping and teaching me in literature search. I am fortunate to have many friends who have supported me with their warm friendship. To my old MacBlo gangs, Glenn Swanlund, Eric Lofkrantz, and Ted Carlson, I thank for their continuing x friendship. It is comforting to know that my old friends in Japan are always ready to help me. I would like to thank Dr. Toshio Kojima and Toshinori Watanabe for their encouragement. I would like to thank two of my friends, Bruce and Dr. Alf Isaksson who have been inspiration for me to pursue the study of control engineering. Bruce rekindled my enthusiasm for control theory when he implemented digester level control at Port Alberni more than 20 years ago. Alf's enthusiasm and dedication to control engineering is contagious. I would like to thank my supervisor at the University of Tokyo, where I studied 35 years ago as a master student, Prof. Mitsuru Terao who instilled the spirit of control engineering which I have kept ever since. Prof. Karl Astrom, after his retirement, told us that he had the best times in his life during his graduate study and retirement. No wonder I enjoyed so much because I was a graduate student and a retiree at the same time! However, one drawback of being an old graduate student is that I cannot share the joy of completing the research with my parents, whose memory sustained me during tough times. xi Chapter 1 Introduction and Overview 1.1 Background and Motivations The research presented in this thesis started with a goal of devising a practical and effective auto-matic control scheme for broke storage tank level, which was identified as the most influential and difficult process for stabilizing a paper machine wet end. The importance of stabilizing a paper machine wet end and the role of the broke storage tank are detailed in Chapter 6. It was soon realized that the broke storage tank level control problem is a special case of the averaging level control problem depicted in Figure 1.1. Here, the incoming flow to the storage tank Fd(f) changes Fd(f) (disturbance) Downstream Process Figure 1.1: Averaging level control problem randomly and the tank level y(t) must be keptbetween the high limit yH and the low limit yi. The goal of averaging level control is to manipulate the outlet flow Fu(t) to minimize its adverse ef-fects on the downstream process(es), while at the same time keepingy(t) between^ a n d ^ . The yH yL m r Storage Tank I -A Fu(t) 1 1.2. Overview and Structure of the Thesis 2 broke storage tank level control problem has a very specific input disturbance, but the downstream processes are diffuse and difficult to pin down by a simple mathematical model. Before embarking on the task of designing a broke storage tank level controller, available tech-niques for averaging level control, stochastic and deterministic, were reviewed. Averaging level control can benefit many chemical processes that include surge tanks. Also, many areas in the pa-per machine wet end can be improved by proper applications of stochastic averaging level control. Chapters 2, 3 and 4 treat the stochastic averaging level control problem in a unified approach with optimization methods emphasizing the importance of the input disturbance characteristics and the downstream processes. The fact that the input disturbances to the broke storage tank can be modelled by a continuous time Markov process (Khanbaghi, 1998) motivated the study of Markov jump systems. Although it was found that jump linear quadratic theory as applied to the broke storage tank level control problem has no advantages over linear quadratic optimal control theory, a new theorem and nu-merical methods for evaluating the state variances and the state distributions are devised. The new methods are utilized in the synthesis and analysis of control schemes in Chapter 7. With the techniques of stochastic averaging level control and Markov jump system theory, linear and nonlinear controllers for the broke storage tank level control problem are designed and analyzed. The new control schemes are tested in simulations with industrial data. 1.2 Overview and Structure of the Thesis The thesis has three main parts: (1) stochastic averaging level control, (2) Markov jump systems, and (3) broke storage tank level control. Trie first two parts provide technical tools necessary to attack the problem treated in the third part. Averaging Level Control Three consecutive chapters 2, 3, and 4 cover the stochastic averaging control problem. In Chapter 2, the stochastic averaging level control problem is defined. The importance of the input distur-bance characteristics and downstream processes is highlighted. Existing techniques for determin-istic and stochastic averaging level control are reviewed. The popular random-walk assumption on 1.2. Overview and Structure of the Thesis 3 the input disturbances is scrutinized. A simple first-order low-pass disturbance model is preferred. Chapter 3 treats the problem with linear quadratic control methods. A state-space approach with noise-free observer is utilized. The small dimension of the problem makes it possible to obtain solutions in closed form. Thus, straightforward algorithms and good perspectives are obtained. In Chapter 4, the downstream process is included explicitly in the problem formulation. The performance of the resulting controller is compared with ones designed without downstream pro-cesses. A convenient parameterization is devised to streamline the performance comparison. Markov Jump Systems Chapter 5 presents new theorems on state probability distributions and their applications for the variance evaluation and the state distribution calculations. This chapter is independent from the previous chapters and can be read independently from the rest of the thesis. The new numerical method for calculating the state distribution is utilized in designing controllers for broke storage tank level. Broke Storage Tank Level Control Chapter 6 provides motivation and background for broke storage tank level control. Difficulties arise from the input disturbance characteristics, which are quite different from those encountered in usual averaging level control problems, and defining appropriate outlet flow smoothness re-quirements. A new nonlinear control scheme called the minimum overflow probability controller (MOPC) is introduced. Also, the optimal linear controller is designed for comparison. In Chapter 7, industrial data collected from a paper machine where the broke storage tank level is controlled manually, is used to estimate the incoming disturbance Fd(f) from a six-month history of the broke storage tank level y(t) and outlet flow Fu(t). The two controllers designed in Chapter 6 are tuned for Fd(t) and tested by simulation. Comparisons are made among the two controllers and manual control. 1.3. Contributions of the Thesis 4 1.2.1 Format and Notation As the three parts of this thesis deal with subjects more or less independent from each other, overviews and literature reviews are placed at the beginning of each part, not at the beginning of the thesis. Detailed mathematical derivations are included at the end of each chapter rather than at the end of the thesis. This is a compromise between two goals: that the thesis should be comprehensive and that mathematical details should not impede the smooth logical flow. Most of the mathematical symbols are standard in the literature. The same symbol is used with different meanings, depending on the context, if there is no danger of confusion. Economy in notation is obtained by using the same symbol with different arguments to represent a time signal in continuous time, discrete time, and in the Laplace transform. For instance, y(t) denotes a process variable in continuous time, and y(s) the Laplace transform of the same variable. For discrete time, the standard integer symbols k etc. are used, i.e. y(i) denotes y(t), t = iAt, At being a sample period. A notation ":=" is used to signify "is defined to equal". Blackboard bold letters: F[event] means the probability of event and E[X] means the expectation of a random variable X. All the controllers in the thesis are meant to manipulate the outlet flow Fu(f) to control the tank level. The formal controller output is denoted with u(t), which may be u(t) = Fu(t) or u{t) = Fu(t) - E[Fd(t)]. In either case, Var[Fu(r)] = Var[w(r)] and Fu(t) = ii(t). The rate of change in the outlet flow Fu{t) appears often, and is given a special symbol v(t): d „. d , „ This convention is used throughout the thesis, but the definition of v(r) reappears a few times as reminder. Sometimes, u{t) is used as the manipulated variable of an abstract problem, such as, x(t) = Ax(t) + Bu{t). The meaning of u(t) should be clear from the contexts. 1.3 Contributions of the Thesis Averaging Level Control • The stochastic averaging level control problem is put on a solid theoretical basis by iden-tifying the importance of input disturbance characteristics and the outlet flow smoothness 1.3. Contributions of the Thesis 5 definition. Exploiting the continuous time state-space approach, clear perspectives are ob-tained for synthesis and analysis of the linear optimal controllers. • The random walk input disturbance assumption is given a detailed analysis. Advantages of stationary disturbances are presented; The optimal linear controller for a first-order low-pass input disturbance, which is a phase lag network, is synthesized. • A new problem formulation to include the downstream process explicitly is presented. This helps to clarify the meaning of outlet flow smoothness. With a simplified downstream pro-cess model, practical controllers are obtained. • A proper state-space formulation is presented when the input disturbance is a random walk. The key is to make the state vector x(t) stationary, so that Var[x(r)] exists. Markov Jump System • A new theorem on the state probability distributions is given in Chapter 5, which leads to new methods for evaluating the mean and the variance of the state of a linear jump system, and a new reliable numerical method to calculate the state distributions. Unlike approaches based on simulations, the new method does not have time-discretization approximations nor randomness due to computer generated random numbers. These results are utilized to evaluate the overflow probabilities of controllers. • Linear quadratic optimal theory for jump systems, where the state has discontinuity, is ex-tended to include load-change disturbances in a natural way. This formulation provides a more flexible means to include Fu(t) in the performance index. Broke Storage Tank Level Control • The actual incoming flow to a broke storage tank Fb(t) was estimated from the tank level and outlet flow measurements. This is the first published study of its kind. It was found that break and normal durations follow the exponential distribution quite well. This is the second confirmation of the published results of Khanbaghi et al. (1997). However, Fb(t) is not truly binary as was assumed by the previous studies based on material balance. 1.3. Contributions of the Thesis 6 • A practical nonlinear controller with hard constraints on the outlet flow Fu(t) and its rate of change Fu(t) is presented. The controller minimizes the overflow probability and guarantees that the level will stay above the low limit yi. The new method for calculating the state dis-tribution of Markov jump systems makes the design process straightforward without relying on simulations. The nonlinear controller is compared to the optimum linear controller and manual control with the industrial data, and is found to control the tank level better and with a smaller maximum outlet flow. • The optimal change strategy of Fu(t), where Fu(t) is changed from its initial value UQ to the final value u\ in a given time T, is investigated. The performance index is the integrated square error of the downstream process variable. With a simple downstream process model, the problem is solved with the Maximum Principle's singular control. A closed form solution is obtained, which makes it possible to compare linear change strategy to the theoretical optimum. Chapter 2 Averaging Level Control 2.1 Introduction Averaging level control refers to liquid level control of storage tanks, where the objective is to keep the outlet flow as smooth as possible against the fluctuating inlet flow, while at the same time keeping the tank level within high and low limits (Figure 1.1). The term "averaging" indicates that the average of the tank level is significant rather than the level itself. For example, Shinskey (1967) (page 147) states, "This application is often referred to as 'averaging level control,' because it is desired that the manipulated flow follow the average level in the tank". The reason that the outlet flow Fu(t) should be smooth is to minimize its adverse effect on downstream processes. This is quite different from "normal" regulatory control problems, where the goal is to keep the controlled process variable y(t) close to its (constant) target yr. Even for regulatory problems, excessive movements in the manipulated variable u(t) are usually avoided for a number of reasons, including to avoid saturation in u(t) and to make a controller robust for modelling errors. However, smoothing of u(t) is never a principal goal of control. An underlying idea for investigating the averaging level control problem in the thesis is constrained optimization, where the "best" controller is a one that minimizes some performance index J subject to a constraint condition g(x, u) < 0. This idea provides systematic procedures to design controllers for given situations, which are characterized by J and g. Also this methodology makes it possible to compare many published control schemes in a consistent and unified way. 7 2.2. Mathematical Models 8 2.2 Mathematical Models Optimization problems are tackled either analytically or numerically. In either case, mathemat-ical models for the problem are necessary. Mathematical models for practical problems are al-ways a compromise between fidelity to the target objects and mathematical tractability. To judge the fidelity and usefulness of the mathematical models, intimate knowledge of the subject area is necessary. Thus, construction of practically useful mathematical models is one of the important activities that distinguish control engineering from applied mathematics. The three main models needed for averaging level control are: 1. disturbance model 2. process model 3. flow smoothness or roughness model In the following sections, the mathematical models for averaging level control are discussed in detail. This is important because historically mathematical models for the input disturbance and the flow roughness have been treated rather indifferently, although it will be found that they influence controller design greatly. 2.2.1 Process Model Figure 2.1 shows a simplified diagram of a material balance around a storage tank, where V(t) denotes the liquid volume in the tank, sn.dy(t) the tank level. Fj(t) denotes the incoming flow, and Fu(t) the outgoing flow, both expressed in volume!time. Assuming that the liquid is noncompres-F u ( t ) Figure 2.1: Storage tank material balance. V(t) is the volume of liquid in the tank, and y(f) is the level. 2.2. Mathematical Models 9 sive, we obtain the following equation on V(t): V(f)=V(0)+ f[Fd(t')-Fu(t')W (2.1) o Dividing the both sides of (2.1) by the cross-sectional area of the tank, A, yields, y(t) = v(0) + Kp f[Fd(t') - Fu(t')]dt', (2.2) where K p = 1 /A is the process gain. This simple integrator model works well in practice, and is used throughout the thesis. When the cross-sectional area varies by height, K p is not a constant, and (2.2) becomes nonlinear. However, linearity can be kept by treating the problem by (2.1). The level high limit is replaced by the volume high limit. A storage tank for noncompressive material is a rare example in process control of a simple mathematical model derived from the first principles that works well. In a physical process, many dynamical subsystems, including manipulating valves, and measurement sensors with noise and delay, are involved. However, they can be ignored, at least for liquid storage tanks in a paper machine wet end because the execution time of level control can be slower than other loops. For instance, a broke storage tank level can be controlled using control intervals of 30 seconds to 5 minutes, while most pressure or flow loops are executed at 1 second intervals. Therefore, these dynamics can be ignored by the level controller that is executing every 30 seconds to a few minutes. This fact is generalized in Marlin (1995), page The process has no dead time and a phase lag of only 90°, indicating that feedback control would be straightforward for tight level control. This is actually the case in many systems, since the sensor and valve dynamics are usually negligible. The system parameter K p can be obtained with a good accuracy from experiments or physical dimensions of the tank. Thus, robustness problem due to the process uncertainty is not a concern. 2.2.2 Disturbance Model The role of the disturbance model was rather neglected in the early optimum control literature. However, its importance, especially in process control, is now well accepted. Still this is not always the case as stated in Harris and MacGregor (1987): 586: 2.2. Mathematical Models 10 All optimal controllers (i.e., those optimizing a performance index) must be based on a specified disturbance (or set-point change) model. Often in the design of such controllers the disturbance model is not explicitly stated, but rather indirectly implied in the problem formulation. In averaging level control, very few papers have paid detailed attention to the input disturbance. The term "averaging" itself suggests some stationary stochastic process because non-stationary process cannot be averaged. However, many papers accept the need to bring back the tank level to a setpoint to maximize the surge capacity. This suggests that the incoming flow is non-stationary since if Fdif) is stationary there is no need to bring back the tank level to its setpoint as stated in Shinskey (1988). Probably the authors of these papers envisioned the input disturbance as a combination of stationary rapid fluctuations, which can be averaged out, and less frequent step-type load changes. This kind of disturbance poses interesting control problems that involve statistical detection of abrupt changes or hidden Markov processes. However, this is left as one of the future endeavours. Most of the performance analyses are based on responses to step disturbances, probably be-cause of mathematical tractability. If the input disturbance is step-like and stays at a constant value long enough for the closed loop system to settle down to a quiescent state, then the analyses based on the step responses would be valid. However many real disturbances do not have this characteristic and care must be taken in interpreting analyses based on the step assumption. In the paper machine wet end, disturbances for storage tanks are mostly stochastic. Most of the studies in a stochastic context assume that the input disturbance is a random walk, which is non-stationary. The random walk is a popular choice for process control because it introduces integrating action in the controller. However, in reality, all disturbances are bounded and most of them can be deemed stationary if observed for a long enough time. In this thesis, stationary disturbance models are preferred because there is no strong need to eliminate offset for averaging level control. A simple stationary process with first-order dynamics is used as a disturbance model. Let d(t) denote the input disturbance deviation variable d(t) := Fj(t)-E[Fd(f)] and w(t) white noise. It is assumed that d(t) is the output of a first-order low-pass filter with the cutoff frequency d(s) = -^—w(s). (2.3) S + 0)d 2.2. Mathematical Models 11 This model was used because most disturbances have some kind of roll-off frequency. Also the input disturbance model for broke tanks, which is a binary-valued random process, can be charac-terized with respect to the frequency power spectrum with this model. 2.2.3 Flow Smoothness/Roughness Model All the literature in averaging level control states that one of the control objectives is to keep the outlet flow smooth. However, the mathematical definition of flow smoothness varies from author to author. Furthermore, some papers contain no mathematical definition of flow smoothness at all, so it must be guessed from context. To apply optimization based methods, which minimizes some performance index it is more convenient to express flow "roughness" or "variability", denoted with r(u(t)). In deterministic cases, the controller tries to minimize r(u(t)) subject to the level constraints. In stochastic cases, a natural choice for the flow roughness is the mean value of r(u(t)),E[r(u(t))]. Since the purpose of outlet flow smoothing is to minimize its adverse effect on the downstream processes, the flow roughness indicator should be selected taking into account the downstream process. As shown in Section 4.4, for downstream processes that can be approximated by a first-order system, Var[«(r)] = Var[v(r)] works, well if the downstream process time constant is much shorter than that of the closed loop tank level system. If the downstream process mixes various material streams, and it is desired to keep mixing ratios constant, then Var[w(r)] rather than Var[v(/)] might be the better choice. However, if the input disturbance is non-stationary, Var[w(r)] does not exist, and cannot be used. For broke tank level control, Bonhivers et al. (1999a) and Dabros et al. (2003) expand the idea of downstream process to a certain paper quality index. Crisafulli and Peirce (1999) contains an excellent concrete example of the downstream process and industrial data. The process is a raw cane sugar factory, where crushed cane fibre and water, called "mixed juice", is heated and fed to a clarifier to extract clear juice. This clear juice is subsequently crystallized to produce the final product, sugar. The mixed juice goes through two surge tanks connected in series before it is fed to the clarifier. The paper reports a new feedforward control scheme to maximize the combined surge capacity of the two tanks resulting in smooth inlet flow and reduced turbidity of the clear juice. The paper says that, "the main technical control 2.2. Mathematical Models 12 objective was to minimize the rate of change of the secondary mixed juice flow." Since the inlet flow to the surge tank fluctuates randomly, E[|v(f)|], or the more mathematically tractable E[v(/)2] = Var[v(r)], would be an appropriate mathematical measure for flow roughness in this case. MRCO The maximum rate of change in the outlet flow (MRCO) is mathematically defined as MRCO := max |v(0|. 0<Koo Although the definition implies an infinite horizon, the MRCO makes sense in a finite horizon context because, for a step input disturbance, the maximum occurs early in the step response. The MRCO is used as the flow roughness indicator in many papers since it was used for PI tuning analysis with step responses (Cheung and Luyben, 19796). For step disturbances, the MRCO is mathematically tractable for P and PI controllers to the extent that it is expressed in closed form. Also intuitively, the MRCO seems to represent some kind of flow roughness. However if the MRCO is taken at its face value (not as a control parameter for step responses by the PI controller), from the point of view of the downstream processes, its effectiveness as a flow roughness indicator is limited. For most chemical processes, the effect of the outlet flow is cumulative, not a function of the value at a single time point. There may be some processes that the MRCO must be kept below some limit value to prevent catastrophes from happening. Then the problem should be posed as the minimization of some other flow roughness indicator subject to the MRCO limitation, not the minimization of the MRCO itself. Another drawback of the MRCO is that it is not cumulative, thus many analytical and numer-•'XX-. > ical optimization techniques, such as dynamic programming or quadratic programming, are not applicable. This non-cumulative nature of the MRCO may cause impractical (very oscillatory) solutions by linear programming as reported in Sidhu (2003). Therefore, the MRCO is not used in this thesis. 2.3. Historical Perspective and Literature Review 13 2.3 Historical Perspective and Literature Review It was recognized from the early stage of process control technology development that surge tanks require different control objectives from normal regulation loops, and many papers and sections of books have been written for averaging leve^cpntrol since the 1960s. They may be categorized into the following three groups: (1) P/PI Detuning!;(2) Deterministic Constrained Optimization and (3) Stochastic Constrained Optimization. Each approach is discussed in detail in the sections below. 2.3.1 P/PI Detuning Studies published up to the early 1980's constitute this group. The ideas expressed in these papers reflect control hardware prevalent at that time, pneumatic or electronic single loop controllers. They are summarized as follows: 1. no mathematical definitions of the input disturbance nor outlet flow smoothness are given. 2. performance analyses are mostly non-mathematical being heavily depending on visual in-spection of responses for a step disturbance. 3. controller synthesis is mostly ad hoc or heuristic. This lack of clear definitions both of the input, disturbance and of the meaning of flow smoothness lead to a number of different schemes claiming good performance. Proportional Only Control This is advocated by Shinskey in Shinskey (1988) and Shinskey (1994), although originally he preferred nonlinear controllers (Shinskey, 1967). Since there is no integrating action, the tank level tends to be away from its setpoint (off-set). To some authors, this off-set is a shortcoming of the P controller. However, when the input disturbance is stationary or bounded, the off-set is not disadvantageous. This point is clearly stated in (Shinskey, 1988): "Not only is there no need to return the level to some set point such as 50 percent, but this practice actually reduces the effective capacity of the vessel." Although, the meaning of flow smoothness is not stated clearly, it appears to be defined as a smaller variation of u{f). It is shown later that the P controller is close to optimal when the input disturbance is a first-order low-pass stationary process and the flow roughness is represented by Var[w(r)]. Thus, if me input disturbance is such that the downstream process 2.3. Historical Perspective and Literature Review '.;)';[ 14 operates better with a smaller Var[w(/)], then the P controller performs well. However, care must be taken since for some downstream processes, Var[v(0] is abetter indicator of the flow roughness. Then a phase-lag network or low-pass filter is a better choice. PI Tuning Rules In Cheung and Luyben (19796), tuning rules for PI controllers are given based on analyses of step responses. The design parameters are the maximum peak height (MPH) and the MRCO. No justifications for using the MRCO as a flow roughness indicator was given. It seems that the MRCO was selected for its mathematical tractability in step response analysis. The tuning procedure in the paper determines the PI controller settings from two design specifications, the MPH and the MRCO. This procedure potentially leads to a very lightly damped system depending on how the MPH and the MRCO are selected. A better procedure would be to fix the damping factor n to some reasonable value between 1 and 0.7, specify the desired MPH, and accept the resulting MRCO as a price to pay for obtaining the target MPH. This is done in Marlin (1995) with n = 1. Nonlinear PI Controllers In Shunta and Fehervari (1976), two kinds of nonlinear PI controllers were proposed, one with two settings switched by the level error (when the level error is small, slower tuning makes the outlet flow change more smoothly, and faster tuning kicks in when the level error exceeds a limit to prevent level limit violation), and the other called "wide range" which changes the proportional gain and reset time continuously. No tuning guidelines were given in the paper. To estimate the behaviour of a nonlinear controller against random input disturbance is a difficult problem even if the input disturbance distribution is available. Tuning of these controllers must be done on a trial-and-error basis. It seems these controllers did not gain much popularity. The basic weakness of nonlinear schemes based on the level error is that a large step-like input disturbance is not detected until the level error exceeds the threshold. It is best to act immediately against a large input disturbance to minimize any flow roughness indicator based on v(t) - u(t). As discussed later, the model predictive control (MPC) methods address this problem properly. 2.3. Historical Perspective and Literature Review 15 PL Level Controller A Proportional-Lag (PL) level controller was proposed in Luyben and Buckley (1977) and further studied in Cheung and Luyben (1979a). It is claimed that their approach maintains most of the desirable features of proportional-only control but eliminates off-set. It is shown below that the PL control is a PI controller when measurement noise or delay is negligible. In the PL controller, the outlet flow u(f) is a sum of a proportional term of the level error and a filtered inlet flow w(t): u(i) = KLe(t) + w(t), (2.4) where K L is the proportional gain, e(t) is the level error (y(f) -y r). The filter is a first-order low-pass filter KF/(TFs + 1) (Cheung and Luyben, 1979a). In the following analysis, KF is set to 1 to avoid an off-set in the level. Also there is no advantages in setting KF ± 1. Without losing generality, it is assumed that the level setpoint.yr is zero, so e(t) = y(t). Using the same symbols for Laplace transforms of time signals, we get u(s) = KLy(s) + -^-. (2.5) TfS + 1 The process is y(s) = ^[w(s) - «(*)]. (2.6) s Then, the transfer function from w to y, GpL, is y ( s ) _ KpTfS , v^(sj " (TFs + \)(KLKP + s) K~ls2 + s(KL + (KpTF)-1) + KLT Now, we close the loop with a standard PI controller Kc(l + 1 / T r s ) , where Kc is the proportional gain and Tr the reset time. Then, the transfer function from w to y, Gpj, is Q .= y W = £ (2S) P I' w(s) K~ ls 2 + K cs + KcT~l' K ) From (2.7) and (2.8), for a given KL and TF, we can make GPL = GPI by selecting Kc and Tr as follows: Kc = KL + -r\=-, (2.9) Kpl F and Tr = TF^ = TF + -^— (2.10) Thus, the PL controller has no special advantages over the PI controller. The performance may be different from the PI controller if dead-time or large measurement noise is involved. 2.3. Historical Perspective and Literature Review^'V 16 2.3.2 Deterministic Constrained Optimization Supervisory computers and computer based distributed control systems (DCS) became widely available in the 1980s, and computation intensive control algorithms such as model predictive control (MPC) became practical as a result. Papers reviewed in this section are characterized by their handling of level constraints in a deterministic context (step input disturbance). DMC In Cutler (1982), the Dynamic Matrix Control (DMC) algorithm was applied for averaging level control. The paper does not give mathematical details, but the comments here can be inferred. In normal DMC, the future trajectory of the controlled variable in discrete time, y([k, k + Np]) (a positive integer is a prediction horizon) is compared to the desired trajectory yr([k, k + Np]) and the squared sum of the error E^'CKO ~>v(0)2 ^ a function of control increment, Au(k), is used in the performance index to be minimized. For averaging level control, the desired trajectory is set as a funnel shaped region centred around the level setpoint. The tank level error is deemed to be zero if it is within the funnel shaped region. Thus, the controller makes no move (Aw(r) = 0) if the level error is inside the funnelled region, or a smaller move just enough to bring back the level into the region if the level error is outside. The paper does not mention level constraint handling, which is one of the strength of DMC. It reported good filtering characteristics of the algorithm when the input disturbance is small so that the predicted level stays inside the funnel shaped region. Although a number of simulation results were included, actual applications of this DMC algorithm to averaging level control problems have not been reported (McDonald and McAvoy, 1986). Optimal Predictive Controller In McDonald and McAvoy (1986), the averaging level control problem was studied in continuous time for a step input disturbance of known amplitude. The control objective is to minimize the MRCO while keeping the level within constraints. Due to the simplicity of the plant dynamics and the MRCO criteria, analytical solutions, which are nonlinear, were obtained. The amplitude of the input disturbance is estimated from the level and the outlet flow, and the implemented control algorithm takes a form of predictive control, and thus, was named the optimal predictive controller 2.3. Historical Perspective and Literature Review 17 (OPC). The OPC does not have the weakness of the level error based nonlinear control schemes. The authors thought that it was necessary to bring the level to its setpoint, and combined a PI algorithm with the nonlinear solutions. The matter of how quickly the level should be brought back to the setpoint can only be discussed with stochastic input disturbances. The MRCO is difficult to use in a stochastic scenario. However, the paper treated the input disturbance as steps of undefined length. This makes the tuning procedure ad hoc, heavily depending on simulations. The contribution of the paper is to put the averaging level control problem in the MPC framework. MPC The problem discussed (above) by McDonald and McAvoy (1986) was treated in discrete time in Campo and Morari (1989). As in continuous time, if u(t) or Au(t) is not constrained, analytical solutions are available. The need to bring the level back to its setpoint was acknowledged without discussion. The problem was posed as a finite horizon constrained optimization. To bring the level to its setpoint, a terminal condition that the predicted level at the end of the control horizon, y(te), is equal to the setpoint yr was included. This makes the control algorithm cleaner than adding integral action in ad hoc way. The horizon length becomes a tuning parameter. The algorithm is implemented as a receding horizon controller. With the finite control horizon, the problem can be solved by linear programming when \u(t)\ and |Aw(r)| are constrained. As was discussed for the case of OPC, the proper selection of the control horizon requires stochastic arguments, but the paper treats the problem in a purely deterministic context. Neverthe-less, the paper provides the most complete treatment of the problem for step-like input disturbances and the MRCO criteria. When \u(t)\ and |Aw(f)| are not constrained, the control algorithm is imple-mentable on any DCS as was reported in Allison and Khanbaghi (2001). 2.3.3 Constrained Minimum Variance Controller This methodology was introduced by Kelly (1998) for tuning PI controllers, and extended to more general cases including dead time by Foley et al. (2000). Both papers used a discrete-time for-mulation to investigate averaging level control as a stochastic optimal control problem. The input disturbance is assumed to be a random walk. The controller is synthesized to minimize the outlet 2.3. Historical Perspective and Literature Review 18 flow roughness (usually expressed as the variance of the outlet flow rate change Var[v(f)] in con-tinuous time or Var[Aw(&)] in discrete time) keeping the variance of the tank level Var[y(0] to a target value. The performance index J in discrete time, is set up as follows: J = Var[y(k)] + pVav[Au(k)] In discrete time, p can be set to zero to yield the minimum variance controller, where Au(t) is not constrained. The term "constrained" stems from the fact that p > 0 constrains Var[Aw(r)], and does not imply that the level is constrained. Constrained optimization of stochastic systems is a very difficult problem, and no practical algorithm for present day control hardware is available (Batina et al, 2002). The level is indirectly constrained by specifying its variance. Then, there is some probability that the level will violate its limits. For most practical applications, this method pro-vides satisfactory result because input disturbances in the real world are bounded. The gaussian assumption, which makes the input disturbance unbounded, is a fiction to expedite the mathemati-cal derivations. However, if the input disturbance happens to be larger than the design specification, the level will exceed the limit because there is no hard constraint mechanism. This aspect is differ-ent from OPC and MPC of the previous section, where closed loop hard level constraints are built into the algorithm. When there is no deadtime in the process, the optimal controller is a PI controller that makes the closed loop system second order with the damping factor 77 = V2/2. This fact was made clear by treating the problem in continuous time (Ogawa et al, 2002). This thesis extends this methodology to include stationary input disturbance and to include the downstream process in the performance index. Chapter 3 Linear Optimal Controllers 3.1 Introduction In this chapter, linear optimal controllers for averaging level control are studied. Although linear controllers cannot guarantee that the level constraints are met for random walk disturbances, they are considered for the following reasons: 1. For some problems, it is sufficient to detune existing linear controllers (PI or PID) in order to make the outlet flow smooth. The study here provides tuning guidelines. 2. The performance of linear controllers can be calculated analytically (and for some simple problems, in closed form). Usually this is not the case for nonlinear controllers, where simulation is the only reliable way to assess performance. The linear optimum controllers provide a reliable benchmark, which any non-linear controller must surpass. The input disturbance Fd{t) is assumed to be either a non-stationary or stationary stochastic pro-cess. The performance index is chosen to express the smoothness of the outlet flow Fu{t) in math-ematically tractable form, usually the variance of Fu(f) or a weighted sum of the variances ofFu(t) and Fu(t). Note that the latter choice is only available for stationary input disturbances. The output of a controller is denoted by u(t) and either u(f) = Fu(t), or u(i) = Fu(f) - Fm (Fm = E[Fd(t)]). In either case, Fu(i) = ii(t) = v(t), and Var[Fu(t)] = Var[w(0] and Var[Fu(r)] = Var[v(0]. Thus in the sequel, u{t) is used for the variance expression. Letj>(r) denote the tank level andyr its setpoint. A 19 3.2. State-Space Method 20 linear controller is sought that minimizes the following performance index J, J := E[(y(0 - yrf] + p2E[v(t)2 + n(u(t) - E[u(t)])2] = Var(y(0) + p2[Var(v(f)) + fiVar(u(t)), (3.1) where // = 0 when the disturbance is non-stationary. The weight p2 is a Lagrange multiplier (Newton et a l , 1957) so that Var[y(r)] is minimized subject to a constraint where a constant c is the upper limit of the flow variability. In practice, c is not specified, but p2 is adjusted to obtain an acceptable Var[y(r)]. Then the resulting optimum controller provides the minimum flow variability (or maximum smoothness). This approach was pioneered by Kelly (1998) and Foley et al. (2000) using discrete-time polynomial methods. In the thesis, the prob-lems are treated in continuous time using state-space methods, so that some characteristics of the resulting controllers, such as the closed loop damping factor, become evident. One advantage of the discrete-time approach is its ease of handling dead-time. However, the target applications are such that dead-time can be safely ignored. The problem is solved using the state-space linear quadratic (LQ) optimization methods combined with noise-free observers. Since most of problems in the thesis are single-input-single-output, the transfer function approach (Wiener-Hopf method) is a viable method, since there is no need for observer design. However, the state-space approach was adopted for the following reasons: 1. The algebraic Ricatti equations (ARE) can be solved analytically for this size of problem, and controller parameters, therefore,, can be expressed in closed form. The procedure is more straightforward than the factorization which is necessary for the Wiener-Hopf method. When the problem dimension increases, stable computer algorithms and software for solving the ARE numerically are widely available. Presently, this is not the case for the Wiener-Hopf approach. Var[v(0] + //Var[w(r)] < c, (3.2) 3.2 State-Space Method 3.2. State-Space Method 21 2. From the nature of the problem, noise-free observers are practical. Thus, the requirement of an observer is not really a disadvantage for the state-space approach. Rather it may, in some cases, add flexibility in controller design by providing a framework for selecting proper (or ad hoc) observers when the measurement noise cannot be ignored (Allison and Khanbaghi, 2001). Of course, by adding a noise model, the transfer function approach can produce the optimum controller that includes the observer. However, then the problem size becomes too large to accept easy factorization by hand and the approach loses the edge. The basic equations and notation which willibe used throughout this chapter are introduced here. The plant is expressed in the following state-space form. where x(t) is the state vector, u(t) is the manipulated variable (not Fu(t) here), w(t) is a Wiener process, and w(f) white noise. A, B, D are constant matrices of appropriate size. The white noise is expressed as a formal derivative of the Wiener process w(t). The derivative of the Wiener process does not exist in a strict mathematical sense, and the ordinary differential equation (ODE) (3.3) must be written as a stochastic differential equation (SDE). However, the SDE does not provide any advantage over the (formal) ODE in the development below. Thus the SDE notation is not used here. Generally, w(f) can be a vector process, but in the thesis w(t) and w(f) are always scalar processes. The performance index / is where Q is a positive semi-definite matrix and R a positive definite matrix. It is well known that in stochastic LQ problems where the white noise enters the system equation additively, the linear state feedback solutions are identical to those for the deterministic LQ problem (for example see Fleming and Rishel (1975)). The optimum control u(t)* that minimizes J is a feedback form, —x(i) = Ax(t) + Bu(t) + Dw(t), (3.3) (3.4) u(ty = -Kx(t). (3.5) The feedback gain K is related to a positive-definite solution P of the following algebraic Ricatti equation (ARE): PA + ATP - PBR-lBTP +Q = 0, (3.6) 3.3. Evaluation of Variances 22 and , K = RrlBTP. (3.7) A number of stable computer algorithms and programs are available for numerically solving the ARE. However, the low dimensionality of the problems here makes it possible to obtain closed form solutions easily. Thus, most of the AREs are solved analytically, so that the equivalence of the state-space approach and the Wiener-Hopf method can been clearly seen. 3.3 Evaluation of Variances It is necessary to calculate the variances of process variables to evaluate and compare the con-trollers. When a variable is the output of a rational transfer function driven by a signal with a rational spectrum, its variance can be calculated analytically. Ify(0 has a rational power spectrum density ._ N(s)N(-s) where Gv(-sr) := v 7 v , (3.8) N(s) = b„-is" 1 + • • • + bis + b0 D(s) = a„s"+an-is"~l + • • • + a\s + ao, its variance Var(y(r)) is calculated as Var[y(01 = ^ - £ S H ^ . (3-9) N(s)N(-s) _JOO D(s)D(-s)1 Let /„ denote the above integral when D(s) is an «-th order polynomial. Newton et al. (1957) contains a table that lists /„ up to n = 10. The first three /„ are listed below. h 2 h = 2a0ai / 2 = % i l % : ( 3 . 1 0 ) 2aoaia2 bja0ai + (b\ - 2bob2)a0a3 + b\a2a^ 2a0ai(a\a2 - a0a3) 3.4. Stationary Input 23 It becomes unwieldy for n > 3 for hand calculations. Astrom (1970) includes a recursive algorithm and a FORTRAN program for /„. In this chapter the following short-hand notation is used. IA — ) - = ^ j J j 7^ 7777^ 7 ~ds- (3.H) D{s)j 2nj J _ 7 0 0 D(s)D(-s) 3.4 Stationary Input In this section, the input disturbance Fd{i) is assumed to be the output of a first-order low-pass filter driven by white noise plus a bias. The state-space formulation is straightforward when the input disturbance is zero-mean. Let F m denote the mean of Fd(i), F m := E[Fd(t)]. The deviation variable d(i) for the input disturbance is defined. d ( t ) : = F d ( t ) - F m . (3.12) Then, d(s) = -^—w(s), (3.13) S + lOd where o)d is the cut-off frequency and w(s) the Laplace transform of white noise. 3.4.1 System Equation and Performance Index Since the process is an integrator, the mean of the outlet flow must be equal to Fm to keep the level stationary. Thus, u(t) is set as a deviation from its mean Fm. x\(t) = y(f) - y r tank level error xiit) - d(t) input disturbance u(t) = Fu(i) - F m outgoing flow deviation (manipulated variable) As shown below, x2(i) and u(t) appear as the input disturbance and outlet flow for the system equation. Xi(t) = Kp[Fd(t) - Fu(t)] = Kp[x2(t) + F m - (u(t) + Fm)] = Kp[x2(t) - u(t)]. (3.14) Since u(t) is stationary, not only the variance of u(t) = v(t) but also u(i) can be included in the performance index. / := E[(y(r) - yr)2] + p2E[v(t)2 + ^ u{t)2} (3.15) = Var[y(0] + p 2 (Var[v(0] + pVar[M(0]), 3.4. Stationary Input 24 where p is a weight for u(t) variation. The control weight is set to p 2 rather than p to simplify the calculations for the optimal gain below. When p > 0, the resulting controller tries to minimize the variances of both v(r) and u(t). Those cases are studied later in Chapter 6. Here, cases where p = 0 are investigated to compare with PI controllers for random-walk disturbances. To include u{f) = v(f) in the performance index in a straightforward way, v(r) is treated as a formal control input, and u(f) as an element of the state vector x(t), x(t) := x2{t) u(t)]T. Then the performance index J is expressed with the following two matrices Q and R, J = E{x(tf 1 0 0 0 0 0 0 0 up2 x(t) + v(t) ~R (3.16) (3-17) As X2(f) is the output of a low-pass filter with the cut-off frequency o>d, x2(s) = S + CJd Then the following system equation is obtained w{t) => X2(f) = -(OdX2(t) + 0Jdw(t). d_ Jt Xl(t) x2(t) = u(t) 0 Kp - K p 0 -o>d 0 0 0 0 3.4.2 State-Space Solution Analytical Solution The solution of (3.6), P, is expressed with its elements as P = P\\ Pn Pn ph P22 Pn Pl\ P32 P33 Xi(t) 0 0 x2(t) + 0 v(0 + u(t) 1 I T 0 (3.18) (3.19) (3.20) 3.4. Stationary Input 25 Substitution of A,B,Q, and R of (3.19) and (3.17) into (3.6) yields six second order polynomial equations on pij, from which a positive definite P can be obtained in closed form (algebraic details are included in Section 3.8.1). Since K = BTP/p2 = [pu p 2 i p 3 i ] / p 2 , (3.21) pn, pn and 7733, which are listed below, are required for K. P\3 = ~P, P23 = ; — — ^ — P i 3 = p ^ 2 p K p (3.22) Thus, K = [k\ k2 £3] is obtained as follows. k\ = —, k2 = , , k3 = J2KP p (3.23) P coj + ojd^KpT^ + Kp/p' V 3.4.3 Noise-Free Observer In most applications, x2(t) is not measured, and an observer is necessary. A simple noise-free observer is constructed by differentiating x\(f). Since = Kp[x2(f)-u(t)] => sx\(s) = Kp[x2(s)-u(s)], x2(s) = sKpXx (s) + u(s). (3.24) Usually, constructing an observer by differentiating the measured variable is not a good practice due to measurement noise. However, the problems here are formulated with a formal manipulated variable ii(t) and differentiating once does not introduce real differentiation when the real manip-ulated variable is u(i). Even if real differentiation is necessary, usually level measurements have low noise. Also, most applications can be executed with long sampling periods, say 5 seconds to a few minutes, which make it possible to filter the level measurement adequately if noise is not neg-ligible. Thus a noise-free observer is practical. See Allison and Khanbaghi (2001) for a Kalman filtering approach for noisy cases. The equation u(t) = v(t) = -Kx(t) is expressed in its Laplace transform, su(s) = -kiXi(s) - k2x2(s) - kiu(s). (3.25) Substitution of (3.24) into (3.25) yields, u(s) (s + k2 + k3) = xi (s) (-Jfci - sk2K-1). (3.26) •;5 n n 3.4. Stationary Input 26 Let CL(s) denote the transfer function of the controller, u(s) = CL(sMs)-yr) = CL(s)Xl(s). (3.27) Then where -k\ - sk2K~] s + k\KJk2 s + b CL(s) = = -k2K-1 = Kc , (3.28) w 5 + ^ + ^ 3 p s + k2 + k3 s + a v ' Kc — —k2/Kp a = k 2 + h (3.29) b - k\Kplk2. The sensitivity function S is S = 1 = ^ + f3 Kp s + b s2 + s(a + KPKC) + KpKcb' K' } 1 + — K - c s s + a From (3.29) and (3.23), a + KPKC = k2 + ki+ Kp(-k2/Kp) = h = ^2Kp/p, and KpKcb = Kp{-k2lKp)kxKplk2 = -^Kp = Kp/p. Thus, s _ s(s + a) _ s(s + a) s 2 + s y]2Kplp + Kp/p s 2 + y /2co c s + co 2 c ' where wc = ^Kp/p. The closed loop system is a second order system with damping factor rj = V2/2. The second-order maximally flat (Butterworth) filter has the same damping factor. This fact does not come out clearly in the discrete-time approach (Foley et al, 2000), where the sampling period introduces additional complexity. 3.4.4 Controller Parameterization For controller design, the following parameters are given: K p process gain o)d input disturbance cut-off frequency 3.4. Stationary Input 27 As a tuning parameter, the variance ratio Rv(y) between the input disturbance Fd(t) and the level y ( t ) is selected. Rv(y) := Var[y(0] VartxKO] (3.32) Var[Fd(t)] Var[d(t)] R v ( y ) is specified to meet the level constraint requirements (the standard deviation of the tank level is ^Rv(y) times the standard deviation of the incoming flow rate). A convenient way of characterizing the controller is by the bandwidth of the closed loop system relative to o>d, OJC K := (3.33) Substitution of y/Kp/p = OJC = KCJd into (3:23.) yields k2 and k-$ in terms of K and CJd- Then a and b are expressed with K and as follows. K2 + V2/c a - K2 + V2/c + 1 (jjd < did, , K 2 + V2/C + 1 b = cjd > cod-V2K+ 1 This shows that a < b , thus, the controller C L = K c ( s + b ) / ( s + a) is a phase lag network l . Figure 3.1 shows a and b for K between 0.01 and 100 for o j d = \. Figure 3.2 shows the gain of C L for Figure 3.1: Plots of a and b of C L = K c ( s + b)/(s + a) for K between 0.01 and 100. K = 0.1,1, and 10 for Kp = 1. Figure 3.3 shows the magnitude of the sensitivity function S for 'This is why the controller is denoted with Ci, where L stands for "lag". 3.4. Stationary Input 28 20 -10 -CD 10~2 10"1 10° 101 102 to Figure 3.2: Gains of CL for K = 0.1,1 and 10. (Kp = 1 and cod = 1) K = 0.1,1 and 10. In addition, for comparison purposes, the magnitude of S for a PI controller, of which shape is a function of the damping factor n only and does not change by K (see Section 3.4.5), is plotted. The plot marked as PI is for n = 0.7 and K= I. 10~2 10~1 10° 101 102 (0 Figure 3.3: Magnitude of the sensitivity function S for CL for K = 0.1,1 and 10. \S \ for CPI with n = 0.7 is also plotted. (Kp = 1 and o>d = 1) 3.4. Stationary Input 29 The controller performance is represented with the variance ratio between d(i) and ii(t), Rv(u), Var[w(0] Rv(u) := (3.34) Var[d(r)] Thus, for a given Rv(y), smaller Rv(u) indicates a better controller. With the use of the formulae in (3.10), Rv(y) and Rv(u) can be expressed in terms of K (algebraic details are in Section 3.8.2), K 2 K 4 + 3 V 2 V + 9K 2 + 6V2K + 3 Rviy) = — or, V2K(K 2 + V2K+ l ) 3 Rv(ii) = or y/2~K(K 2+ V2/f+ : il). ( K2(\ + V2K) (K 2 + V2K) + KA (3.35) [\K2+ yj2K+\) It is observed that Rv(y) ~ K2/iod and Rv(u) ~ o?d. The variance ratio of the level is proportional to the square of the process gain (this is easy to understand), and is inversely proportional to the square of the input disturbance cutoff frequency. When cod is high, more energy in the input disturbance is contained in higher frequencies and it is easier to filter out them with the same closed loop system, which is a low-pass filter. When o>d is high, the input disturbance is "busy", and increases Rv(u). To subsume the two parameters Kp and o)d in the design parameter, the "standardized" variance ratio Rv(y) and Rv(u) are defined as follows. Rv(y) := Rv(u) :•-Ry(yWd K 4 + 3 V 2 V + 9K 2 + 6V2~/C + 3 *1 R M = 0J2d V2/c(/c2+ V2K+ 1) V2/d>2 + V2K + l ) 3 1 U 2 + V2*+ I ; (3.36) Rv(y) and ^v(") are convenient in comparing controllers as they are independent of Kp and a>d-Figure 3.4 shows Rv(y) and ^v(") for CL and a PI controller that makes the closed loop damping factor n = V2/2. 3.4.5 Compar ison to P I Contro l ler A phase lag network necessary to implement the optimal controller CL = Kc(s+b)/(s+a) is usually not readily available in present day control system hardware, and it must be base loaded with the mean incoming flow Fm. On the other hand, PI controllers are a standard feature of every process control system. Also they don't have to be base-loaded with Fm. In short, the PI controller has a 3.4. Stationary Input 30 Figure 3.4: Rv(y) (decreasing graphs) and Rv(ii) (increasing graphs) for CL and a PI controller with r\ = 0.7. x-axis is K = (oc/<o<j. practical advantage over Cx. In this section, the performance of a PI controller is compared to CL. Let Cpj denote a PI controller, Trs + 1 ^ s + T~rx Cpr — Kc Trs = Kr (3.37) where K c 2 is the proportional gain and T r reset time. The sensitivity function S is S = (3.38) s 2 + K p K c s + K P K C T; 1 s 2 + 2noj c s + co 2 ' where OJc = ^KpKcT~l and n = KpKc/(2coc) is the damping factor. As before, OJC is parameterized as OJc = KiOd- After some calculations (details are in Section 3.8.3), the standardized variance ratios Rv(y) and Rv(u) are obtained. 1 ~>r>K + 1 ^ (3.39) Rv(y) = Rv(ii) = 2T}K(K2 + 2rj  + 1) y 3[(l+4?y 2> + 87?3] . . 2T)K(K2 + 2r]K + 1) 2The same symbol Kc is used for both Q and C>/ as this does not cause any confusion. 3.4. Stationary Input 31 If the optimum tuning of the PI controller is sought, it becomes the minimization of Rv(u) subject to Rv(y) = c\ (a positive constant c\ is the design specification). Rv(y) = c\ means 2T]K(K2 + 2TJK + 1) = 1/cj. Since 2VK(K2 + 2nK + 1) is the denominator of Rv(u), the minimization of Rv(u) becomes minimization of its numerator K3[(1 + 4n2)/e + 8n3]. So the optimum PI controller is characterized by a pair (K*, if): (K*,V*) = argmin{/c3[(l + 4n2)K + 8n3]) subject to 2TJK(K2 + 2TJK + 1) = l/c{. This is a mathematically tractable problem at least numerically. However, since (K*, n*) is a func-tion of Rv(y), the solution may be impractical. Practically more meaningful cases would result when the PI controller is tuned to a "moderate" (meaning that the closed loop system is not too oscillatory nor too slow) mode, say n between 0.5 and 2. When the input disturbance is a random walk, the optimum controller is a PI controller with n = V2/2 as shown in Section 3.6. Also the critical damping rj = 1 is used often. These two values of n and n = 2 for slow mode, are used for comparison. Figure 3.5 shows the ratio of ^v(") by the PI controller and Rv(u) by CL- The performance of the PI controller deteriorates when Rv(y) = Rv(y)a)j/K2 is high. Thus CL is worth consideration for cases of high Rv(y). 2.8 -2.6 -10"3 1<f2 10"' 10° 10' r„(y) Figure 3.5: Comparison of CL and the PI controller with n = V2/2, 1, and 2. The x-axis is Rv(y) and the y-axis is Rv(u) of the PI controller divided by ^ v(") of CL-3.5. Minimum Variance of Outlet Flow 32 3.5 Minimum Variance of Outlet Flow For some types of downstream processes, the adverse effect of Fu(f) is proportional to the deviation of Fu(t) from its mean value. Then the natural choice of the performance index is J=E[Xi(t)2+p2u(t)2]. (3-40) Now, u(i) is no longer included in J, so the system equation is constructed with u(t) as a control input. d_ dt Xl(t) 0 Kp X\(t) -Kp u(i) + 0 = + x2(t) 0 -cod Xlit) 0 0)d w(t) (3.41) A B The optimal control w*(r)is a linear state feedback with the feedback gain K =[k\ k 2], u(s) = -kixi(s) - k2x2(s) With the noise-free observer as before (x2(s) = KpSX\(s) + u(s)), (3.42) u*(s) = x\(s)[-ki - k 2 K p S ] - k 2u(s) u*(s) = - sk2K-1 + k 2 where c 0 = - --x\(s) = (c0 + cis)xi(s), (3.43) (3.44) l+fe ' " l + k 2 This is a proportional plus derivative (PD) controller. Let P = {pij} denote the solution of the ARE. Then, the feedback gain K is K = [kx k2] = BTP/p2 = [-Kppn/p2 - KpPl2/p2] The ARE is solved analytically with the following result. Pn = -7T, P\2 = Kp pcod + Kp (3.45) (3.46) Then Kpp +iod 1 co = , Ci = (3.47) pOJd PUd The controller CPD(s) = Co + C{S is not proper, and is implemented with a low-pass filter a/(s + a) to make the controller proper as a s + b CPD{s) - ci(s;+ CQ/CI ) = K r . , (3.48) s + a s + a where K c = ac\ and b = c 0 /ci . Thus, the PD controller is implemented as a lead-lag controller. 3.5. Minimum Variance of Outlet Flow 33 3.5.1 Comparison to P Controller In Shinskey (1994), the proportional only (P) controller is recommended for averaging level con-trol. Although flow smoothness was not mathematically defined, it can be understood from the context that the author sought to reduce the variance of u(t). In the following, the PD controller is compared to P controllers. Although, the PD controller is implemented as a lead-lag network, its performance is evaluated as CPD(s) = CQ + c\s to simplify the calculation and to provide the best case scenario for the PD controllers. PD Controller Performance The sensitivity function S is S = 1 = s = 1 . s ( 3 49) 1 + (c0 + cxs)Kp/s 5(1 + Kpci) + KpC0 1 + Kpcx s + \+Kpc\ It is convenient to define a parameter o)c as From (3.47), So coc = -^. (3.50) P OJc+OJd , KpC0 l + K p c \ = , and — = OJC. (3.51) ojd 1 + Kpc\ S = -J* i _ . (3.52) COc + OJd S + 0)c As before, define K as o>c = KOJJ, then S = — — . (3.53) 1 + K S + KOJd After some calculations, Rv(y) is obtained. The variance ratio between Var[w(r)] and Vax[d(t)], Rv(u) is VarKQ] ir^ + Sic+l) 3.5. Minimum Variance of Outlet Flow 34 P Controller The transfer function of the P controller C p(s) is •Cp — Kr (3.56) Then the sensitivity function with this controller is 1 S = 1 + KcKp/s s + KCKP S + OJC (3.57) where coc = KCKP = KOJJ. Then, Rv(y) and Rv(u) for Cp are obtained. Rv(y) = OJIK(K+ 1) Rv(y) = 1 K(K+ l) Rv(u) = K+ 1 (3.58) (3.59) The reduction of R v (u) by C P D compared to C P is plotted against Rv(y) in Figure 3.6. The maximum reduction is 0.8638 when Rv(y) = 0.7077. This reduction is 0.93 in terms of standard deviation. So the PD controller reduces the standard deviation of u{t) by 7 % compared to the P controller at best. In reality, the reduction is less thanr7 % due to the low-pass filter. This small reduction probably means that for practical purposes, the PD controller is no better than the P controller. Figure 3.6: The ratio of the variance of u{f) between Cpo and Cp. x-axis is Rv(y). 3.6. Random Walk Input 35 3.6 Random Walk Input Although stationary processes are preferred as a model of the input disturbance, the random-walk case is treated here for two reasons: (a) to show that the optimal controller is the same PI controller for deterministic cases where the input disturbance is a step, and (b) to show a proper way to choose the state variables for a non-stationary disturbance. Doss et al. (1983) considered the control of liquid levels in three tanks in series by a state-space LQ approach without explicitly including the input disturbance in the state , but resorted to other more intuitive approaches that they felt led to a simpler, more physically acceptable scheme. Harris and MacGregor (1987) reworked the same problem, and showed that the same satisfactory result can be obtained by the transfer function method with a disturbance model. If the input disturbance had been included in Doss's work as shown below, the state-space approach would have produced a satisfactory result. A continuous-time random walk is a Wiener process w(t). Since w(f) is non-stationary, the difference between the incoming flow and the outgoing flow, which is stationary, is chosen as a state variable. With a Wiener process w(t) as the input disturbance, the manipulated variable u(i) = Fu(t), and the tank level error x\(t) = y(t) -yr, the process is expressed as = xi(0) + K p f (w(r)-u(T))dT. Jo With X2(t) = w(t) - u(t), w(f) = dw(t)/dt, and v(t) = u(t), the system equation is (3.60) d_ dt X\(t) 0 Kp *i(0 0 v(0 + 0 = + x2(i) 0 0 x2(t) -1 1 w(t). (3.61) B The performance index J is J := E[xi(t)2 + p2v(02] = E{x(t)T x(t) + v(0 v(0} (3.62) Note that the outlet flow u(t) itself cannot be included in the performance index because it is non-stationary and may not have a well-defined variance. The optimum control v(t)* that minimizes / is a feedback form v(0* = -Kx(t) (3.63) 3.6. Random Walk Input 36 As before, the solution of the ARE, P = {prf is obtained analytically. Substitution of A,B,Q, and R of (3.61) and (3.62) into (3.6) yields the following three equations. -p22/p2 + \ = 0 PnKp-pnP22/p2 = 0 (3.64) IpnKp - p\2lp2 = 0 The above three equations and the positive definiteness of P determine ptj as follows. P n = p , P 2 2 = P y j 2 p K p p n = ^ 2 p / K p And (3.65) (3.66) K = R-lBrP = [-pl2/p2 -P22/P2] = [-l/p - y/2p~Kp/p] Since u ( t ) = v(t) = - K x , su(s) = xi(s)/p + x2(s) yjlpKpIp (3.67) Usually x2(t) is not available for measurement and an observer is necessary. A simple noise-free observer is designed below. Since x\(t) = Kpx2(t), in Laplace transform, sx\(s) = Kpx2(s). Thus, x2(s) = SX\ (s) (3.68) Substitution of (3.68) into (3.67) yields su(s) = X\(s) Thus fl ^2pKP\ ~ + s—-p— = xi(s) u(s) = X\(s) 1 — + (3.69) (3.70) sp ypKp This is a PI controller with the proportional gain Kc = yf^, and the reset time Tr = ^2p/Kp. The sensitivity function S = 1/(1 + PC) is 1 S = l + * z ( ± + / X \ s 2 + Syl2Kp7^ + Kp/p S2 + 2^U>CS + G J I (3.71) where OJC = •\fKp~/p. So the optimum controller is a PI controller that makes the closed loop system a second order low-pass filter with a damping factor 0.5 VI. 3.6. Random Walk Input 37 3.6.1 Effects of Damping Factor Some PI tuning rules recommend the damping factor n = 1 for integrating processes , for example Morari and Zafiriou (1989). Effects of n on Var[v(/)] was investigated when n is deviated from the optimum value rf = 0.5 V2. The PI controller is tuned so that the sensitivity function becomes s 2 S = -^—z =• (3.72) s1 + 2T]0)cS + (x>Lc Then Var[y(r)] and Var[v(f)] become (denoting the variance of w(t) with cr^), Var|M0] = S . Var[v(0] = ^ I ' c (3.73) For a given Var[y(J)] target, o)c is adjusted for each n ± n* to keep Var[y(0] to its target. Then Var[v(£)] is compared to the optimum value Var*[v(0], Var[v(0] (4n2 + l ) 2 V 2 ( 1 x l / 3 Var*[v(0] 4n l V 2 ^ J (3.74) This is plotted for 77 = [0.5 1.0] in Figure 3.7. Var[v(f)] increases by 5% when 7/ = 1. 1.06H r\ (damping factor) .1 r Figure 3.7: Comparison of the optimum controller (n = V2/2) and non-optimum (n = [0.5,1]). X-axis is for the damping factor n and Y-axis for the variance of u divided by the optimum variance. 3.7. Summary 38 3.7 Summary The averaging level control problem is formulated as a minimization of Var[w] subject to the target tank level variance Var[y(f)]. The problem is solved by the state-space method combined with a noise-free observer. When the input disturbance d(t) is modelled as a stationary random pro-cess (output of a first-order low-pass filter with a cutoff frequency ojd driven by white noise), the optimum controller CL is a phase lag network, Ci = Kc . s + a When the input disturbance is modelled as a random walk, the optimum controller is a PI controller, Cpi = Kc . s For both Q, and C/>/, the closed loop system is a second-order system with damping factor V2/2. A single parameter Rv(y) that characterizes the controller is introduced, Var[y(Q] (ojd\2 For a stationary input disturbance, Ci and CPI are compared by the ratio of Var[zi(r)] with the two controllers set up to give the same tank level variance. CL is not much different from CPI when Rv(y) is small (below 1), but CL performs better for higher Rv(y). For instance, for Rv(y) = 10, CL's Var[w(0] is 1/2.6 of that of CPI. When the flow smoothness is represented by the variance of u(t) in stead of u(f), the opti-mum controller, which minimizes Var[u(t)] with a target tank level variance, is a proportional plus derivative (PD) controller, Cpo, CpD = CQ + C \ S . However, the performance difference between CPD and a proportional only controller CP is small (7 % in standard deviation of u(t) at best), and for practical purposes, Cp is recommended. 3.8. Mathematical Details 39 3.8 Mathematical Details 3.8.1 Solution of Section 3.4.2 The A R E yields the following seven equations. («) -P2n/P2 + 1 = 0 (b) pnKp - pX2Ab - pnpnlp1 = 0 (c) - pnKp - pnpislp2 = 0 ( d ) 2 p n K p - 2 p 2 2 A b - p \ 3 i p 2 = 0 ( e ) - p l 2 K P + p l 3 K p - p 2 , A b - p 2 3 p 3 3 / p 2 ( J ) - 2 p l 3 K p - p 2 3 / p 2 + p p 2 = 0 First from (a), p \ 3 = ± p , but from (J), p \ 3 must be - p to make p 3 3 real. (a) - p 2 1 3 /p 2 + 1 = 0 => p l 3 = -p (/) - 2p l 3Kp - p\ 3lp 2 +pp 2 = 0 => p i 3 = p yJ2pK p +pp 2 Then using (c), (c) -pnKp-pnp 3 3/p 2 = 0 => - /?nA> + ^2pK p + pp2 = 0 => /?n = ^2plK p +p(p/K p) 2 Using (6), ( 6 ) p u K p - p n h - p n P i - i l p 2 = 0 = > p l 2 = ^ l ( p n K p - p u p 2 3 / p 2 ) = A ; \ ^ j 2 p K p + p p 2 + p 2 i / p ) Then, substitute the known to (e), (e) - p l 2 K p + / ^ .K /? - p n h - P2iPnlpV=> - p u K p -pKp - p 2 i A b - p 2 3 y J 2 p K p + p p 2 / p So Kp(p + A~bl yj2pKp+pp2) Kp(Abp + y/2pKp + pp2) p 2 3 = = Ab + ^/2Kp/p+p + KpA-blp-x A2 + Ab ^/2Kp/p+p + Kpp~l 3.8.2 Rv(y) and R v (u) in 3.4.4 As the sensitivity function 3.8. Mathematical Details 40 the transfer functions from d(s) to x\(s) GDY and from d(s) to ii(s) G D y are K p ( s + d) GDY = GDU = s2 + y/2o)cs + co2 K p K c ( s 2 + bs) s2 + ^j2cocs + co2 Thus Var[x!(0] = /v 1 Kp(s + a) = K2pIv s + a [s + cod s2 + y/2cocs + co2c) We express coc as Kcod, then a, b, and Kc are expressed as follows: ,,2 U 3 + *2(o>rf+ V2o»c) + s( yf2coccod+co2) + £4/o>2, a = >C + A/2/C K2 + V2A: + 1 COD = KA(OD < 0)d K2 + V2/c + 1 b = iOd = KbO)d > 0)d 1 + V2K Var[x,(0]=^/v s + a = K: [s3 + s2(l + K V2)w r f + S(K V2 + /c2)^2 + K 2 ^ , 2 AC 2 W 3 +^( l - r / cV2)^ 2K2co6d ( ( 1 + K V2)(/f V2 + /c2) - K2) 1 + ' P2w 3 V2/rt>2-t- V2/c+ 1) Since Var[d(t)] = Iv(l/(s + cod)) = \/(2cod), _ Var[x,(Q] _ 2 Rv^ ~ Verity] ~ K l + J ^ ^ ^ l + K ^ _K2P * 4 + 3 V 2 « 3 + 9tt2 + 6 V 2 « + 3 * co2. V2K(K 2 + V2K+ 1) w2 V2K(K2 + V2K+ l ) 3 So *v(v) := Rv(yWd K4 + 3 V2K3 + 9K2 + 6 V2K + 3 *2 V2K(K2 + V2/f+ l )3 (3.75) Similarly for Var[w(0] and Rv(u), ( 1 KpKc(s2+bs) ) Var[zi(r)] = 7V + 5 2 + yJ2coCS + C02) s 2 + bs U 3 + s2(cod + V2fa>c) + s( ^[2coccod + co2) + codco2 Since K P K C = - k 2 , KpKc — y 2(l + V2/c) V2/c+ 1 3 Here, the variance ratio is concerned, so the input variance can be set to any convenient value. 3.8. Mathematical Details 41 So, Var[zi(0] = ( K\\ + V2K) ^ OJd 2 + Vz* + 1 ) s 2 + bs U 3 + S2(\ + K V2)it»d + S(K V2 + K2)CJ2 + K2OJ. Substituting b = ^j-^^^d, and after some algebra, we get, Rv(u) = 0>A V2K(K2 + V2/f + 1) + V5#c + u *v(u) := /?V(M) 1 w2, ~ V2K(K2 + V2* +1) V ( l + V2*) ,K2+ V2K + l j (K2 + V2K) + (3.76) 3.8.3 R v(y) and i?v(w) in 3.4.5 _1 frS s + T~L . s + b Then Thus, CPI(s) = K c \ \ + ^ - = K c ~ — ^ - = K c -S = b = T:1 1 + Ikri+b s 2 + sK p K c + K p K c b s 2 + 2rjojcs + KPKC = 2na>c, and KpKcb = OJ2 Substitution of K P K C - 2T]OJC into K p K c b = o 2 yields So, G n v = GDU = Kps Kps s 2 + 2no)cs + OJ2 s 2 + 2r)KO)dS + K2to2d KpKc(s2 + bs) 2 n o j c ( s 2 + sojc/2n) 2nKOJds2 + K2o?ds S2 + 2r\(x>cS + QJ>2 S2 + 2j]KCOdS + K2W2d S2 + 2r}K0)dS + K2OJ2 Var[xi(0] = /v ( — , , „ ^ , , , ) = K2L s + ojd s2 + 2nKiOdS + K2cod j ^ v \5- 3 + ^2(1 +2r)K)a>d + s(2nK+K2)oJd + K20Jd, 2 ((1 +2i7/t)(2^K+«r2)w5 - K2ufy 2<uj2 i^r(»r2 + 2^+1) 3.8. Mathematical Details 42 Rv(y) = Var[xi(0]2wrf = cod2r]K(K2 + 2nK + 1) Rv(y) := Rv(y)co2d 1 2TJK(K2 + 2TJK + 1) (3.77) Var[w(0] = Iv 1 2nKu>dS2 + K2u)ds s + cod s 2 + 2r\K0)ds + K2O>2 2rjKO>dS2 + K2co2ds • d , \S3 + s 2(\ +2nK)cod + s(2rjK+K2)a>d + K2cod; {2t]Kcod)2{2T}K + K 2) + K4G>4 ^ 1 + 4n2)/c + St]3) 2cod2r]K(K2 + 2T]K +1) 4TJK(K2 + 2TJK + 1) R v ( u ) = Yar[u(t)]2tOd = OJ2K\(1+4?]2)K + 8T]3) 2TJK(K2 + 2nK + 1) Chapter 4 Downstream Processes 4.1 Introduction In the previous chapter, the smoothness of the outlet flow Fu(t) is mathematically captured as the variance of Fu(t) = v(t) and/or Fu(t), and the problem is solved by optimization with the following performance index (3.15), J\ = VarLMO] +p2 (Var[v(0] +pVar[U(0]), wherey\(t) is the tank level and u(t) is the deviation of Fu(t) from its mean value Fm. The reason for requiring outlet flow smoothness is to minimize adverse effects of tank level control on the downstream processes. In this chapter, the downstream process variable is explicitly included in the problem formulation. The averaging level control problem is considered as the minimization of the downstream variable variance while keeping the tank level between high and low limits. Let y 2 ( t ) denote the downstream process variable. The minimization of Var[y2(r)] makes the problem mathematically tractable, and also it makes practical sense for some cases. For example, y 2 ( t ) may represent a product quality parameter, and y2(t) must be in a specified range for the product to be acceptable. Minimizing Var[y2(0L then maximizes the product yield. In this chapter, the following performance index is used. J := Var[y,(0] +P^ Var[y2(0] +p2Var[v(0] (4.1) As in Chapter 3, pi and p 2 serve as Lagrange multipliers so that Var[yi(0] is minimized subject to Pi Var[v(0] + plVar\y2(t)] = c, where c is some constant. Actually, c is not given, and pi and p 2 are 43 4.2. Downstream Process Model 44 adjusted to obtain a required Var[yi(r)], which is specified as one of the design parameters. Since the problems are treated in continuous time, non-zero p\ is necessary to keep the feedback gain K and Var[v(r)] finite. 4.2 Downstream Process Ivlodel Figure 4.1 depicts a model of a combined system of the tank level control loop and the downstream process, where the outlet flow Fu(t) enters at the output of the downstream process P2. The c2 Fd(i)-Ci Fu(t) 6-r*y2(t) Figure 4.1: The tank level control loop (Ci and P\) and the downstream process P2. The outlet flow Fu(t) enters at the output of the P2. purpose of the model is to investigate how the downstream process influences tank level controller design when the statistical relation between Fu(t) and y2(f) is taken into account. The model is not meant for designing a multivariable system, where Fu(t) is determined from both y\(t) and y2(t). Thus, the interaction between Fu(t) and y2(t) are simplified excluding any deadtime or gain units, which may exist between Fu(t) and y2(t). Therefore, it is assumed that the downstream process variable y2(t) is not available to the tank level control loop. If Fu(f) affects more than one downstream processes, y2(t) is a hypothetical "representative" process variable that has no real physical entity. Downstream processes are represented here by a first-order plus deadtime model. K2 Pi = -e~TS = K2-T2S +1 S + 0>2 where K 2 is the steady state gain, T 2 the time constant, co 2 = T 2 the cutoff frequency and T the 4.2. Downstream Process Model 45 deadtime. It is assumed that P2 is closed with a PI controller C2, C2 = Kc ( l + -U = KCS-^-, (4.2) \ Trsj s where Kc is the proportional gain, Tr the integrating time, and cor = l/Tr. Then the sensitivity functions'2 is o 1 <*+<Ol) ( 4 3 ) 1 + P2C2 s(s + io2) + K2to2Kc(s + cor)e-TS The IMC tuning rules (Morari and Zafiriou, 1989) recommend either Tr = T2 or Tr = T2 + 0.5r and (Ogawa, 1995) Tr ~ T2 + 0.5T for robust performance. Here, Tr is set to T2, that is cor = co2, to simplify the analysis. Then, S2 becomes S2 = . (4.4) .. s,+ K2Kco)2e-TS ' S2 = \ / ( l + P2C2) is the transfer function from u(s) to y2(s) when u(s) enters at the output of the downstream process, The transfer function from ii(s) = v(^ ) toy2 is S2/s, which is a low-pass filter. When the deadtime T is small, S2/s can be approximated by the first-order filter. GV2:=^*—!—, (4.5) v s + co2c where co2c = K2Kcco2 is the closed loop cutoff frequency of the downstream process. For robust performance, Kc is set to make KCK2 between 0.25 and 0.41 , that is, co2c = 0.25 ~ 0.4w2. (4.6) When the deadtime is not negligible, e~TS can be represented by a Pade approximation, _T, 1 - 0.5™ ioT-s 2 1 + 0.5TS COT + S T Then G n L _ L i ! i = s + \ . (4.8) S + K ^ C J ^ j ^ S2 + S(C0T - C02c) + L02cLOT This is a low-pass filter with a cut-off frequency = ^co2ccoT. The gain plot has -40 dB/decade slope between to'p and wT, and -20 dB/decade slope beyond coT. This can be approximated by a first-order low-pass filter with a cut-off frequency cop chosen between co'p and coT. GV2 ~ —-—, yJoj2ccoT <cop< coT. (4.9) S + COp 4.2. Downstream Process Model 46 20 10 0 -10 CD T3 -20 -30 -40 -50 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 log(to) Figure 4.2: The magnitude of Gvi (solid line) and two approximations for io2 - K2 = r = 1 and Kc - 0.3. The dashed line is for OJp = io'p and dotted line is for a>p = l.4a>'p. Figure 4.2 shows the true Gvi and two approximations with different choice of a>p. When u enters at the input of the downstream process, Gvi has an additional turn-over fre-quency a>2 from which \Gvi\ decreases with -40dB/decade. If such a model is used, the resulting controller, which is more complex, will produce outputs that have more high frequency compo-nents. Therefore, in the sequel, the downstream process is assumed to be a simple low-pass filter driven by v(t). Let x3(t) denote the deviation of^ (0 from its mean value: xi(t):=y2{f)-E\y2(t)l (4.10) Since x3(t) is the output of a low-pass filter driven by v(r), x3(i) = -ojpx3(t) + v(t). (4.11) In reality, there is some "gain" involved between v(t) and x3(t) depending on where and how u(t) enters the downstream process. Then the downstream model is x 3 ( s ) = - ^ - v ( s ) , (4.12) S + (Op where Kj is the gain. The downstream process variable x3(t) is assumed not available for the tank level control loop, and recovered by a noise-free observer from xi(t) and w(r) using the downstream process model (4.12). Therefore, Nc(s) 'The exact value depends on the plant uncertainty and T (Ogawa, 1995), but the IMC rules recommend similar values (Morari and Zafiriou, 1989) • ' • 4.3. State-Space Solution 47 where CD(S) is the controller transfer function. CD is calculated to attain required Var[yi(0] and Var[v(r)] with the minimum Var[y2] by adjusting pi and p 2 in (4.1). Let CD denote the optimum controller for Kj = 1, and C'D for Kd ± \, both with the same Var[yi(r)] and Var[v(f)]. Let cr2, denote Var[y2(r)] f ° r Kd = 1 and o\ for Kd ± \. From (4.12), cr2 = K^a2 if the same controller is used for both cases. If CD i1 C'D, cr2 ± Kda2. This contradicts the fact that CD and C'D are the optimum controllers for each case. Thus, CD must be equal to C'D. This means that Kd does not influence the optimum controller. However it does influence the value of p 2 to obtain the desired Var[yi(r)] and Var[v(r)]. Thus, Kd is set to a convenient value for calculation. 4.3 State-Space Solution The downstream process variable x->,{f) is added to the state variables of Chapter 3. x \ (0 = y\ (0 ~~ yr tank level error x2(0 = Fd(i) - Fm input disturbance deviation *3(0 = yi(f) ~ E[y2(0] downstream deviation variable u(t) = Fu(t) - Fm outgoing flow deviation (manipulated variable) In this chapter, the input disturbance is assumed to be a stationary process with the first-order low-pass spectrum. x2(t) = -ojdXiit) + toaMt). (4.13) Figure 4.3 shows the structure of the combined system of the tank level control loop and the white noise —H C bid s + a>d J ~ T " T ' J i T cop s + a>n u x2 Kn Figure 4.3: The system structure for the tank level control loop and the downstream process. 4.3. State-Space Solution 48 downstream process, where the downstream process variable x3 is controlled in an open-loop man-ner from v. The controller, C, calculates the formal manipulated variable v from x\ , x 2 , x 3 and u, but x2 and x3 are not measured. Thus, they are not shown as inputs to C. With (4.11) and (4.13), the following system equation is obtained. Xl(t) d x2(t) dt u(t) 0 Kp 0 -cod 0 0 0 0 0 0 -ojr o p'" 0 X\(t) 0 0 x2{t) 0 v(t) + COd + xiii) 1 0 u(t) 1 0 Mt) (4.14) The performance index is / = E[xi(02 + p2v(02 + P2*3(02] = E[x(t)TQx(t) + Rv(t)2], (4.15) and R = [p]]. (4.16) where r 1 0 0 0 0 0 0 0 0 0 p \ 0 0 0 0 0 In continuous time, p\ must be greater than/zero for the problem to have a proper solution. As ' j \ ' usual, the positive-definite solution of the ARE, P and the feedback gain K = [k\ k2 k3 k*\ are !1 • • obtained. For a problem of this size, it is difficult (may not be impossible) to obtain closed form analytical solutions. Thus, all the analyses in this chapter rely on numerical solutions. The optimal control is a linear state feedback form, su(s) = -k\X\(s) - k2x2(s) - faxiis) - k4u(s) With the noise-free observer, (4.17) xi(s) = x2(s) = Kp sx\(s) + u(s) S + 0), •u(s) su(s) = - & i X i ( s ) - k2[K^sxiis) + u(s)] - k3 S + 0)B -u(s) - k^u(s) (4.18) (4.19) 4.4. Performance Comparison 49 Thus the controller denoted with CD(S) is -(s + iop)(k2K-x s + ki) CD(s) = s2 + (ojp + k2 + h + k4)s + ojp(k2 + k4) ^ b2s2 + b\s + bo s2 + a\s + ao Since CD(s) includes a factor (s + OJP), which is due to the noise-free observer, in its numerator, it cancels the pole in the downstream process 2 . The sensitivity function S is s = S(s2 + alS + ao) s3 + s2(a{ + Kpb2) + s(a0 + Kpbx) + Kpb0' The transfer functions from x2(s) to x\(s), GDYX , to v(s), GDV and to x3(s), GDY2 are Kp(s2 + a\s + ao) GDYX = 5 3 + s2(ai + Kpb2) + s(a0 + Kpb{) + Kpb0 Q = sKpjb2s2 + bxs + 6p) D V s3 + s2(ai + Kpb2) + s(a0 + Kpbx) + Kpb0 ' sKp(b2s2 + bis + b0) 1 GDY2 s3 + s2{ai + Kpb2) + s(a0 + Kpb\) + Kpb0 s + a>p From the above transfer functions, the variances of y\(i), v(r), and y2{t) driven by Fd(i) are obtained using the integrals of (3.11). Since Var[Frf(f)] = a>d/2, the variance ratios Rv(-) are obtained by dividing the variances by a>d/2. For example, R ^ i ) = I(_^£ Kp(s2 + alS + a0) \ 2_ 1 v\s + CL)d s3 + s2(a\ + Kpb2) + s{ao + Kpbx) + Kpbo) ojd' 4.4 Performance Comparison In this section, the performance of Co (4.20) is compared to that of a controller designed with no consideration given to the downstream process. The controller, denoted with d, that minimizes the performance index JL = Var[yi(r)] + pfVar[v(0] is compared to CD. To compare on an equal basis, CD and CL are tuned so that the two, controllers yield the same tank level variance. The main concern of comparison is the reduction of Var[y2(0] by CD compared to CL- We denote the variance 2 In the Wiener-Hopf method, the reason (s + OJP) shows up in the controller transfer function numerator is obscured by complicated algebra - another advantage of the state-space approach. 4.4. Performance Comparison 50 ratios under a particular controller by subscripts, for instance, Rv(y2) under CD(S) is denoted with R-v(yi)cD- The reduction of the 3/2 (0 variance due to CD is denoted with rv(y2)5 that is, r M : . * ^ S . . (4.24) 4.4.1 Effects of v(0 weight Since CD is designed with the performance index J = Var[yi(Y)] + p2Var[v(/)] + p^Varlj^COL for a given target Rv(yi), different manipulated variable variance, Rv(v), is obtained by adjusting pi and P2. As CL minimizes Var[v(r)] for a given Var[yi(r)], V a r K O k ^ V a r K O ] ^ . The Var[v(r)] increase by CD over CL is denoted with rv(v), r " ( v ) : = p / \ • ( 4 - 2 5 ) ^v(v) C i Intuitively speaking, Co reduces Var[y2(0] compared to Cx by increasing Var[v(/)]. The down-stream process is a low-pass filter driven by v(t). Thus, high frequency components of v(t) perturb 72(0 very little. Therefore, Co shifts control energy to higher frequencies to reduce Var[>>2(0]-Var[y2(0] would be reduced more for larger rv(v) with large high frequency components. The gains of Co and Q are compared as shown in Figure 4.4. CD has higher gains at high frequency ranges compared to CL and lower gains at lower frequencies. Figure 4.5 shows the magnitude of the sensitivity function S, \S\, for Co with r„(v) = 1.2 and 2, and for CL. \S\ for CD has lower magnitudes around the cut-off frequency cop = 1 of the downstream process, and small increases at lower frequencies as a result of CD shifting the control energy towards higher frequencies. This is the case as shown in Figure 4.6, where rv(y>2) steadily decreases to 0.7068 for rv(v) = 100. In reality, the reduction is less because the dynamics of a valve actuator has a limited bandwidth, and high frequency components of v(f) cannot influence the tank level y\(t). The valve actuator dynamics is not included in the mathematical model to simplify the derivations. As rv(y2) = 0.75 at rv(v) = 2, a good portion of possible Rv(yi) reduction is obtained with a modest increase in Rv(v). 4.4. Performance Comparison 51 CO Figure 4.4: Gains of Q> and CL- Both controllers are tuned so that Var[yi(r)] are the same. The graphs marked with 2 and 1.2 are the gain for CD tuned tp make Var[v(0]cD equal to 2Var[v(t)]cL and 1.2Var[v(t)]cL respectively. The system parameters are: Kp =.cod = o>p = 1. The controllers are tuned for Rv(y\) = 1. 5 m 13 CO Figure 4.5: Magnitude of S for CD and CL- The controller tunings are the same as for Figure 4.4. 4.4. Performance Comparison 52 1''' Rv(du) ratio Figure 4.6: Rv(y2) reduction rv(y2) in y-axis and rv(v) in x-axis. The system parameters are: Kp - ojd -(jjp-\. The controllers are tuned for Rv(y\) - 1. 4.4.2 Var[y2j reduction versus o>p The i? v(v2) reduction, rv(y2), by CD compared to the minimum Var[v(r)] controller CL is plotted for ojp/ojd between 0.01 and 100 in Figure 4.7, parameterized by five different normalized tank level variance, Rv(yi), between 0.001 and 100. The downstream process cut-off frequency OJP Figure 4.7: rv(y2)for5 different Rv(y\) of 0.01,0.1, 1, 10, and 100, which are denoted on the graph. 4.4. Performance Comparison 53 is measured relative to the incoming disturbance bandwidth cod, thus, jc-axis of the graph is for cjp/cod- The manipulated variable variance increase, rv(v) is set to 2. When rv(y) =10, the graphs are very similar with rv(y2) ranges between 0.6 and 0.85 for iop/o)d = 0.01. The graphs shift to left for increasing Rv(yi). Since \/Rv(y\) is proportional the bandwidth of the level control loop (the loop must be fast to make the level variance small), the downstream cut-off frequency top is measured against -\J\/Rv(yi) (inverse of the standard deviation). Let o>n denote the normalized downstream cut-off frequency relative to o>d and l/i?v(vi), (On •= — J&v(yi). Figure 4.8 plots the same graphs of Figure 4.7 with the x-axis for the normalized downstream frequency u>„. Figure 4.8 demonstrates that con is a convenient parameter for evaluating the effect 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 10"3 10"2 10~1 10° 101 102 103 Figure 4.8: The same graphs as Figure 4.7 with x-axis for o)„ = u>p/(Od V^vCji)-of the downstream process bandwidth o>p on rv(y2)-1. The asymptotic value of rv(y2) towards cjn = 0 is smaller for higher Rv(yi). Thus, when the downstream process is slow (relative to the incoming disturbance and the level control loop), and the level loop is slow, Q> reduces the downstream process variable variance as much as 35% from Q . 4.5. Summary 54 2. When the downstream process is fast, say OJ„ > 10, rv(y2) > 0.9 for all values of Rv(yi). Thus, CD is not effective compared to Q . Therefore, when the downstream process is fast, Var[v(r)] is an appropriate measure of the outlet flow roughness. 4.5 Summary The downstream process is modelled as a first-order low-pass filter from v(r) to its process variable yi(t), v yi(s) = ——v(s). S + dip The optimum linear controller CD that minimizes Var[y2(0] f° r a target tank level variance Var[yi(r)] is obtained when the input disturbance is a stationary process with the cutoff frequency at o)d- CD is a second order system, (s + (op)(s + b) s1 + a\s + ao It is noted that CD cancels a pole in the downstream process. The best case downstream process variable variance reduction compared to CL ranges 0.8 to 0.6 depending on the tank level loop tuning Rv(yi) and the relative downstream frequency ojn = ojp/ojd ^Rv(y\). CD is not much better than CL when the downstream process is fast or the tank level control loop is tight. Chapter 5 Markov Jump Systems 5.1 Introduction Jump Markov systems, where sudden changes in system parameters and disturbances, or com-ponent failures are modelled by continuous-time Markov processes, can model a wide range of systems not only in engineering but also social and economic spheres. Sheet breaks on paper ma-chines can be modelled by a continuous time Markov process as shown in (Khanbaghi et al, 1997) and later in Section 7.2.4 of this thesis. The input disturbance to a broke tank level control loop can be modelled as a finite state (2 or 3) continuous time Markov process. This makes it possible to ap-ply the jump linear quadratic control method to broke tank level control (Khanbaghi et al, 2002). It is shown later that the jump linear quadratic control method applied to broke tank level control has no special advantage over the normal linear quadratic control approach. However, the investi-gation of jump Markov systems leads to a reliable numerical method to calculate the distribution of state variables without simulations. This is useful for broke tank level control because the prob-ability of level violation for various controllers, linear or nonlinear, can be evaluated reliably. In this context, jump systems are included in the thesis. Thus, the emphasis is on a new theorem for state distributions. The investigation of jump Markov systems has been active since the late 1960s ((Sworder, 1969) and (Wonham, 1970)), and many studies have been published. Two areas, optimal control (unconstrained and constrained (Costa et al, 1999)) of linear systems and stability and control-lability aspects ((Mariton, 1990) and (Ji and Chizeck, 1990)), are the main subjects of published 55 5.2. Jump Markov Systems 56 studies \ However, very little has been published on the evaluation of performance indices or state variances, which is important for practical applications. From the nature of disturbances, the state distributions are most likely non-Gaussian. Therefore, the state variance may not provide enough information on control performance. To the author's knowledge, Sworder and Choi (1976) is the only study that investigated the distribution of performance indices under optimal control. Since a system undergoes sudden discontinuous changes (mode or regime change), it is desirable to know the performance index or the state and control variances for each mode so that the controller can be tuned for each mode. This chapter presents a new theorem that makes it possible to calculate the state variances and distributions for eaehihode. 5.2 Jump Markov Systems Let x(t) G RN x be the continuous state of a system, u(t) e R^" the control input. The system is subjected to sudden changes (jumps). Between jumps, the system is described by a time-invariant differential equation jx(t) = f(x(t),u(t),m = fm(x(t),u(t)), (5.1) where £(t) is a finite-state continuous-time Markov chain. Without loss of generality, the Markov states are expressed as integers between 1 and N%. Thus, £(t) e S , where a set S = {1,2,..., A^}. The system (5.1) is in mode £ when £(t) = €.ln this chapter, any entity that is a function of the finite state continuous time Markov process £(r) is expressed with a subscript I , e.g. f t := /(£(0)lfM=f-Selected elements of continuous-time Markov chain theory in Tijms (1994) and Stewart (1994), which are used in the thesis, are stated here. Continuous-time Markov chains can be characterized by the probabilities of state transitions in a short time period At, P[£(f+Af) = = ?']=1+ quAt + o(At), q„ < 0 (5.2) F[((t+M) = M(f) = i\ = qijAt + o(Ar), q l } > 0, where qu = - }Zj±i Qij and °(0 denotes a quantity such that lim^o o(t)/t = 0. As it is assumed that Markov chains are homogeneous, the <?,/s are not functions of time t . The matrix Q , whose i j -th element is qtj, is called the infinitesimal generator. It is assumed that is irreducible and has a 'The literature is vast, and only a small fraction of representative papers are cited in the thesis 5.3. Holding Time Density 57 unique stationary distribution n = [n\ 712... 7IN(] (A^-row vector), which satisfies the equation nQ = 0. (5.3) Note that a continuous-time Markov chain'with a stationary distribution is ergodic (Grimmett and Stirzaker, 1992). Another way of describing continuous-time Markov chains is by the transition probability and holding time distribution. To distinguish points in time and time durations, follow-ing Feller (1968), the word epoch is used to denote points on the time axis and time for durations. Let t„,n = 1,2,... be epochs when £(i) changes its state. We call tn the n-th jump epoch and Tn '•— tn - tn-\ (to = 0) the «-th switching epoch, that is, jump epochs have a fixed time origin at 0, and switching epochs have time origin shifted to the previous jump epoch. Time durations between jump epochs {r, 0 < T < Tn) are called holding times. • The switching epoch of £(i) in state /, Ti} is exponentially distributed with intensity Let Xi > 0 be - q u , then P [ 7 / < t ] = l - e q i i t = 1 - e~Xit. (5.4) • If £(t) leaves the present state i, it jumps to state j £ i with probability pij which is related to the qi/s as follows. Pij = and, * L = * i j * i -QU ^ (5.5) 0 j = i zZPij = 1. The matrix P whose (/'-th element is ptj satisfies conditions to be a transition probability matrix (rows sum up to 1, and all elements are non-negative). A discrete-time Markov chain whose transition probability matrix is P is called the embedded Markov chain of the continuous-time Markov chain £(t) (Stewart, 1994). The embedded Markov chain of £(t) plays an important role in the development of the new method. 5.3 Holding Time Density The holding time density hc(r) is introduced to describe the long-run average "distribution" of short time intervals included in holding times of state I. The material here is not new, but described in 5.3. Holding Time Density 58 detail because the concept is important for a new theorem introduced in the next section. The holding time density would be best understood by an example of its usage. Suppose that a random variable X(t) taking its value in R is a function of £(/) and how long £(t) stays at the present state, X(t)=fm,t-ts)=fr(T) (T = t - t s ) , (5.6) where t s is the latest jump epoch before t . The conditional distribution of X(t) given £(t) = I is obtainable if we know h i , where B is any Borel set and 1 denotes the indicator function of the inverse image of fi(B) = fi t , B) . Then, for instance, the conditional mean of X(t) given £(/) = t is calculated as E[*(0I£(0 = t] = f xfieidx). JR Conceptually the holding time density is similar to the histogram for holding times. Suppose a histogram of holding times for state t is built with bin width At, where relative frequencies that the holding times spend in each time intervals (bins), ([0, At), [At, 2At), [2At, 3At),...) are collected. The bins are numbered from 1, so the /-th bin covers a time interval [(/ - I)At, iAi). If a switching epoch of a sample is longer than (/ - 1)A/, bins 1 to / are incremented by 1. Let c,(0 denotes the /-th bin count up to epoch t. The relative frequency of the /-th bin is Since £(/) is ergodic, limr,(/) = E( r j= I he(f)dt a.s. (5.7) '- > 0° J(i-1)A/ For a short time interval [t, t + At) to be included in a holding time, the switching epoch must be longer than / + At. The probability that the switching epoch is longer than / + At is e~A^'+&t\ Thus h((t) is proportional to e~Ae', he(t) = ae~Xct, (5.8) where a is a normalizing constant. The constant a is determined such that ht(t)dt = 1 5.4. State Distribution 59 holds. Since JT°° e~Xe'dt = A~(x, h((t) = XteXlt. (5.9) From the definition of continuous time Markov chains (5.4), the switching epoch T(£) for state £ is exponentially distributed with intensity At, F[T(-C) < t] = 1 - e~X(t. Let st(t) be the probability density function of the switching epochs of state £, seiO := ^[T(£) <t]= Ate'*". (5.10) at With (5.10) and (5.9), the following equation is proved. se(f) = ht(t) = Aee-Xtt, (5.11) where At is the exponential intensity for the state. That is, in a continuous time Markov chain, the switching epoch probability density function st(t) and holding time density function h((t) for state £ are the same. This is a variation of Poisson arrivals see time average (Kao, 1997), but presented in a form that is convenient for the later development. The exponential distribution is the only distribution that has this characteristic. When IAT(f) of (5.6) is the holding time r, the mean holding time fe for mode £, is A~ex as shown below,r •; , The(T)dT= \ TAee-XeTdr = A-(x, o Jo which is equal to the mean switching epoch ft. However, on each mode £ segment with a switch-ing epoch T, the mean holding time is 0.5T. Therefor, it seems unlikely that the mean holding time equals to the mean switching epoch since the mean holding time for every segment is one half of the switching epoch. This paradox is explained that holding times from longer segments contribute proportionally more to the mean holding time than those from shorter segments while all the switching epochs contribute equally to the mean switching epoch. The memoryless property of the exponential distribution is fundamental to this paradox including the well-known "waiting time paradox" (for instance, see Feller (1971) page 11). 5.4 State Distribution The probability distribution of the continuous state x(t) is an important issue in jump systems because most likely it is not Gaussian and the variance of the state may not provide enough in-5.4. State Distribution 60 formation for utilities of control (Sworder and Choi, 1976). In this section, a system of integral equations for the state probability distributions is derived. The new method is not limited to linear systems, and is applicable for non-linear and state-constrained systems. We consider jump systems under state feedback control (it is assumed that both x(t) and £(t) are measured without noise), u(t) = fu(x(t),at)). (5.12) Then (5.1) becomes jx(f)=f(x(t),fu(X(t),at)),at))= fc{x(t),m\ (5.13) Let X c RNx denote the region where x(t) is defined, and S denote the Borel field on X. It is assumed that x{t) has a unique bounded stationary distribution LI(B) := F[x e B] and that the system function / and the feedback function fu are such that the closed-loop system function fc satisfies a Lipschitz condition in X and that the closed-loop system is stable in each mode. Therefore, the closed-loop system (5.13) has a unique solution g for a given initial condition xo at each mode segment where £(t) is constant. X{(XQ, r) = g(x0, r, £(T) = €) = ge(x0, r), (5.14) where r = t - ts is the holding time and X{{XQ,T) denotes x{r)\^T)={TX(Q)=XQ. Since fc satisfies a Lipschitz condition, ge{x$,T) is a continuous function of xo and T. For a fixed xo, let ge,X0(T) := gt(xo, T), which is a continuous function of r. Let B e & and xo be given, and T{(x0, B) be a set of T e [0, co] such that giXO(T) e B, that is, ,\.{xi Te(x0,B)=g^X0(B). (5.15) Since ge,Xo(T) is a continuous function of r, T((xo,B) is measurable. Let K{(xo,B) denote the probability that X{(x0, r) e B for a fixed initial value x 0. With the holding time density h(, Jr»oo \TeM{r)hiT)dT. (5.16) o The conditional probability distribution L>{(B) := P[x(f) e B\£{t) = £] is obtained by HAB)= f iueo(dx0)K((xo,B), (5.17) Jx 5.4. State Distribution 61 where //«>(•) is the probability distribution of the initial value xo for mode I. Now, we consider the probability distribution of x{t) at jump epochs tn. Let t~ denote the epoch just before t„ and t+ just after t„. The sampled values of x(t) at t~, where = I, are denoted by xe\(i), i = 1,2,..., that is, xt\{i) is the value of x(t) at the end of time segment where = I. Let Afi (JCO, B) denote the probability that X((xo, T) e B for a fixed initial value JCO and a switching epoch 7\ Since xt\(i) = ge(xo, T), a set of switching epochs that xe\(J) e B is T((xo,B) of (5.15). Therefore, with the switching epoch density Se(T), J f^OO \T(M(T)Sl(T)dT. (5.18) o Since st(t) = ht(t) by Equation 5.11, from (5.16) and (5.18), Kt(x0, B) = Kn(x0, B). The probabil-ity distribution of xe\(i), /*t\(B) is fia(B)= f Li(o(.dx0)Ka(xo,B)= f fim(dx0)Ke(x0,B). (5.19) Jx Jx Comparing (5.17) and (5.19), we conclude that ne(B) = Lie\(B). Theorem 5.1 For a stable Markov jump system, the conditional probability distribution of x(f) given £(t) = £, Li((-) and the probability distribution of x(t) at the end of mode I segments, pfi(-) are the same. By the above theorem, statistical characteristics of x(f) of a continuous-time Markov jump system can be obtained by investigating a discrete-time system (sampled at t~). Since Theorem 5.1 relates a continuous-time Markov process to a discrete-time Markov process, it may appear that there are some common aspects between Theorem 5.1 and uniformization. The basic idea of the uni-formization of a continuous-time Markov process £(t) is that under certain conditions £(f) can be represented as the composition of a discrete-time Markov process Z„, n > 0, and a Poisson process {N(f)}, i.e. £(t) = ZN(I) (Liu and Tu, 1999). However, the state space of a Markov process must be discrete for uniformization is applicable, whereas x(t) of Theorem 5.1 is continuous. Uniformiza-tion is a much stronger concept than Theorem 5.1 as it is valid not only for stationary distributions but also for transient analysis. Theorem 5.1, is only applicable for stationary distributions. In summary, Theorem 5.1 is independent from uniformization. In the following, a system of integral equations for the conditional probability distributions Ht, I e S is derived. Suppose that a state transition j -> I takes place. In this case, the initial 5.4. State Distribution 62 value xo for mode t is the value of x(t) at the end of mode j segment, which has the probability distribution Lij\ = pj. Let rjt be the probability that the prior state is j given that a state transition to I happened. In the thesis, is called the entry probability to state t, which satisfy, 2Zj&s,*e rji - 1 • Then, m = £ rjefij (5.20) Substitution of (5.20) into (5.17) yields jeS,*e J x (5.21) which is an Nf set of integral equations. Since it is assumed that x(t) has a unique bounded sta-tionary distribution p, when the stationary distribution of £(f), n is such that 7r, > 0, i e S, a unique bounded pe exists for all I e S. In principle, we can obtain pt, I e S (at least numerically) by (5.21) because lQ(xo, B) can be calculated from fc and Q for £(f). The task is much easier when p( and Kg have the densities (pe and Q>e respectively, that is, Pe(B)= I fa{x)dx, and Ke(x0,B)= I ®e(xo,x)dx. (5.22) JB JB Substitution of (5.22) into (5.21) yields, I <pt(x)dx= Z rjt I <f>j(xo) I ®e(xo,x)dx JB j€S,±e J x [JB = | | ^(x 0 ,x) 11 rje(pj{xo)dxo JB [ J x jcs,±e o dx. (5.23) Equation (5.23) must hold for all B e S, thus, the following system of integral equations for the densities is obtained. 4>e(x) = I ®e(xo,x) 2 rje(pj(x0)dx0, t e S . J x j€S,*e (5.24) Equation (5.24) is a system of integral equations for unknown functions <pe(x), € e S. Except for very simple cases, (5.24) must be solved numerically. However, this approach has two advan-tages compared to simulations: (a) no approximation due to time discretization is involved, and (b) randomness due to computer generated random numbers does not exist. Methods for solving (5.21) or (5.24) numerically are an important practical issue, but are not treated in the thesis. Note 5.4. State Distribution 63 that O^(x,xo) can be calculated from the knowledge of (5.14) and the distribution of the holding time T for state I. The overall probability density for x(r), <p(x) is obtained by taking the expectation with respect to £, X i i 0(*) = ZjiW<(*)- (5-25) i I , _ AO-Note that the probability density for the sampled state x(t~) is different from <p(x) as it is pt(pt(x), where pc is the £-\h element of the stationary distribution of the embedded Markov chain of £(f). 5.4.1 Entry Probabilities Entry probabilities are all dependent on state transitions, and are not influenced by how long £(f) stays in a state. A new Markov chain is constructed from the events of state transitions. This new Markov chain Zn is a discrete-time chain. Z„:={(tn), «=1,2 , . . . , (5.26) where tn is the «-th jump epoch. Z„ is the embedded Markov chain of g(t) mentioned in Section 5.2. Z„ is governed by the transition probability matrix P = [py] in (5.5). If y(t) has the stationary distribution n, Z„ has the stationary distribution p = \p\ pi •• • PN(] which satisfies pP = p and is related to n by the following equation (Tijnis, 1994), * = v ^ T T T ' ( 5 ' 2 7 ) From the definition of the conditional probability, _ P[Z„ switches from j to £] V j t = P[Z„ switches to £] ' ( ' ' As P[Z„ switches from j to £] = PjPjt, rH=J>IPJE , jeS,±t. (5.29) ItwPiPit 5.4.2 Examples Filtered Symmetrical Binary Random Signal As an example of the integral equation (5.24), filtering of a random symmetrical binary signal d(t) = ±b (b > 0) is considered (Figure 5.1). This is a rare case where analytical solutions are available. 5.4. State Distribution 64 b-i r - i 0--b d(t) a s + a H ) Figure 5.1: Symmetrical binary signal d{i) is filtered by a low-pass filter. The signal d(t) is a function of two-state continuous-time Markov chain £(r) such that d(t) = b when £(t) = 1 and d{i) = -b when r(t) = 2. With the following infinitesimal generator for r(t) Q = -A A A -A the stationary distribution n = [0.5 0.5], and d(t) is symmetrical, that is, the probabilistic proper-ties for d(t) = b and d(t) - -b are the same. Let x(t) denote the output of the filter. Since the filter is first-order, -b < x(t) < b. Suppose that £(t) just switched to state 1 (d(f) = b) at time to, then with T = t - to and xo = x(to), X(T) = (xo - b)e~aT + b (5.30) Let xi be a constant such that b > x\ > xo. From (5.30), X(T) = x\ when - i i X l ~ b T — T\ — —(X log -. x0 - b Since X(T) is monotone increasing when d(t) = b, x(r) < x\ for 0 < r < T\. AS T distributes exponentially with intensity A, W T ) ^ ] = 1 - . - n = 1 .gz|) (5.31) With B = Aa~\ the density of X(T) for state: 1, Oi(x0,x) is obtained by differentiating (5.31) with respect to xi. Oi(xo,x) = i 0 for x < xo Similarly the density O2(xo, x) for state 2 (d(t) = -b) is obtained 0 for x > xo O2(x0,x) = I 5.4. State Distribution 65 As £(r) has only two states, the entry probabilities r \ 2 = r 2 \ = 1. So (5.24) becomes (pi(x) = cp2(x) = J ®i(x0,x)(f>2(xo)dxo B(b -xf-x J (b - x0yp<p2(xo)dx0 ^ O2(x0, x)<p\(xo)dx0 B(b + xf~l £ (b + x0TpMxo)dx0 (5.32) (5.33) From the symmetry of d(t), (p\{x) = tp2(-x). Indeed, substitution of 0i(x) = <pi{-x) into (5.32) produces (5.33). So the system of integral equations (5.32) and (5.33) is equivalent to a single integral equation <px(x) = B(b - xf~x J (b - x 0 ) - V i ( -^o )^ 0 , which has a solution 0i(x) = Cl(b + xf(b - xf~x, (5.34) where c\ is determined so that £b<p\(x)dx = 1. The overall density of x(t) from (5.25), (p(x) = O.50i (x) + O.50i (-x) is «x)=Cib(b2 - x*t1=ci^-1 [ i - (5.35) When 8=1, that is A = or, (p{x) = c\b =constant. So x(t) distributes equally between -b and b. M \ • /• \ \ \ • -0.8 -0.6 -0.4 -0-2 0 0.2 0.4 06 OS -OA -0.6 -0-4 -02 0 02 0.4 0.6 0.8 Figure 5.2: Probability density <p(x) of filtered binary signal for B = 0.5 (at left), and B = 1 and f3 = 2 (at right). When 6 < 1, 0(x) is concave with the minimum at x = 0 and increasing towards x = b and x = -b 5.4. State Distribution 66 as shown at the left graph of Figure 5.2. On the other hand, when B » 1, approaches the Gaussian distribution as shown below. <Kx) = C l ^ - > ( i - ^ + ^ ^ - - ) Compare the above to the Gaussian density, •VSr 2 ^ 2 4 ^ / When B » \ , ( B - 1 ) 0 8 - 2 ) « 08-1)2 and so on, thus 0(x) « ^ e " * 2 ^ , where 2cr2 = b 2 / ( B - 1 ) . Figure 5.3 shows that (p(x) for B = 10 is already very close to the Gaussian distribution. For B > 1 of moderate values, 0(JC) assumes a shape between the flat distribution and the Gaussian distribution as shown in the right graph of Figure 5.2. The signal d{i) has a first-order low-pass spectrum with Figure 5.3: <f>(x) (solid line) for 8 = 10 and b = 1. The corresponding Gaussian density is shown in dashed line. the cut-off frequency 2A. So the filtered signal x(t) has a second-order low-pass spectrum with cut-off frequencies at 2A and a . Thus, it is possible to accommodate both the spectrum characteristics and the amplitude distribution by adjusting A and a . The filtered signal x(t) may be useful as a noise or disturbance model because it is bounded (the Gaussian noise is unbounded theoretically) and can have a wide range of distribution characteristics with a low-pass frequency spectrum. 5.4. State Distribution 67 Filtered Three-state Random Signal As an example of numerical solutions of the integral equation (5.24), a case of filtering a three-state random signal d(i) is considered. In this example, £(/) has three states, and d(t) takes one of three values depending on the state of £(t), that '-is;- d(i) = d(£(i)), and d(\) = dx = -\2A, d(2) = d2 = -0.24, d(3) = d3 = 0.76. The Markov chain £(/) has the following infinitesimal generator Q and transition probability matrix P : Q = P = -M qi2 qi3 q2\ -A2 q23 q-s\ qn ~h 0 pn pn P2\ 0 /?23 P31 P32 0 -2.0 1.0 1.0 0.5 -1.5 1.0 0.8 0.2 -1.0 0 0.5 0.5 0.333 0 0.666 0.8 0.2 0 The stationary distribution n of£(t) and the stationary distribution p of the embedded Markov chain Z„ are, n = [0.26 0.24 0.50], p = [0.3768 0.2609 0.3623]. The entry probabilities ry in (5.29) are r2l = 0.2308, r,i = 0.7692, rl2 = 0.7222, r 3 2 = 0.2778, n 3 = 0.5200, r 2 3 = 0.4800. Since the filter is a unity gain first-order low-pass filter, x(t) ranges between d\ = -1.24 and d3 = 0.76. Let xmin = -1.24 and xmax = 0.76. The continuous state space [xmin,xmax] is discretized for numerical methods. The grid size for discretization Ax is set to (xmax-xmin)/N, where an integer TV is a parameter to be chosen. Let xt,i = 1,2,..., N+l, be the /-th grid point, i.e. x, = xmin+(/-l )Ax. Let fa be an N + 1 vector that approximates fa, and fa(i) the /-th element of fa. Then, fa(x) « fa(i) forx/ - 0.5Ax < x < x, + 0.5Ax. (5.36) 5.4. State Distribution 68 Then from (5.24), for i = 1,... N+1 fa(i) = I Q>t(xo,Xi) rke(f>k(x0)dx0 (5-37) »^<be(Xj,Xi) ^ rU(j)k{j)&x. j=\ kzs,±e Let be an N+l x N+l matrix whose iy'-th element is 0^(x/, x,)Ax. The equation (5.37) is written as fa = Vt'Lkttn&k, for £eS. (5.38) To calculate we fix X y and evaluate O f(x7,x,) for the varying x ;. From (5.30), since the initial value is Xj, for mode I and sojourn time r,t . x(r) = e-aT(Xj - de) + de. (5.39) Let T(X) be r such that x(r) = x, then T(X) =-or 1 log (5.40) xj - de Since T is exponentially distributed with intensity Xe, the //'-th element of ^ is calculated as ( ¥ , ) , , = P[xi -0.5Ax < X(T) < xi +0.5Ax] _ e - A { T ( x , - 0 . 5 A x ) _ g -A f T(Xft0 . 5A ; r ) (x¥e)ij = 0 when (5.40) is negative or not defined. In this example, an approach that does not use the Jacobian is shown because the same principle is applicable when x(t) e R^*, Nx > 2. With ^e known, (5.38) can be solved for fa under the constraint that 'Zifa(i)Ax = 1. Figure 5.4 shows fa, fa and fa separately for N = 300. In this example, the results are stable for a wide range of N's. Figure 5.5 compares <f>(x) = 2Zentfa(x) calculated with N = 300 and N = 20. Except for the narrow peak at the center, two results are very close with each other. Figure 5.6 shows differences between <p{x) calculated with N = 300 and the histogram from a simulation run with 10,000 switches. This example shows the robustness of numerical solutions of the integral equation (5.24) even with a very simple numerical integration method. It is conceivable that the lack of time-discretization errors of (5.24) contributes to this robustness. 5.4. State Distribution 69 Figure 5.5: Comparison of the state density function <p(x) for two state discretizations. The one with a sharp peak is the result from N = 300, and the other from N = 20. 5.4. State Distribution 70 1 1 ' _ - J 1 1 -1.5 -1 -0.5 0 0.5 1 xffl Figure 5.6: The difference between <p(x) for N = 300 and the histogram from a simulation run of 10,000 switches. 5.5. State Moments of Linear Jump Systems 71 5.5 State Moments of Linear Jump Systems The state moments of linear jump systems can be obtained without explicitly calculating the state probability distribution. As a simple example, the first moment (mean) is treated in this section. A linear jump system is described by the following differential equation. ^x(o=A(mm+B(t(t)Mt)+dm) d t (5.42) = Atx(f) + Beu(t) + de, ((f) = t, where x e RNx is the state , and u e RNu the input. A, B and d are matrices of appropriate dimen-sions. A mode-dependent constant disturbance de is included to accommodate "load change " type disturbances common in process control applications. The state may have a discontinuity (jump) at jump epochs (Sworder and Rogers, 1983). As the state discontinuity can be handled easily, it is assumed that the state x(t) is continuous to simplify the presentation. It is assumed that x(t) and ((f) are measured, and linear state feedback control with a bias is applied (Sworder and Rogers, 1983), u(t)= -Ge(x(t) + ae), (5.43) where the linear gain Ge and bias cte are mode dependent. Then, the closed-loop system becomes d —x(t) = Ac{x(t) + d(, (5.44) where A°( := Ae - B{G{ is the closed-loop system matrix and de := de - BeGeote the closed-loop constant disturbance. The system (5.44) is a hybrid system whose state is composed of a continuous part x(t) and a discrete part ((f). It is assumed that the system (5.44) is stable in each mode and de is finite for all (,. Let me denote the conditional mean of x(t) given that ((f) = £, that is, me = fx X(f>e(x)dx. The overall mean of x(f), m is obtained from (5.25) as follows. m = I xcf>(x)dx - Y?elxne I x(pe(x)dx = Y^ixnemt. (5.45) Jx Jx Let xe\ denote the value of x(f) just before ((f) leaves state I. By Theorem 5.1, me is the mean value of xe\. With a switching epoch T and the initial value xo := x(ts), xei AC<Tx0+ I u< • Jo = eA T 0  eAC(sds-d(. (5.46) 5.5. State Moments of Linear Jump Systems 72 Since xe\ is a function of T and xo, the conditional expectation of xe\ for a given initial value xo, ni((xo) is calculated with the switching epoch density function S((T) = Aee~X{T, m((xo) = Qoxo + C(\d(, where Qo and Q i are mode-dependent constant matrices, e W e - W d T , C n : = A e \ \ e t f ' d s e X t T d T . o Jo Jo Then me is obtained taking the expectation with respect to xo m t = E x o [me(x 0 )] = Q 0E(x 0) + C t \d t = Ceoxeo + C t \de, (5.47) (5.48) where xeo is the mean value of x(ts)\ t^t)=(, that is, the mean value of x(r) when £(t) just switches to state I. When the prior state is j, XQ = Xj\, so with the entry probabilities ry/, xeo = H nMxj\] = 2 rjtntj. Substitution of (5.49) into (5.48) yields ™e = E rjeCeomj + Ceide, £&S. (5.49) (5.50) So me,t e S is obtained as the solution of a system of linear equations constructed from (5.50). When the overall mean m is finite and the stationary distribution of £(?), n is such that 7T, > 0, / e S, the conditional mean me is finite for all t e S, and Equation (5.50) has a unique solution. The overall mean m is obtained from n and /w '^s as shown in (5.45). As an example, the linear equation for = 3 is shown below. I -r 2i Cio -r*\C\o \ -rnCio I -7*32^20 -rnCzo I m\ m2 W3 C n d i Cndi (5.51) 5.5. State Moments of Linear Jump Systems 73 The mode dependent second moment E[(x(t)x(t)T\((f) = i] can be obtained by the same approach with more complicated algebra, and detailed in the Appendix of this chapter. With the knowledge of me and the mode dependent second moment of x(f), the quadratic performance index J can be evaluated, J := E[x(t)TRlX(t)] + E [ M (O r W)J = IMk + IMk- (5.52) t i With the new method, ||x||«, and \\U\\R2 are evaluated separately. This is advantageous because R2 is used as the Lagrange multiplier to limit the control input energy. By evaluating ||x||^, and \\U\\R2 for different R2, the cost of control input for the control performance in terms of the weighted state variance is known. This is impossible by evaluating just J. Another advantage of the new method is that the mode dependent ||x|k and \\U\\R2 can be evaluated. The optimum control theory for linear jump systems allows R\ and R2 to be changed by mode (Sworder, 1976). Evaluating \\x\\R] and \\U\\R2 for each mode allows R\ and R2 to be adjusted to obtain desired compromises between control input energy and control performance. 5.5.1 Example The filtered 3-state signal of Section 5.4.2 is studied to compare the variance calculation of the new method and the one based on the autocovariance. The cut-off frequency a is set to 2. The mean input is zero as Y,entde = 0. From (5.30),, , X(T) = e~aTx0 + (1 - ear)dt. Thus, aXt_ Ceo= f" eatXeeXttdt= — Jo <* + C « = C (l-e-at)A(e-Ae'dt= — Jo a + Xe The system of linear equations (5.51) becomes, 1.0 -0.3095 -0.1733 1.0 -0.1600 -0.3846 m\ -0.6200 -0.1190 m2 - -0.1371 1.0 m3 0.5067 (5.53) (5.54) (5.55) t 5.5. State Moments of Linear Jump Systems 74 which yields mx = -0.5024, m2 = -0.2474 and m3 = 0.3800. Entities for the variance calculation (see Appendix of this chapter for notation): G ® G = a 2 e - 2 a t GGt G ® H = H ® G = a ( e a t - e~2at) = > H G t = G H t = a 2a + Ae' — — I Xe Ae a + Ae 2a + Ae/ H®H = {\- 2e~at + e~2at) => HHe = 11 2Ae + Ae Vl 0.4153 v 2 = 0.0432 v3 0.2310 a + Ae 2a + Ae} Let \e denote the conditional second moment. The system of linear equations (5.76) becomes, 1 -0.0769 , -0.2564 -0.1970 1.0 -0.0758 -0.1040 -0.0960 - 1 which yields vi = 0.5048, v2 = 0.1653 and V3 = 0.2994. Since m = 0, the variance T is equal to the second moment v. Thus T = v = neve = 0.3206. (5.56) t Another approach for the evaluation of the variance is to calculate the autocovariance yx{i) of x{t). Then the variance is yx(0). As a first step, the autocovariance of d(t) is calculated from the probability row vector p{i) = [pi(t) p2(t) Pi(t)], where pe(t) is the probability that £{t) = I. The following differential equation governs the evolution of p{t). -p{t) = p(t)Q, (5.57) which is solved under a constraint that 2ZePe(t) = l.Vr > 0. The conditional autocovariance •ydte(t) = E[d(0)d(t)\d(0) = de] is obtained by solving (5.57) with the corresponding initial condi-tions. For instance, let q(i) = [q\(t) #2(0"'#3(0] be the solution of (5.57) with the initial condition [1 0 0]. Then jd,\{t) is obtained as 7rf.i (0 = q\(t)d\ + q2(t)dxd2 + qi(t)dxd3. (5.58) The (unconditional) autocovariance jd(i) is obtained by taking the expectation on d(0), that is 7 d ( t ) = J ] n e 7 d A t ) - (5-59) 5.6. Linear Quadratic Optimal Control 75 In our case, we get yd(i) = c x e a t + c 2 e b t for t > 0, (5.60) where a = 2, b = 2.5, ci = 0.1520, c2 = 0.5504. 2 Since x(t) is the output of GF(S) = - , yx{t) is calculated from the weighting function w(t) = s + 2 . ( 2e-2tofGF(s): yx(f) = JT wd ){ j^° w(fe)yrf(f + f i - fe )*2J<&i. (5.61) Noting that the autocovariance is an even function of time lag, we get 7 x ( 0 ) = 4 ^ ° <T2" | ^ ° ° e- 2 ' 2 [c,e- 2 1"- ' 2 1 + c 2 e ~ 2 5 ^ } d t 2 ^ dtx Ac\ 4c? = -7T + -7T = 0-3206 8 9 (5.62) , which is identical to (5.56). 5.6 Linear Quadratic Optimal Control The subject of this section is the optimal control of jump linear systems with mode dependent constant disturbance, which was introduced in Section 5.5. The system equation and definitions are repeated here. x(t) = A(x(t) + Beu(t) + d(, (5.63) where x e R N x is the state , and u e R^" the input. A t , B ( and d ( are mode dependent matrices of appropriate dimensions. The performance index is the long term average of quadratic functions of the state and input: J := lim | JT {x(t)TQtx(i) + u{t)TReu{ij)dt^, (5.64) where R( is a positive-definite, and Qt semi-positive-definite matrices. Note that Qt and Re may be mode dependent. The state may have a discontinuity at mode changes x(ts) = x(t~) + Sn, (5.65) 5.6. Linear Quadratic Optimal Control 76 where switches from / to j at epoch ts, and dy is a state jump from mode i to j. In Sworder and Rogers (1983), the problem is solved for de = 0. The optimum control u*(t) is a state feedback (mode dependent) with bias: u*(t) = -R-txBTePe(x{i) + ae), when = € (5.66) where Pe is the solution of the following coupled Ricatti equation: PeAe + A\Pt - PeBtR?BTtPt + ZqejPj + Qe = 0, for t e S, (5.67) j where qej is the entries in the infinitesimal matrix Q of £(r). The bias term ae is obtained by solving the following system of linear equations: (At + P? Qe)ae + Y.qtjP~elfj(at - a} - 6t]) = 0, for I e S. (5.68) j In Khanbaghi et al. (2002), the constant disturbance de is absorbed in the state discontinuity: x(t) = Aex(t) + Btu(t) + d{=A{ [x(t) + A~( 1 de] + Beu(t). (5.69) With a new state, z(t):=x(t) + A^de, (5.70) the system equations becomes, z(t) = A(z(t) + Beu(t), (5.71) where the state jumps dy includes A~(xde. Then the method of Sworder and Rogers (1983) can be applied. However this approach makes it difficult to express x(t)TQex(t) in the performance index in terms of the new state z(t). In Khanbaghi et al. (2002), x(t) = \y(t) u(t)]T and A]1 did not affect y(t) part and Qe matrix has only non-zeros element for y(t), so this problem did not surface. If it is ^' :'• desired to include u(t) term in the performance index, this approach does not work. The problem was solved with the original state x(t) and non-zeros disturbance de. The detailed mathematical derivation is included in Appendix of this chapter. The coupled Ricatti equation is the same as (5.67), and the equations for ae become as follows: (Ae + Pe-lQt)at-dt + 2Zq*jPt'lPAat-aj-6tj) = 0, forteS. (5.72) 5.7. Summary 77 5.7 Summary It is shown that continuous-time Markov jump systems can be investigated as discrete-time systems sampled at the end of each mode. From this fact, two practically important applications are derived: calculations of state probability distributions of jump systems, and calculation of the state moments of linear jump systems. The probability density function of a filtered binary random signal is obtained analytically. Unlike approaches'based on simulations, the new method does not have time-discretization approximations nor randomness due to computer generated random numbers. 5.8. Appendix - Mathematical Details 78 5.8 Appendix - Mathematical Details 5.8.1 Covariance of x(t) In this section, the subscript I for mode is omitted if no confusion exists. Since T := E[(x(r) - m){x(f) - m)T] = E[x(t)x(t)T] - mmT, and m is known, the second moment v := E[x(t)x(t)T] is calculated. We write (5.46) as where G(f) := e*1 and H(t) := j f e*'ds. x(t)x(t)T = G(t)xQxT0G(t)T +H(f)ddTH(t)T + H(i)dxlG(f)T + Git)x0dT H(t)T. Taking vec operations on both sides (dropping t argument for G and H) noting that vec(ABC) = (CT ®A)vec(B), vec (x(t)x(tf) = (G <8> G)vec(x0xT0) + (H® / / ) v e c ( J < f ) + (G® H)vec(dxT0) + (H® G)wec(x03T). For a fixed xo and mode I, vQc{$[x{t)x{t)T\x0]) = E { G ® G ) v Q c ^ The following expectations are calculated. vec(E[x(r)x(Or|*o]) = GGvec(x0xS) + HHvec(ddT) + GHvec(dxT0) + HGvec(x0dT). The conditional second moment of x(t), is obtained by taking the expectations on xo, vec(v )^ = EXo(vec(Et[x(t)x(t)T\xo])) = GGvec(EX0(x0xT0)) + GHVGC(3EXO(XT0)) + 5Gvec(EX o(x 0)^) + HHwecidd7) = GGvec(Vfl)) + GHvec(dxT0)) + HGvec(x0dT) + IMvecidd7), dt. HH, HG and GH are calculated similarly. So 5.8. Appendix - Mathematical Details 79 where veo := E ( X 0 X Q ) and x0 : = E (x 0 ) . It is shown in Section 5.5 that xeo = Itjes^jemj, (5.73) where it is understood that rjj = 0. If the prior state is j, XQ = J C 7 I S O *o*oV=./ = Xj\xTjX By the same reasoning as for the evaluation of the mean of x(t), E(xnxTjX) = Vj. So jeS,*e where is the entry probability from j to € (5.29). Finally vec(v )^ = GG^vec (Eyastyv,-) + GHevec (d E / G S O * W ; ) + HGevec (iZjesn^A) + HHfvec(d{d]). The above is rewritten as vec(v )^ = GG^vec ( Z ^ / V , ) + ct, (5.74) where ce is a vector made by combining terms without v, and known after the calculation of AM,-. ce = GHevec (d Zjesntrf) + ^ v e c (ZjesneMj^e) + HHe^c(d{dTe). (5.75) A system of linear equations is constructed for \ t . An example for = 3 is shown below. I -r2iGG\ -mGGi -rX2GG2 I -r32GG2 -rl3GG3 -r23GG3 I vec(vi) vec(v3) c2 C3 (5.76) 5.8. Appendix - Mathematical Details 80 5.8.2 Linear Quadratic Control The following jump system is considered. -X(t)=Rx, u, t, m=Mx, u, t), at) = t (5.77) Subscripts I denote the state of a continuous time Markov process (it) with the infinitesimal gen-erator Q = [qij]. The problem is first solved for finite horizon, and later steady solution is sought. The performance index J is Jr'f 1 T r'/ Z(x, u, t, at))dt = E L((t)ix, u, t)dt o J LJo (5.78) Note that the penalty function L may be a function of £(0, and its dependence on at) is denoted with the subscript on L. The to-go value function F(x, t) when at) = ^ is denoted with Vtix, t): Vtix, t) - min E I LAx, u, t)dtUt)=e % Jt (5.79) Also in our case, there are discontinuity in the state when at) changes its state. Let (5y denote the deterministic jump of the state when at) shifts from / to j. xit) = xiQ + 6u iat) = j\at-) = i) (5.80) The following coupled Hamilton-Jacobi equations for I e S are obtained following Wonham (1970) and Sworder and Rogers (1983): —rVtixj) = min at u dVe L(ix, u, t) + -r-ftix, u, t) ox + %qtj[Vj(x + 6tj,f)\, i6ee = 0) (5.81) Linear System with Bias x(0 = A(xit) + Btuii) + dt dt U = xitf&xit) + u(t)TReu(t) It is assumed that Vt has the following form:. (5.82) (5.83) Vtix, t) = ix + atit))TPtit)ix + a tit)) +Be(t). (5.84) 5.8. Appendix - Mathematical Details 81 Differentiation of the right-hand side of (5.84), (time argument t is dropped to streamline notation) with respect to t and x yields, and So LHS = -—V{ = -xTPex - d/Pedt - aTePta( - aTePeat ot - (aT{P{ + aT(Pt)x - x\Peae + Ptat) -0e, ^-Ve= —(x + ae)TPe(x+:ae) = [Pe(x + ae)Y + (x + ae)TPt. ox ox a- vtft = Vt{Atx + B(u + dt) ox ox = [Pe(x + at)Y{Atx + B(u + de) + (x + ae)TP{(A{x + Btu + di) = (Atx + Beu + dtfPeix + ae) + (x + ai)TPt(Aex + Btu + di) = uTBTePt(x + ae) + (x + at)TPtBtu + (xTAT( + dTe)Pt(x + a() + (x + ae)TP((Aex + dt). (5.85) Then U + s-Vefe = xTQtx + uTReu + uTBTtPt{x + at) + (x + ai)TPtBtu ox + (xTA\ + d^Ptix + ae) + (x + ai)TPt(Atx + de), which takes its minimum at u* = -Rf BTePt(x + at). And m i n u Le+^-Veft ox = -ix + ae)TP(BeR^BTePe(x + at) + xTQtx + (xTAT( + dfiPtix + at) + (x + at)TPt(4tx + de). A matrix St is defined as, Se '•= P(B(RexBTtP(. (5.86) (5.87) (5.88) 5.8. Appendix - Mathematical Details 82 Substitution of (5.87) and (5.84) into the RHS of (5.81) yields, RHS = -(x + at)TS{{x + a{) + xTQex + (xTATe + dT{)Pt(x + at) + (x + aATP((Aex + dt) + %qtj [(* + (*j + Stj)TPj(x + aj + 6tj) + Bj] = x PtAt + A\Pt - St + ZqejPj + Qt x - xTSiat - aTtStx - aTeStat + xTAT{P(at + d(Pe(x + at) + aTePe(Atx + dt) + xTPtdt + %qtj [xTPj((*j + 6(j) + (aj + StjfPjX + (aj + 5tj)TPj(aj + 6tj) + Bj] = x PtAt + A\Pt -St + zZqtjPj + Qt x + xT[-Stat + A\Ptat + Ptdt + }ZqtjPj(aj + 5ej)] + [-aTtSe + aTePeAt + dTePe + zZq(j(aj + 6tj)TPj]x j - aTeSeat + dTtPtat + aTtPtdi + Zqej [(aj + Stj)TPj(aj + 6{j) +Bj] JliiS (5-89) From the Hamilton-Jacobi equation (5.81), LHS (5.85) = RHS (5.89). This must hold for any x, so equating xT [ • ]x term in the both sides, the following coupled Ricatti differential equation is obtained. -Pi = PtAt + A\Pt -St + ZqejPj + Qt | (5.90) Next, from xT[ • ] and [ • ]x terms, -Ptdt - Peat = -Stat + ATtPtat + Ptdt + ZqtjPMj + hj) Pi is substituted with (5.90). -Peae + \ptAi + A\Pt -Se + zZqtjPj + Q)J ae = -Seae + ATtPtai + Ptde + ZqtjPjiaj + 5{j) So -Ptdt = -PtAtat - Qtai + P(dt + 2ZqtjPj(aj + 8tj - ae) From the above, the differential equation for at is obtained. at = (At + PilQi)at -de + zZqtjP~tXPj(at ~ aj - o~ej) (5.91) 5.8. Appendix - Mathematical Details 83 Equating the terms without x, -d/P e ae - a T ( P t a t - a T t P t de - f ie = -a T e S { a e + <f ( Peae + a T ( P ( d { + [(aj + 6ij)TPj(aj + <%) + Bj\ From the above, a differential equation for fie is obtained, of which derivation is skipped as fie is not used for control. Stationary Solution If the stationary solution exist, P i , de and fie are all zero. From P( = 0, we get the coupled algebraic Ricatti equation. * i ; 1 (5.92) PtAt + AT(Pt -S< + ZqejPj + Qe = 0 For ae, (At + P~txQt)at -de + iZqtjPf Pj(a{ - aj - 6TJ) = 0 Since ae - a t - See = 0, So ( A e + r { x Q e ) a e -d(+ 2ZqtjP?Pj(fXt - a } - Stj) = 0 (Ae + p-elQe+2ZqtjPe~lPj)(Xt-zZqtjPePjaj = de + XqejP? Pj5{J j*e j*e j*e (5.93) (5.94) (5.95) A system of linear equations for a\,..., a^( is constructed by stacking a e ' s . Ax+P?Qx ''(t :'.!*« +zZj*mP\xPj -qnP[XPi ~q^P-xXPz A 2 + P - 2 l Q 2 -qixP^Px +lZj*2q2jP^1 Pj -qi3P-2lP3 Ai+P-^Qi -qnP-3lPi -qiiPjPi +2ZjriqyP~3[Pj ai oc2 <*3 di + ZjtiqijP^PAj d2 + 2Zj^q2jP2-1PjS2j d3 + 2ZmqijP^Pfrj Chapter 6 Broke Storage Tank Level Control 6.1 Introduction Sheet breaks cause the biggest upsets in paper machine operation. When sheet breaks occur, the paper sheet cannot be rolled into a jumbo roll and is routed to the broke storage tank after it is repulped with the addition of water. In the thesis, this incoming flow to the broke storage tank is called "broke flow" and denoted by Fb(f). The broke storage tank is the first process to receive the upset, and therefore to manage the broke tank level properly goes a long way towards stabilizing the paper machine wet end. Presently, very few broke tanks are under automatic control. Effective automatic broke storage tank level control would improve the paper machine operation greatly. In this chapter, the nature of the input disturbance to the broke storage tank is investigated. Then two control scheme, one linear and the other nonlinear will be presented. Utilizing the results of the last chapter, the probability distribution of the tank level can be evaluated reliably, thus each control scheme can be assessed without simulations. 6.1.1 Paper Machine The paper machine produces paper continuously from a very thin pulp slurry (0.3 - 0.6% consis-tency) by first removing water by a rotating mesh of plastic fabric (wire), then pressing and finally drying with steam. Approximately 98% of the water is removed through the wire, and is called "white water". The pulp slurry contains, in addition to pulp fibres, fines and various additives, 84 6.1. Introduction 85 such as clay, dye, starch, and retention aides. Non-negligible amounts of fines and additives pass through the wire and get into the white water. It is important to put these lost fines and additives back into the paper sheet to maximize resource usage and minimize discharges to the environment. This is achieved by a system of disk filters, called a "saveall". Also, some of the white water is used for dilution in the pulping process, and thus eventually comes back to the paper machine. A very simplified material flow in the wet end is depicted in Figure 6.1. During normal operation, the paper machine produces some waste paper in the form of trim-mings, off-grade paper, and paper produced during sheet breaks, which cannot be incorporated into a paper reel. This kind of paper is called "broke", and is put back into the paper machine after it is reintroduced into the pulp slurry by adding water. ft WIRE J W.W. DRYER SAVEALL MACHINE CHEST Fresh Pulp BLEND CHEST BROKE TANK Figure 6.1: Simplified diagram of a paper machine wet end. Fines and short fibres in white water (W.W.) are recovered at the saveall, and mixed with fresh pulp and broke stock at the blend chest. 6.1.2 Benefits of Wet End Stabilization Modern paper machines produce finely graded papers from a complex mixture of various pulps and additives at high speed. Different kinds of pulps including ground wood pulp, chemi-thermal mechanical pulp, kraft pulp, and recycled (deinked) pulp are utilized to attain the desired paper quality economically. In addition to these materials, saveall stock and broke are added to minimize raw material usage as shown in Figure 6.1. These pulps have different physical and chemical characteristics, including average fibre length, strength and freeness. Dewatering and retention characteristics of the pulp slurry on the wire are 6.1. Introduction 86 in turn determined by the pulp slurry composition. For instance, kraft pulp, which is rich in long fibre, tends to be retained on the wire, and fines help to keep additives on the wire by clogging up spaces between pulp fibres. As the pulp slurry dewatering characteristics change, the amount of water carried over to the press and dryer sections changes. This causes moisture variations in the final paper until the moisture control system compensates, which takes longer than several minutes. Any composition change in the pulp slurry affects machine runnability and the quality of the final paper. At the same time, the amount and composition of the white water change. The white water is used in pulping processes and in the stock preparation area. Therefore, changes in white water composition come back to the pulp slurry through many routes with different time constants. Thus, a pulp slurry composition change starts a complex and long-lasting chain of events which affects the paper machine operation and the final paper quality. A number of studies on the white water network using dynamic computer simulation have been published ((Bussiere et al., 1987), (Croteau and Roche, 1987), (Kropholler et al, 1994), (Orccotoma et al, 19976) and (Orccotoma et al, 1997a)). In Orccotoma et al (1997a), fines distribution in a white-water network was investigated for a single line 600-tons/day newsprint milk Simulation results showed that the effects of a step change in the fines content in the fresh pulp lasted more than 5 hours. Stabilizing the wet end operation brings about many benefits (Gess and Kanitz, 1996). First it stabilizes the pulp slurry composition, which keeps the final paper quality constant. Also, it is known that stable pulp slurry composition reduces sheet breaks, thus increasing paper machine efficiency. This not only improves the economy of the paper making operation, but also minimizes raw material usage, energy consumption, and off-grade products, which otherwise have to be re-processed with additional energy. Variations of flows in the stock preparation area are minimized. This reduces overflowing and emptying of the storage tanks, which means reduced fresh water usage and reduced discharge of effluent. In Bonhivers et al. (19996), it is shown that coordinated control of a white water network may reduce fresh water usage and effluent discharge drastically. 6.1. Introduction 87 6.1.3 Broke Handling During normal operation, approximately 3% of production is routed to broke storage as trimmings. During a sheet break, 100% of production is routed to broke storage. Therefore, assuming the same consistency for trimmings and broke, the incoming flow to the broke storage increases 33 times. Thus, the broke storage tank level increases suddenly when a sheet break occurs. As sheet breaks happen at random, it is not easy to keep a constant broke ratio without overflowing or emptying the broke storage tank. When a sheet break happens, not only does the broke storage tank level increase, but also the white water demand increases to dilute broke at the repulper, and the saveall load increases due to the increased cloudy flow from the broke thickener. In short, sheet breaks initiate large upsets in the white water network. In Bussiere et al. (1987), the effect of sheet breaks on the fines content in the fresh pulp was studied by a computer simulation. They found that a 40 minute long break caused 2.4 to 5 % increase in the fines content, and it took 5 hours to stabilize. Usually 10 to 20% of the total production is broke. Broke has different characteristics from other furnish components since it has already been processed by the paper machine. The amount of broke in the furnish is defined by the broke ratio. It is important to keep the broke ratio constant to keep the paper quality constant. 6.1.4 Literature Review There are only a few studies published on automatic control of broke storage tank level. Khanbaghi et al. (2002) treated a network of a broke storage tank and a white water tank. Jump linear quadratic control was utilized because the broke flpw;;was modelled as a finite state Markov process. The flow roughness was measured by the variance of the outlet flow rate change, Var[zi(0] = Var[v(x)]. First passage times of overflow were obtained to assess the performance of the resulting controllers. This thesis owes many ideas, including the broke flow modelling and application of jump system theory, to their study. In Dabros et al. (2003), the optimal way to change u(f) was investigated. A notable feature of the study is that changes of the pulp properties at the head box obtained by simulations, is used as the performance index. However, the time and the final value of u(t) are assigned by the operator. Thus, this study gives some hints on automatic control, but is not enough to construct an automatic 6.2. Characteristics of Broke Flow 88 control scheme by itself. 6.2 Characteristics of Broke Flow The Poisson process, where waiting times distribute exponentially, is a convenient mathematical model for events occurring completely randomly (memory-less). If breaks start and end randomly (in a memory-less manner), broke durations and normal durations follow exponential distributions. Then the broke flow is modelled by an alternating Poisson process as shown in Figure 6.2. Both x2 Figure 6.2: Alternating Poisson process Xt and Yi distribute exponentially with different intensities. Khanbaghi et al. (1997) and Section 7.2 show that this hypothesis holds well for industrial data. Thus, sheet breaks are modelled as a two-state continuous-time Markov process'%($)' 1 normal operation 2 sheet break (6.1) (6.2) It is assumed that the normal operation duration T\ follows the exponential distribution with inten-sity X\ and the break duration T 2 with A2, that is, P[7i < t] = 1 -F [ T 2 < t ] = l - e~Xlt. Then the infinitesimal generator Q of £{t) is A2 —A2 A simple assumption is that the broke flow 'F&(f) is a deterministic function of £(r), Q = (6.3) Fb(t) = F„(at)) = \ Fbi #0=1 Fb2 #0 = 2, (6.4) 6.2. Characteristics of Broke Flow 89 where Fb\ and Fbi are constant. This assumption makes it possible to model Fb(t) as a two-state continuous time Markov process. From a simple material balance point of view, this assumption seems reasonable. Most of the published simulation studies on broke management assume that the broke flow is binary. It turns out that the industrial data used later on this thesis shows non-binary characteristics. However, the binary broke flow assumption is adopted in the thesis because (a) the powerful theory of Markov processes can be utilized, and (b) even though the industrial data is not quite binary, it can be approximated by a binary flow in such a way that the resulting controllers are conservative with respect to overflow probability. 6.2.1 Spectrum of Broke F low To utilize linear optimal control theory, the power spectral density of the broke flow is obtained. Let a row vector p(t) = [pi(t) P2(t)] denote the probability that £(t) is in state 1 or 2 at time t. It is known that pit) evolves according to the following differential equation: jtp{t) = p(t)Q. (6.5) Since p2{t) = 1 - pi(t), P\(t) = -Aipi(t) + A2p2{t) = -A\p\(t) + A2{\ - pi(t)) = -(Ax + A2)pi(t) + A2 (6.6) So p\(t) for the initial condition p\{0) = c is, p x ( t ) = c e - ( A l + A 2 ) t + A 2 f e - ^ + h ) T d r (6.7) Jo Let q\(t) denote the probability that £(t) = 1 when £(0) = 1. Substitution of c = 1 into (6.7) yields (6.8) "1" A2 / l l 1- A2 The following symbols are defined to streamline notations: A\ + A2 A\ + A2 Ai + A2 Ai n A2 Ab '.- Ai + A2, i\ :- —• A] Then, b := i + A2, li := —!—, l2 := (6.9) Ai + A2 Ai + A2 q i ( t ) = t x e X b t + i 2 (6.10) 6.2. Characteristics of Broke Flow 90 Autocovariance of Broke Flow The autocovariance of the broke flow Fb(f) is obtained using qi(t). The mean length of T\ is \/Xx and T2 1 /X2. The mean flow rate Fm is F m = ^ + ^ = i i F b x + { x F b 2 ( 6 - n ) A\ + A2 A\ + A2 Let fb(f) denote the zero-meaned Fb(t), and fb\ and fb2 denote fb(t) when £(t) = 1 and 2 respec-tively. /6(0 := Fb{t) - Fm fti := F M - Fm = (1 - fc)FW - £ , F M = £x(Fbx - Fb2) = W « - / « ) < 0 (6.12) ./M := F M - F M = (1 - £x)Fb2 - l2Fb\ = - F M ) = £2(fb2 - fb\) > 0, where an identity fbX- fb2 = FbX - Fb2 is used. Let y(t) denote the autocovariance function of Fb(t), 7(0 := E[(F,(0) - Fm)(Fb(t) - Fm)] = E(fb(0)fb(t)] (6.13) Here, y(t) is calculated as the expectation of the conditional expectation given £(0): 7(0 = Bm[E[fb(0)Mt)W)]] = mwmMtmo) = i]+mifb(o)Mt)\m = i\ (6.14) where 71 (0 and 72(0 are the conditional autocovariances. For instance, 7,(0 := E[fb(0)Mt)\a0) = 1] = fbiE[fb(t)\m = 1]. (6.15) When £(0) = l,fb(t) is fbx with the probability px{t) and fb2 with the probability 1 - px(t). Thus, 7i(0 =Pi(t)fb\ + 0 -Pi(t))fbifb2=pi(t)(fb2l-fbxfb2)+fbxfb2. (6.16) Substitution of (6.10) for px(f) yields, 7i(0 = (heXbt + t2)(f2bx - fbiM) + fbifbi = £xe-^(fb\ -,/*,,/„) + £2{flx - fbxfb2) + Mb2 (6.17) = e-Xbtfbxex<JbX-Jb2)+Mt2<Jbx - / M ) + / M ] - Jb\e ' 6.2. Characteristics of Broke Flow 91 where between the third and fourth lines, (6.12) is substituted. Similarly, y2(0 = fb22e Xb'- Substi-tution of j\(i) and y2(0 into (6.14) yields, 7(0 = tiyxif) + ^ 72(0 = {hf h x + hf 2 b l)e h t (6.18) As liflx + t\jl2 - E(Fb - fm)2, which is the variance of the broke flow, it is denoted with cr2,. cr2 := Var[Fb] = E[Fb - E(Fb)}2 (6.19) y ( t ) = oie-Xbt, 0 < t (6.20) Since the autocovariance functions are even functions of time lag, the above expression is extended to include negative t, which is necessary for spectrum calculations. y(0 = oie-m, -co < t < oo (6.21) Spectrum of Broke Flow The spectrum of the broke flow Sb is the Fourier transform of y(t). Sb(oj2)= y(t)e-J"'dt= f %J —CO KJ—OO oieXb'e-jwldt + £° oie-Abte-Jojtdt f O— oo J r°° 1 T 1 1 e-(h+Mtdt =(J2 + 0 AB - J'OJ AB + JOJ (6.22) 2oiAb 2oi A2H oo2 + A2B AB to2 + A\ The spectrum of the broke flow is that of a first-order low-pass filter with a cutoff frequency A B . For the linear optimal controller design, the broke flow (zero meaned fb(t)) can be treated as the output of a unity gain first-order low-pass filter G b(s) driven by white noise of the intensity cr 2. G„(s) = s + AB' 2o-2 white noise ——-AB s + AB broke flow fb(s) (6.23) 6.3. Flow Smoothness Requirements 92 6.3 Flow Smoothness Requirements What is the appropriate flow smoothness or roughness indicator for the outlet flow Fu(i) of the broke tank? The immediate downstream process is the consistency control loop of the broke stock flow to the blend chest. However, the effects of Fu(f) reach much further, and the final paper quality and runnability of the paper machine are two major concerns for the paper maker. For paper quality, it is best to keep Fu(t) constant, but for most of paper machines, the size of the broke storage tank and break frequency make this impossible. The second best would be to minimize changes in the broke ratio, which mathematically Var[Fu(r)] can represent nicely. For some downstream processes that are regulated by a feedback controller, for instance, broke stock consistency, Var[Fu(t)] = Var[v(f)] would be a better indicator as shown in Chapter 4. Some production people believe that sudden changes in the flow rate increase the probability of sheet breaks, although the wet end operators routinely make step changes in Fu(i). There has been no statistical study published on this subject. Since there are many causes of sheet breaks, it would be difficult to obtain clear-cut conclusions on correlations between the speed of flow rate changes and sheet breaks. However, to constrain Var[v(r)] would be beneficial for regulated downstream processes and acceptance of the control scheme in the field. For instance,.in Khanbaghi et al. (2002), Var[v(r)] is used as the flow roughness indicator. " , ; Thus, it would be best to constrain both Fu(t) and v(t) := w(r). In the thesis, a weighted sum of Var[«(r)] and Var[v(r)] is used as the flow roughness indicator for linear controllers, and max u(t) and max |v(r)| are constrained for nonlinear controllers. 6.4 Linear Optimal Controller Since the input disturbance (broke flow) has a first-order low-pass spectrum, the method described in Section 3.4 is applicable without modifications. For completeness, key definitions and equations are repeated here. Let Fm denote the mean value of Fb(t), Fm := E[Fb(f)]. The input disturbance is zero-meaned, that is, Fm is subtracted from Fb(t). The state variable x(t) is defined as follows. x(t) :=[x,(0 x2(t) u(t)]T, (6.24) where 6.4. Linear Optimal Controller 93 x \ ( r ) : y{t) - y r tank level error x2(f): zero-meaned input disturbance = Fb{t) - Fm u ( t ) : Fu(f) - F m outgoing flow deviation Actual outlet flow F u(t) = u(i) + F m . Since x2{f) is the output of a low-pass filter with the cutoff frequency Ab, the system equation becomes, d_ dt *i(0 x2(t) = u(t) 0 Kp -Kp 0 -Ab 0 0 0 0 Xl(t) 0 0 x2(t) + 0 v(r) + h w(t), u(t) 1 "IT 0 (6.25) where as before K p is the process gain and Mt) is white noise. For broke tank level control, both the amplitude of u(t) and the rate of change v(t) are constrained. So the performance index J is J = Var[y(0] + P2 (Var[v(r)] + ^ Var[M(0]), (6.26) which is expressed with the two matrices Q and R, 1 0 0 0 0 0 0 0 up2 J = E{X(0' x(t) + v(t) v(t)}. (6.27) Q The optimal control v(r)* is a state feedback, v(f)* = - K x ( t ) = -kiXi(t) - k2x2(t) - hu(t) The algebraic Ricatti equation (ARE): P A + A T P -PBR-lBTP + Q = 0. Then, the state feedback gain K is, K = [h k2 k3] = BTP/p2 = [p l 3 p2i p33]/p2, where py is the ij-Hh element of P. Thus, K = [k\ k2 k3] is obtained as follows. (6.28) (6.29) (6.30) 1 Kplp{Ab + y/2Kp/p + Li) k\ = —, k2 = — P A2 + Abj2Kp/p + L> + Kp/p > (6.31) 6.4. Linear Optimal Controller 94 The equation v(r) = -Kx{t) is expressed in its Laplace transform, v(s) = su(s) = -k\X\(s) - k2x2(s) - k-$u(s). (6.32) Substitution of (3.24) into (6.32) yields, u(s) (j + k2 +k3) = xi (s) \[-ki - sk2Kp). (6.33) Let CL{s) denote the transfer function of the controller, u(s) = CL(s)(y{s)-yr) = CL(s)xi(s). (6.34) Then -k\ - sk2K-1 s + k\KJk2 s + b C L(s) = = -k 2K-; l^Li = K c , 6.35 s + k 2 + h s + k 2 + k 3 s + a where K c = -k 2 /K p , a = k 2 +k 3 , b = k\K p lk 2 . (6.36) The closed loop natural frequency OJc = -\jKplp and the sensitivity function S is S = (6.37) S2 + U>CS ^ 2 + fJ./0J2 + (O2 The closed loop damping factor 77 is n = 0.5 ^ 2 + H / O J 2 . (6.38) Thus, 77 increases from 0.5 VI towards 00 as more weight is given to Var[w(r)] (larger ji). The controller has two tuning parameters, p and /J. They are adjusted to obtain desired Var[y(0] and Var[w(r)]. An example applied for an industrial case is given in Section 7.3.2. 6.4.1 Jump System Theory In Khanbaghi et al. (2002), linear jump quadratic control theory was applied to the broke manage-ment problem. The state vector is x(t) := [xi(t) u(f)]T, and the system equation is, *i(0 0 -Kp X](0 0 v(0 + KpFb(t) = + «(/) 0 0 u{i) 1 0 (6.39) A B 6.5. Nonlinear Controllers 95 The performance index is (6.40) o The optimum control is a state feedback with bias, v(t)* = -K(x(t) + at) = -kxxi(t) - k2u(t) - Kae. (6.41) The feedback gain K is generally a function of £(0, K = K(£{f)). In this case, the system matrices A and B and the performance index matrices Q and R, which may be a function of £(?), are constant. So K does not depend on £(£), but the bias term a does change by £(/)• The two v(t)* of (6.28) and (6.41) are compared. The equation (6.28) has a term £2*2(0 that seems to be a variable. However x 2 (0 is constant for each state of £(t). The bias term ae in (6.41) serves the same purpose as xiif). Thus, the two equations are the same. If the same Q and R are used, there is no advantage in the jump linear system method. With the new method of evaluating the variances for each mode presented in Section 5.5, it is possible to adjust Q and R for each mode to obtain some desired results. Then, the feedback gain K is a function of £(r), and the controller becomes variable gain. This path is not pursued further because for broke tank level control, the nonlinear controllers of Section 6.5 would be more practical. 6.4.2 Implementation Issues Since the controller CL has no integrating action, the controller output u(t) must be base-loaded with the mean value of Fb(t), Fm as shown in Figure 6.3. The statistical nature of Ft,(t) may change with time due to grade changes or some other operational or equipment changes. Thus, occasional adjustments of Fm by some averaging or adaptive mechanism may be necessary. It is an important issue and could be a shortcoming in implementing CL, but is not pursued further in the thesis. 6.5 Nonlinear Controllers In this section, the input disturbance does not have to be zero-mean, and, u(t) = Fu(f). The broke tank level control has the following special characteristics: 6.5. Nonlinear Controllers 96 e(t) v s +B u(i) •c s + g Figure 6.3: Implementation of Q . The controller output w(f) is base-loaded with the mean incoming flow Fm. 1. Fb(f) during breaks F b 2 is much larger than the mean F m , and it is not practical to make the outlet flow u(t) equal to F b 2 to bring down the tank level. Usually, piping and pumps for the broke stock flow are not geared to handle flow rates as high as F b 2 . Also, from paper quality point of view, most operations will not allow such a high broke ratio. 2. The difference between Fm and Fb\ (Fb(t) during normal operation) is small, and to bring u(f) equal to F M can be done with small u(t). Thus, preventing the tank from becoming empty can be done easily with small flow roughness indicator. This means that the best the controller can do is to minimize the probability of tank overflows. In practice, for breaks of exceptionally long duration, normal operation of the paper machine is discontinued, and the broke tank level controller is put into manual mode. 6.5.1 Minimum Overflow Probability Controller One way of dealing with the flow smoothness issue, is to put hard limits on both u(t) and \u(t)\ = |v(r)|, and design control strategies to minimize the overflow probability under such constraints. From the binary nature of Fb(t), it is straightforward to design the minimum probability control strategy. Let u(t) and v(r) be constrained as follows: The key condition is that u m a x < Fb2 and u m i n < F b l , so during breaks the level increases (although we can control the rate of increase) and during normal operation, it is possible to prevent the level from going below the low limit by setting u(t) = F b \ . Usually it would be satisfactory to have symmetrical constraints on v(/) as vm i„ = - v m a x , but non-symmetrical constraints do not complicate Umin < U(f) < U, Vmin < V(t) < V, max- (6.42) the controller design. To minimize the overflow probability, 6.5. Nonlinear Controllers Ci ill). 97 during normal operation: the level should be brought down to the low l imits as fast as possible to maximize the storage capacity at the onsets of breaks, and during breaks: the level increase rate should be minimized. This means u(t) should be increased to its maximum value um a x as fast as possible. At the onset of a break, to make the level increase rate as low as possible, u(t) is increased by v m a x toward u m a x if u( t ) < u m a x . Thus, during normal operations, u{t) should be kept at u m a x to make the overflow probability low. If there is no limit on v(r), u(t) can be kept at umax until y(t) = yi, then u{t) is brought down suddenly to F M to keeper) at^. However, since v(r) is limited, u(t) must be decreased toward when gets close toyz,- The threshold value of y(t) to start decreasing u(f) is calculated below. Set the present time at 0, and denote the present level y(0) with yo and w(0) with «o- Let us assume that u(t) is decreased by the rate of v m i n < 0 to F M , The time to bring u(t) to F M , T is u ( t ) = u 0 + v m i n t _ F M - up Then the level at T, y{T) is, y(T)=yo + Kp f ( F M - u{i))dt = y0 + Kp f ( F M - u0 - vmint)dt KJ o J o (6.43) = + j-Z-iFbi ~ "o)2 To prevent low level limit violation, y(T) > y i must hold. Thus, a quadratic equation y ( t ) = y L - ^ - ( F b l - u ( t ) ) 2 (6.44) . Z V m / w determines if u(i) must be decreased to prevent low level violation. Equation (6.44) is a parabolic graph on the (u(t),y(t)) plane. If (u(t),y(t)) is above this graph, u(t) can be increased without violating the low level limit. When (u(t), y(t)) is below the graph, y(t) may become lower than yL with u(t) being decreased at the maximum rate of v m i n if normal operation lasts long enough. The logic of the minimum overflow probability controller (MOPC) is summarized in Figure 6.4. Figure 6.5 shows the phase plane where (u(t),y(t)) can move under the MOPC, whereyL = 0 and a is the parabola of decreasing u(t) corresponding to a in Figure 6.4. On a vertical line c, u(t) = umax, and 6.5. Nonlinear Controllers 98 normal: ifv(0<^-^-(FM-t/(0)2 decrease u{t) with the rate of v m i n <— parabola a else i f u { t ) < u m a x increase u(f) with the rate of v m a x . else keep u(t) at u m a x end end break: i f u ( t ) < u m a x increase u(t) with the rate of v m a x . else keep u(t) at u m a x . end Figure 6.4: The rninimum overflow probability strategy. b is the parabola of increasing u(i) for a break starting from (Fbi,>x). The point (w(0>.v(0) moves in region ^  and line c . A variation of the minimum overflow probability strategy is to keep u(t) as is when a break ends and the operation becomes normal. This way, u(t) stays constant most of the time. This strategy is called "quiet minimum overflow probability strategy" (QMOPC), and shown in Figure 6.6. 6.5. Nonlinear Controllers 99 break: if u(f) < umax increase u(f) with the rate of vmax. A else keep u(t) at,umax. end Figure 6.6: The logic of QMOPC 6.6. Optimal Change ofu(t) 100 6.6 Optimal Change of u(t) In the MOPC, u(t) is increased to um a x with the maximum rate allowed v m a x to minimize the level increase rate. So u(f) changes linearly with a slope of v m a x . The problem considered in this section is as follows. The initial value of u(f), UQ and the final value u\ are given. Also given is the time to move from UQ to u \ , T. What is the best way to change u(t) from UQ to u\ in a specified time T l Here, "best" means to minimize w(r)'s effects on the downstream process. A similar problem was treated in Dabros et al. (2003), where the performance index is the weighted sum of incremental changes in the pulp properties at the machine headbox. The performance index was calculated by simulations, and the best strategy was sought by adjusting parabolic curves. The downstream process is assumed to be a properly closed self-regulating loop as in Section 4.2. This is a rough model of real processes, but makes it possible to apply analytical tools (Max-imum Principle) and to evaluate various strategies in closed form. Let x\(t) be the downstream process variable. As in Section 4.2, the transfer function from u(s) to xi(s) is a high-pass filter with a cut-off frequency a . The process gain does not affect the optimum control, so for simplicity, the process gain at high frequencies is set to 1. ^4~r = —— => + axi(t) = ii(t) = v(t) u(s) s + a Without losing generality, the initial value is set to zero. X l(0) = 0; w(0) = w0 = 0 (6.45) It is required that at time T, u(T) must be equal to u \ , and stay at u \ after T. It is desired to change u(t) in a manner that minimizes its effects on x\(t). The performance index is J = £ ° x x { t ) 2 d t . (6.46) The control horizon is T, but / is integrated to oo. This is to account for the effects of u(t) after u(t) stops changing. If the downstream process is properly closed, x\(t) should return to 0 within a reasonable time. Thus, practically there is not much difference between a well chosen finite integral span longer than T and the infinite integral span. The infinite span is chosen because it 6.6. Optimal Change of u(t) 101 e-2an°° simplifies the mathematics. To keep u(t) constant after T, v(t) := u(t) = 0 for t e [T, oo]. Then J fT poo " pT poo xi(t)2dt+ xx(t)2dt = \ xx(t)2dt + [xx(T)e-a('-T)]2dt 0 JT Jo JT J f>T poo pT xx(t)2dt + xx{T)2 I [e-a(t)]2dt= f xx(t)2dt + xx(T) . 0 o Jo Jo Y~2a rT I = Xl(t)2dt + —xx(T)2 Jo 2or The problem is solved by Pontryagin's Maximum Principle. The solution involves singular arcs, which demonstrate the full power of the Maximum Principle. Mathematical details are deferred to the end of this chapter, and the main results are shown here. PROBLEM: Change the outlet flow u(t) from 0 to u i in a given time T to minimize J = xx (t)2dt, where X\(t) is the downstream process variable, subject to xi (0 = -axx (t) + v(t), xi (0) = 0 u(t) = v(t). The manipulated variable v(t) is bounded as -vm < v(t) < vm (vm > 0). Theorem 6.1 The optimum outlet flow u*(t) consists of three lines as shown in Figure 6.7. The slope of the first and the last line are the same and they are sgn(u\)vm, and the length r is the root of the following equation: T - e~aT(T - 2T) = — (6.47) Vm Quadratic programming would provide similar results after discretizing the problem, but cannot give the closed form solution given by the Maximum Principle. 6.6.1 Asymptotic Analysis Here the behaviour of r and J is studied when vm —> oo. Since r —» 0 as vTO —» oo, e~aT « 1 - ar. Without loosing generality, it is assumed that wi > 0, then, (6.47) becomes T - e~aT(T - 2T) * T - (1 - ar){T - 2T) = — (6.48) Vm So, T is a root of the following quadratic equation: 2orr2 - (aT + 2)T + ux fvm = 0. (6.49) 6.6. Optimal Change of u{t) 102 • i i i i 7/ *>> \\ \\ / f ^ ^ ^ ^ 10 12 14 16 18 20 time Figure 6.7: The graph of u*{t) at left and two responses of the downstream process variable x\(t); x*{(t) corresponds to u*(t), and the other corresponds to linearly changed u(f). Thus, T = aT + 2 ± yjiaT + 2)2-8auJv„ 4a (6.50) We pick a root that tends to zero as vm —> oo, so aT + 2- •s/(aT + 2)2-Saul/vm _aT + 2-(aT + 2) ^ 1 - g^ gfe 4a aT + 2 (. _ n _ 1 8aui/vm Finally, 4or As v„ 2 (orr + 2)2 Ml oo, r = U i 4or aT + 2 4au\/vm _ 4a ' {aT + 2) 2 " vJpTFTl) 0, and vTOr (6.51) Ml vOT(«r + 2) " """ aT + 2 ( 6 - 5 2 ) Thus, the optimum control v*(r) has a delta function of strength ux/(aT + 2) at t = 0 and r = 7\ The delta function causes the jump in x\(t),Ax\, Axi = lim Vm-»OO,T-»0 e a 'vmdf = vmT = aT + 2 With this v*(t), xi(0 has two jumps of u\/(aT + 2) at r = 0 and t = T. So, xi(0 = 0 "1 <*r+2 L ar+2 t = 0 o < r < r (6.53) Thus, , = , ' (6-54) + 2/ or(ai + 2) The above expression provides the infimum of J, against which all other control can be measured. 6.7. Summary 103 6.6.2 Linear u(t) In the minimum overflow probability strategy, u(f) is changed with a constant v(f) = v m a x . Let uc(t) denote the linearly changed u(f). "c(0 = Y1 And x l ( t ) = - ^ ( l - e - ° t ) (6.55) aT Let Jc denote J[°° xi (tfdt with uc(i), =(IT) 2 JX'0 - 2 e ~ a t + e ~ 2 a t ) d t + h { x ~ e'aT)2} ( 6 , 5 6 ) _ /_wi_\2 a T - 1 + e"a:r Vorr/ or Now compare Jmin and J c . / c _ / «i_\ 2 orr - 1 + e- a rd(o;7' + 2) _ ( a T + 2)(a7T + e ~ a T - 1) 5 /m ,„ var/ a- w2 a 2 r 2 When a T » 1, e - a r « 0, so J c/^,>, w 1 + ^ , and when orr « 0, e _ Q r : r * 1 - + ^(orT1)2, so dd Jmin ~ 1 + 5<z7\ Thus, at either extreme of orT" range, Jc/Jmin approaches 1 and it takes the maximum value 1.139 at a T = 2.7 as shown in Figure 6.8. In any case, the difference between the optimum u(f) and linear u(f) is small, and for practical purposes, the linear u(f) would be quite acceptable. 6.7 Summary The importance of broke storage tank level control for stabilization of a paper machine wet end was explained. The input disturbance Fb(t) can be modelled as a two-state continuous-time Markov process. The spectrum of Fb(t) was obtained as the output of a first-order low-pass filter driven by white noise. Many factors are involved for the flow smoothness requirements of u(t), and it was decided to use both Var[w(r)] and Var[v(r)] or max u(f) and max |v(x)| for the flow roughness indicator. 6.8. Maximum Principle 104 1.15 oT Figure 6.8: Jc/Jmin v s <xT. The optimal linear controller Q designed based on the input disturbance spectrum, was shown to be the same as the controller designed with jump linear system theory. From the nature of the input disturbance, it is not practical to expect a controller to keep the level below the high limit for exceptionally long breaks. Thus, a nonlinear controller that minimizes the overflow probability is proposed. The optimal change strategy of u(t), where u(t) is changed from its initial value wo to the final value «i in a given time T, is investigated. The performance index is the integrated square error of the downstream process variable. With a simple downstream process model, the problem is solved with the Maximum Principle's singular control. A closed form solution is obtained, which makes it possible to compare the linear change strategy to the theoretical optimum. 6.8 Maximum Principle The Maximum Principle (Pontryagin et a l , 1962) can be stated in various forms depending on the terminal conditions and fixed or variable control horizon. A simplified fixed control horizon version, where the end time T is prescribed, is used here following Loewen (2003). The control u(t) must be in a control set U(t), which can be a complicated set, but here U(f) is a simple closed 6.8. Maximum Principle 105 interval, i.e. «i(r) < u(t) < uH(f). The following three equations are given. system equation: x(t) = f(x(i),u(f),t), x(0) = xo (6.58) performance index: J = £o(x(T)) + I Lo(t, x(t), u(t))dt (6.59) J o terminal condition: Cj{x{T)) = 0, j = 1,... ,M (6.60) PROBLEM: find a control u(f) e U(t), t € [0, T] that minimizes J, such that x(t) moves from x0 to x (T) and satisfies the terminal condition (6.60). The adjoint vector p(t) and the Lagrange multipliers Aj are introduced. adjoint vector: p(t) same size as x(t) scalar Lagrange multipliers: Aj,j= 1,... ,M Hamiltonian: H(t,x,p,u) := p{t)Tf(t,x,u) -L0(t,x,u) (6.61) (6.62) (6.63) MAXIMUM PRINCIPLE: Suppose a control-state pair (u*(t),x*(t)) gives the minimum in the problem; assume that u*(t) is piecewise continuous. Then there exist AjS and a piecewise smooth vector function p(t) such that H(t, x*(t),p(t), u*(t)) is piecewise smooth, and one has (a) Adjoint equation: (b) Maximum condition: (c) Transversality condition: -*t)=lTx u(tY e arg max H(t,x* (f),p(t),v) veU(t) xr dUx*(T)) , v , d£j(x*(T)) -P(T)T dx +5>- dx 7=1 In our problem, the state vector x{t) = [xi(t) u(t)Y and the control is v(r). System equation: d X\(t) -axi(t) + v(t) dt u(t) v(0 Performance index: Terminal condition: So J = ± X l ( T ) 2 + 2a u(T) - «! = 0 Jo xx(t)2dt k = ^ - x x ( T f , L0 = X i ( t ) 2 , h = u { T ) - U i ( M = 1 ) 2ar 6.8. Maximum Principle 106 H = px(-ax\ + v) + p2v - x\ = v(pi + P2) - ap\X\ - x\ From (a), -Pi(0 = ~api(t) ~ 2xi(t) -HO = 0 From the transversality condition (c) -p(T) = [- X l (T) A\] => p\(T) = --xi(T) , p2(T) = a , a (6.64) A s p 2 ( t ) = 0 , p 2 ( t ) = - A l . The maximum condition (b) is v*(0 € a r g m a x i v [ p i ( t ) + p2(t)] - apx(i)x{(i) - xx{t): veU 1 With/?2(0 = -Ai, v*(0 G argmax {v|>i(0 - Ai]} veU As shown above, H is linear in v(f), so v(t) must be bounded to be able to construct the finite optimum control. We choose KOI < v m (6.65) Then v m if p i ( t ) - A { >0 v*(0 = - any valued [ - v m , v m ] if p i ( t ) -A i=0 (6.66) - v m if p \ ( f ) - A \ < 0 The state trajectory where p\(t) - A\ = 0 is called the singular arc. To investigate conditions for singular arcs, equations for p\(i) and x,(0 are studied. On a singular arc pxit) = a/?1(0 + 2xi(0 xi(0 = -axi(0 + v(0 pi(f) = A\ ~ constant (6.67) (6.68) (6.69) 6.8. Maximum Principle 107 F o r p \ ( t ) to be constant, x\(t) = -aA\/2. For x\(t) to be constant at -aA\/2, v(f) = -a2A\. There is a freedom to make v(t) to vm or - v m fromthe singular arc because both vm and -v O T are optimum as p\{t) - A\ follows the sign of v(r). v(0 -> vm => xi(t) T=> T=> MO - > 0 v(0 -> - v m => xi(t) |=> pi(0 |=> /?i(0 - Ax < 0 However, once and /?i leave the singular arc, they cannot come back to it. Thus, there is only one singular arc in v*(f). So the optimum control v*(f) consists of three parts: 1. forte [0,/i],v*(0 = vj = ±vm 2. for t e [tut2], v*(t) = v2 e [-vm,vm] 3. forte [t2,T],v*(t) = v3 = ±vm Parameters t\, t2, and v 2 are determined to satisfy the following two conditions: f v(f)dt = vih + v2(t2 + v3(T - t2) = ui Jo pi(T) = —xi(T) ,b Trivial solutions such as v*(t) = ±v m can happen when T and u\ are set up so that they are only feasible solutions. It is assumed that the problem is set up so that there is freedom for v*(r) to assume more than one value. On the singular arc, x\(i) is constant. Let x i c be that constant. As *i(0) = 0, xic := x,(f,) = v, P e a s ds = - ( 1 - e~ a t l ) = x,(/2) Jo a To keep *i(0 constant on the singular arc, from (6.68), v 2 = v([tut2]) = axXc = v,(l - e^' 1), and to keep p \ ( i ) constant, from (6.67), Pi([tut2]) = -2xlc/a = - e"a") Since x\(t2) = x\c, for t e [72, T] X l (0 = e" a ( ' - ' 2 ) x l c + v3 f * . . e ^ = e~ a { t- t 2\x l c - v3 /or) + v 3 /a (6.70) Jo ' So for t e [fe, 7/], p x(t) = apxit) + 2 [e- a ( t~ h\x l c - v3/or) + v3/or] 6.8. Maximum Principle 108 (6.71) Since p\(t2) = -2x\c/a, p i ( T ) = - 2 e a ( T - ' 2 ) x l c / a + 2 f e a ( T - s ) [ e - a { s - t 2 \ x l c - v 3 / a ) + v 3 / a ] d s ( 1 _ e-2a(T-t2) i _ e-a(T-t2) = -2e a { T - t 2 ) x X c la + 2e a ( - T - , 2 ) \ (x l c - v 3 /a) + v3/a {2a a The transversality condition (6.64) is p i ( T ) = - - X i ( T ) = > a P l ( T ) + x , ( T ) = 0 a Let T = T - h, then from (6.71), ; a p i ( T ) = -2e a T x l c + (eaT - e~ a T )(x X c - v3/a) + 2(e a r - l)v3/or And from (6.70), xi(T) = e~aT{xXc - Vila) + v3/or So ap x(T) + xi(T) = -e a Tx l c + v 3/a(e a T - 1) = 0 Substitution of xic = (1 - e~°"l)vi/a yields, - e a T { \ - e~ a t i)vi/a + (e a T - l)v 3/a = 0 We have two cases, v3 = vi and v3 = - v i . F o r v 3 = + v i - ( e a ( T ~ h ) - \ ) = 0 = > T = h a F o r v3 = -v, - ( - 2 e a T + e 0 ^ + I ) = 0 = > r = - - log(2 - e ~ a " ) a a Since 2 - e~a h > 1, log(2 - e ~ a ' ' ) > 0. This makes r non-positive. Therefore, for r to be positive, t\ = r, and vi = v3 (6.72) From the terminal condition u(T) = u\, V\t\ + V 2(T -ti-T)+ V 3T = U\ Substitution of v 2 = v t ( l - e~ a h) and r = t\ and v3 = vi yields, V l [ h + ( 1 - e a t x ) ( T - h - T ) + T ] = v x [ T - e a t \ T - 2 * 0 ] = ux. 6.8. Maximum Principle 109 So t \ = T is obtained as a root of the following nonlinear equation. T - e - m \ T - 2 t x ) = — v i (6.73) Since 0 < tx < 0.5T, T - e~atx(T - 2t{) > 0, thus wi/vi > 0. This means, Vi = v 3 = -v„ U\ > 0 Mi < 0 (6.74) Now with (6.72), (6.73) and (6.74), Theorem 6.1 is proved. Chapter 7 Simulation With Industrial Data 7.1 Introduction In this chapter, two broke tank level controllers, the quiet minimum overflow probability controller (QMOPC) and the optimal linear controller Q which are introduced in the previous chapter, are tested against an industrial data set from a paper machine via simulation. First, the data is analyzed for break and normal durations. Then, the two controllers are designed based on the break flow model estimated from the industrial data, and tested by simulating with the same data. 7.2 Data Most modern paper mills have multiple paper machines sharing more than one broke storage tank, some with complicated sharing arrangements. For this simulation study, a paper machine with a dedicated broke storage tank was selected to make it straightforward to estimate the incoming flow from the history of the storage tank level y(f) and the outlet flow u(i). 7.2.1 Process A simplified diagram of the broke management system of the paper machine where the data was collected is shown in Figure 7.1. Broke is diluted with water at the repulper, and fed to the broke storage tank. The consistency of the broke flow Fb(t) tends to be lower than the target consistency 110 7.2. Date 111 to the blend chest (around 3%). The consistency is increased by the thickener, and the thick broke is stored at the broke chest. From the broke chest, the consistency is controlled (not shown) and the controlled flow, u(t), is fed to the blend chest. Another stream from the thickener is water that is rich with fines and short fibres, called "cloudy water". The cloudy water is fed to the saveall, where fines and fibres are recovered and fed to the blend chest. The outlet flow u(t), which is controlled by the valve V I , directly influences the broke chest level. The broke chest level is controlled by the recirculation valve V2, through which u(t) in-fluences the broke storage tank level y ( t ) indirectly. There are some dynamics involved between u(f) and y ( t ) due to the broke chest level controller, but for a sampling rate of 5 minutes, it can be deemed that u(f) directly influences y ( i ) as, The paper machine produces various grades of light papers, and the broke contains substantial amounts of short fibres, fines and additives, which are recovered by the saveall. Thus, Ff,(t) in (7.1) is not a real incoming flow, but a hypothetical incoming flow seen from u(t) disregarding the material flow through the saveall and water taken out at the thickener. Figure 7.1: Simplified diagram of the broke management system. Broke is diluted with water at the re-pulper, and pumped to the storage tank. Water is removed from low consistency broke stock at the thickener. The high consistency broke stock is fed to the broke chest, which level is controlled by a recirculation valve V2. (7.1) Fbif) water brok< 7.2. Data 112 7.2.2 Data Set Three process variables, the broke storage tank level y(t), the broke stock flow u(f), and the break detector signal, were collected for a 6 month period starting in February 2002 with a sampling period of 5 minutes. The pattern of breaks changes a great deal from day to day as shown in Figure 7.2. In both graphs, it may be seen that the state of the break detector is not a good indicator for estimating the incoming flow to the broke storage tank, which is not measured. During the 6 month period, there were a number of shutdowns, scheduled and non-scheduled. The operation log was not available. Therefore, shutdown-like periods were detected when u(f) stays close to 0 for a prolonged period of time, and were subsequently taken out of the data set prior to the following analysis and simulations. Figure 7.2: Plots of tank level y{i) in % (bottom), outlet flow u{i) in GPM (top) and break detector (middle) for two 10 day periods. There were no breaks registered between day 50 and 53. 7.2. Data 113 7.2.3 Incoming Flow Estimate The incoming flow Fb(t) to the broke storage tank is estimated from measurements of the level y(t) and outlet flow u(f). From a theoretical material balance point of view, Fb(t) should be binary (one value for breaks and the other for normal operation). All the published studies assumed a binary broke flow, and the same assumption is used in Chapter 6. However, a cursory study of the data reveals that Fbif) is nowhere close to a binary flow. The binary flow assumption is discarded, but it is reasoned that Fbif) has distinct phases depending on the state of the paper machine, normal operation or break. Therefore, it is assumed that Fbif) is piecewise-constant, but not necessarily binary. The piecewise-constant Fb(t) assumption makes it sensible to obtain break duration statistics from Fbif). The sampled data value of Fb(t), y(t) and u(t) are expressed with time argument / or similar integer expressions, such as y(i) when there is no confusions. The piecewise-constant Fbif) is re-constructed from_y(0 and u(i) as follows. Suppose Fb(t) consists of & piecewise-constant segments. Let Fkb be the constant value for the A:-th segment: F b { f ) = F k b , for / = i k , ik+l, . . . , i k + n k - \ , where ik is the beginning of the k-th segment and nk the length. The continuous time model (7.1) is discretized with a sampling period At for each segment where Fb(i) is constant. ym(h +j)=/o + KpAtY/m=0[Fkb - u{ik + m)], j = 0 , . . . , nk-1, where ym{i) is the model level and y^ is an initial value of the level for the Ar-th interval. The two parameters y^ and Fkb and the length of the intervals nk are determined by the least square method. Let e,- be the model error for given F b and y^: e-x := ym(ik + i) - y(ik + i)-Then, yl + KpAtF* = y(ik) + KpAtu(ik) + e0 y* + 2KpAtFk b = y(ik+\) + KpAt}Zm =ou(ik + m) + e l y \ + n k K p A t F k = y ( i k + n k - \ ) + Zmli^ik + m ) + e „ k - i 7.2. Data 114 The two parameters y \ and F \ are determined to minimize the squared error £ e2. Thus, they are the solution of the following normal equation. 1 KpAt 1 2KpAt A' y(ik) + KpAtu(ik) y{ik+1) + KpAt £o w0'*+™) 1 nkKpAt . 1 . y(ik+nk-1) + £pAJ £0 "(**+»») When the least square solution of y \ and F b are obtained, the level errors e, can be calculated. et := y(ik + i) - KpAtYlm=Q[Fkb - u(ik + m)] ' i The length of an interval is increased until the maximum absolute level error exceeds the limit eum. The procedure is as follows: 1. k = 1 and i\ = 1 and «i = 2. 2. Obtain the least square solution ofy^ and Fk b and calculate the maximum absolute level error max,-3. If emax < eum, set n k = n k + \ and go back to step 2. 4 . If emax > eum, roll back nk asn k = n k -\ and Fb is the prior solution. 5. k = k+l,i k = 4_i + nk = 2, and goto step 2. The tank level is measured in feet. The tank height is 60 feet and inside diameter is 44 feet. The level ranged between 0 feet and 52 feet during the data collection period. Different e«w were tried. The value of 0.5 feet appeared to provide the best compromise. Figure 7.3 shows Fb(f), u(t) and y{i) for the same time period of Figure 7.2. There are clear abrupt changes in Fb(t), which indicate breaks. Figure 7.4 compares y(t) and y m(t) for a one-day period. It shows the piecewise-constant Fb(f) can explain y(f) quite well. 7.2. Data 115 10002 Figure 7.3: Graphs of estimated incoming flow Fb(f) (top) and level (bottom) for the same period as Figure 7.2. 7.2.4 Statistical Analysis The disturbance to the level controller is not breaks detected by the break detector, but changes in incoming flow. Thus, the break durations and intervals that are obtained from the estimated incoming flow Fb(t) are more meaningful. Figure 7.5 shows the histogram of Fb(t). It has a peak around 300 GPM, which is the centre value of normal operation. Any Fb(i) above 500 GPM or so would likely be for breaks. However, Fbijt) is widely distributed for breaks. A number of possible 7.2. Data 116 Figure 7.4: Comparison of the level measurement^/) and reconstructed level ym(i) from Fb(i) and u(i) for a one-day period. The y-axis is level in feet. explanations are: • Off-grade rolls are put into a repulper from time to time. The increase in Fb{i) due to off-grade rolls can vary depending on the quantity and the way they are repulped. • The paper machine produces different grammage (grades) of paper. • The operation of the thickener may change, which affects the amount of fibres going through the saveall. 10000 6000 4000 2000 Figure 7.5: Histogram of the estimated incoming flow Fbif). The n-th bin covers 100(w -1) < Fb{i) < 100« inGPM. 7.2. Data 117 The break durations and intervals (normal durations) are analyzed with the following rule: Fb(i) > B : break (7.2) Fb(f) < B : normal In the following analysis, B is set to 600 GPM. For break durations, equivalent durations were calculated assuming that Fb(f) is at some constant value Bf during breaks. The equivalent break duration r b is obtained as follows. Suppose that Fb(t) > B for t = [ t a , t b ] , then 1 r'b Tb : = — F b ( t )d t . That is, the equivalent break duration is a length of a hypothetical break which produces the same volume of incoming flow if the flow rate is constant at Bf. In this study, Bf is set to 3000 GPM. As is shown in Figure 7.5, 3000 GPM is the high end of Fb(f) for breaks. So, for most of the breaks, r b becomes shorter than actual ones. This makes the equivalent break characteristics more difficult to control than real incoming flow. Thus, controllers designed for the equivalent breaks are conservative, and less likely to cause overflows. 7.2.5 Probability Plots Probability plots are graphical tool to assess the fit of data to a theoretical distribution (Rice, 1988). Let x be a random variable with a strictly increasing probability distribution function FX(X) := P[x < X]. Let X,- be the i-th element of order statistics of N samples of x . That is, the N samples are arranged in ascending order, so that x \ < x 2 < • • • < x ^ - \ < x ^ . A random variable .y = Fx(x) distributes uniformly on [0,1]. Let v,- = FX(XJ). It is known that E[y,] = This suggests that a graph of Fx(xi) versus i/(N +1) should look linear if indeed x has the distribution function FX(X). In our case, FX{X) = 1 - e'™. Then 7, = 1 - e ' ^ ' . So, 1 - e~Ax' x —-— => logfl —j « -Axi . N + \ *\ N + l J Thus, the plot of x t versus log(l - i/(N + 1)) should look linear if x, distributes exponentially. The slope of the plots gives an approximation of A. Figure 7.6 shows graphs plotted with break or normal durations on the x-axis and log(l - i/(N + 1)) on the v-axis. Both graphs are close to a straight line. The least squares is used to obtain the slopes and offsets of the lines. The intensity 7.2. Data 118 minutes minutes Figure 7.6: The probability plots for break durations at left and for normal durations at right. Straight lines are obtained by least squares. parameter A's are obtained from the slope of the lines. These graphs support the exponential distribution hypothesis. 7.2.6 Chi-Square Test The chi-square test is a numerical method to accept or reject the hypothesis x ~ F x based on the difference between the sample distribution and its expectations. The range of the random variable x is divided into n cells C„ i = \ , . . . ,n, where C, := [c„ c M ) . The number of samples that lie in the j'-th cell, Oj is counted. Let be the expectation of 0,. If x ~ F x , E( = N \ Fx(x)dx. How far the sample distribution deviates from the expectation is measured by the following quan-tity, 2 _ y (Pr-Ejf If the hypothesis is correct, the sampling distribution o f ^ 2 is approximately the chi-square distri-bution with degrees of freedom of n - 2. Then, the expectation of x2 is approximately n - 2 . If X2 is much larger than n - 2 , the hypothesis is discarded. Figure 7.7 shows the 0,s in bar graph and plots of the corresponding model exponential distributions. The chi-square analysis is shown in Table 7.1. The break durations pass the chi-square test with P[y 2 > 12.15] = 0.5151, but the normal durations do not with P|jf2 > 37.98] = 0.0039. Although both break and normal data seem good in the probability plots. 7.2. Data 119 100 150 duration (minutes) 1000 1500 duration (minutes) Figure 7.7: The sample distributions (bar graphs) and corresponding model exponential distributions. At left is for break and at right for normal. Table 7.1: Chi-square test of break and normal durations break durations normal durations cell expected observed x1 cell expected observed x2 0 3 238.08 241 0.04 0 100 134.32 168 8.45 3 6 130.73 146 1.78 100 200 100.15 109 0.78 6 9 71.78 60 1.93 200 300 74.67 62 2.15 9 12 39.41 35 0.49 300 400 55.68 33 9.24 12 15 21.64 19 0.32 400 500 41.51 33 1.75 15 18 11.88 10 0.30 500 600 30.95 24 1.56 18 21 6.53 8 0.33 600 700 23.08 17 1.60 21 24 3.58 3 0.09 700 800 17.21 19 0.19 24 27 1.97 1 0.48 800 •900 12.83 18 2.08 27 30 1.08 2 0.78 900 1000 9.57 12 0.62 30 33 0.59 1 0.28 1000 1100 7.13 8 0.11 33 36 0.33 1 1.40 1100 1200 5.32 3 1.01 36 39 0.18 1 3.77 1200 1300 3.97 4 0.00 39 42 0.10 0 0.10 1300 1400 2.96 3 0.00 42 45 0.05 0 0.05 1400 1500 2.20 2 0.02 1500 1600 1.64 5 6.85 1600 1700 1.23 1 0.04 1700 1800 0.91 0 0.91 1800 1900 0.68 1 0.15 1900 2000 0.51 1 0.48 YX2 = 12.1538P = 0.5151 Yjc2 = 37.9827 P = 0.0039 7.2.7 Summary From the strength of the probability plots, the following model for the broke flow will be used for controller design. 7.3. Controller Design 120 i • t Ai 1/(398/60)= 1/6.633 (mean normal duration is 6.633 hours) A2 1/(26.2/60)= 1/0.437 (mean break duration is 26.2 minutes) FbX 311 G P M (normal flow rate) Fb2 3000 G P M (break flow rate) 7.3 Controller Design Although the estimated incoming flow Fb(t) is not quite binary, for controller design purposes, it is assumed that Fb(t) takes two distinct values, one for normal operation and the other for breaks. The level range of 0 to 50 feet is mapped to 0 to 100 % for the process gain Kp calculation. 7.3.1 QMOPC The MOPC and QMOPC of Section 6.5.1 have two design parameters, the maximum outlet flow Umax and the maximum rate of change of u(t), v m a x . They are specified with the mean incoming flow fm as a unit. The mean incoming flow fm is calculated as fm=Fl A2 M' + F> bl- = 477GPM. A\ + A2 A\ + A2 The overflow probabilities for both the MOPC and QMOPC are calculated from the state probabil-ity distributions obtained by the integral equation method of Section 5.4 with various combinations of the maximum outlet flow um a x and the maximum rate of change v m a x , and is shown in Table 7.2. The conditional probability densities of the level y(t), 0i(y) for normal operation and (f>2(y) for Table 7.2: Overflow probabilities for QMOPC and MOPC P[y(r) > 100] in % during breaks Umax Vmax QMOPC MOPC 2fm 2fm 0.0319 0.0297 2fm fm 0-0421 0.0358 2fm 0.5fm 0.0801 0.0500 1.5/; 2fm 0.3124 0.2887 1.5/* fm 0.3531 0.3065 1.5/M 0.5fm 0.4491 0.3421 7.3. Controller Design 121 breaks are numerically obtained. Then, overflow probability is calculated as P[y(r)>100]= f°° <f>2(y)dy. Jioo During a break, u(t) < Fb2, thus, y{t) is increasing. Therefore, y(t) has local maximums at the end of breaks. When v(r) at the end of a break is below 100%, y(t) stays below 100% during that break. By Theorem 5.1, the conditional probability.distribution of.y(r) for breaks, is the same as that of y(t) at the end of breaks. Thus, for instance^ if P[y(f) > 100] = 1%, the level y(t) stays below 100 % on average 99 times out of 100 breaks. The overflow probabilities tabulated are the conditional probabilities on breaks. The average number of breaks per day is, 24/(6.633+0.437)=3.39. For instance, if P[y(0 > 100] = 0.8%, there will be, on average, 3.39 x 365 x 0.008 = 9.9 overflows per year. There are not much differences between the MOPC and QMOPC, and it was decided to use the QMOPC, which keeps u(t) constant most of the time during normal operation, and is superior to the MOPC in stabilizing the paper qualities. More combinations of um a x and v m a x are tried for the QMOPC as shown in Table 7.3. For simulation studies, a (umax, vmax) combination with the Table 7.3: Overflow probabilities for QMOPC P[y(0 > 100] in % during breaks 0.611 0.659 , 0.801 . ~J , 1.237 ;':\">.1 2.920 1.377 3.532 overflow probability near 1 %, that is (\Afm,0.5fm) was selected. Suppose u(t) is changed from Fbi to umax with the rate of Q.5fm, then it takes 7= 1.496 hours. Assume that the consistency loop downstream is tuned to a time-constant of 3 minutes, so it has a cut-off frequency a = 20 hour -1. Substitution of these values to (6.57) yields A = ( < * T + 2 ) ( a T + e - ° r - l ) = { m $ Jmin O2 T 2 So, changing u{t) with a constant rate increases the consistency squared error 5 % over the optimal case. The probability distribution of the tank level is shown in Figure 7.8. The low level limit yL is set to zero. During normal operation, there is a probability mass at y = yL (y(t) stays at yL for ttmax Vmax lAfm 2fm lAfm fm 14fm 0.5fm 1.4/„ 0.25/m 1.4/M o.i/. 1.3/. 2fm 1.2/m 2fm 7.3. Controller Design 122 Figure 7.8: Probability distribution of XO under the QMOPC. At left, a plot marked B is the conditional dis-tribution for break, and N is for normal without a probability mass at 0. At right, the probability distribution for normal with the probability mass at 0. durations longer than zero). Thus, the graphs are not for the probability density functions, but the probability that y(t) e [iAy,(i + I)Ay] (Ay<= 1). 7.3.2 Linear Controller The optimal linear controller of Section 6.4 is designed to compare to the QMOPC. The perfor-mance index J = Var[y(0] + p 2 (Var[w(0] + pVax[u(t)]), where p and p are two design parameters, to constrain both u(t) and u(t). The closed loop natural frequency o)c is determined by p and K p as o»c = yjKpIp and p determines the closed loop damping factor r\-\ yjl + p/aj 2.. W h e n p = 0, n = 0.5 A/2, and the controller minimizes Var[«(0] with the largest Var[w(0]. Increasing n with higher p decreases Var[«(0] but increases Var[w(0]. Figure 7.9 shows changes in the standard deviation of u(t), cr[u] and that of u(t), cr[u] for r\ between 0.5 V2 and 5. They are shown relative to the values for p = 0, keeping cr[y] constant. Let o-0[u] and <TQ[U] be the standard deviations of M(0 and w(0 for p = 0. The plots are cr[u]/cr0[u] (lefty-axis) and o-[u]/o-0[u] (rights-axis). It is observed that cr[u] increases linearly within but o~[u] reduction tapers off quickly. As the best compromise, 77 = 2 was chosen since it yields most of the cr[u] reduction with a modest increase in cr[u]. The graphs were drawn for cr\y] = 20%, but they are almost identical for other values of cr[y]. It was found from simulations with the industrial data, cry = 20% is necessary to obtain 7.3. Controller Design 123 damping factor Figure 7.9: cr[u] decrease and cr[u] increase versus the damping factor JJ for a constant cr[y] comparable performance with the QMOPC. The controller has the following transfer function, CL = KC . s + a Table 7.4 shows the controller parameters and o~[u] and o~[u]. Table 7.4: The controller parameter for cr[y] = 20 n Kc a b <r[u] cr[w] 0.5 V2 0.0807 0.1316 2.444 24.48 2.344 1 0.1121 0:2204 2.446 22.35 2.521 2 0.3232 0.7524 2.452 20.61 3.950 5 1.8145 4.4790 2.471 20.10 9.101 The probability density of y{t) and u(f) calculated by the method of Section 5.4 are shown in Figure 7.10. 7.4. Simulation 124 Figure 7.10: Probability density functions of y(t) at left and u(t) at right under Q, 7.4 Simulation The QMOPC and the linear controller CL specified in previous sections are tried with the estimated incoming flow Fbif) obtained from the mill data. The controllers are designed assuming the binary incoming disturbance, but subjected to the estimated Fbif), which is not binary. The original data set is 180 day long, but after taking out shutdown days, the total length became 155 days. For both the QMOPC and CL, yif) and u{t) are reset to the original manual control value after each shutdown period. Then, u(t) are calculated according to each controller, and y(t) are updated based on the process model, which is a discrete time version of (7.1). Table 7.5 lists statistics of the three cases, manual, QMOPC and CL. The linear controller CL has a very small Var[v(x)] compared to the other two. Figure 7.11 plots the outlet flow u(t) of the three cases for a short period (two days) to clearly show the flow smoothness of each case. Table 7.5: Statistics of simulations y(t) (%) i/(0 (GPM) <r[v(0] control specification „ max min max min GPM/H manual 104.4 0.06 1004 1.8 334.8 QMOPC Umax ~ Idfmy ^max = ^*^fm 102.5 1.92 668 0 88.91 CL n = 2, o-[y(0] = 20 105.1 -23.7 875 154 13.7 7.4. Simulation 125 1000 Figure 7.11: Comparison of the QMOPC, CL, and manual control 7.4.1 Manual Control The histograms of XO and u(t) under manual control are shown in Figure 7.12. The histogram of XO has a clear peak at a 35 % to 40 % range. Probably this range is preferred because there is no imminent danger of either overflow or emptying the storage tank. Thus, the operator does not have to pay close attention to the broke storage tank level. The histogram of u{i) has a broad peak ranging between 350 GPM and 600 GPM. However, u{t) spends substantial time at both high and low ends reflecting the operators' reaction towards high and lowXO- The manual control yields the highest o"[v(0] as a result of the operators making step changes in u(t). 7.4.2 Linear Controller The histograms by the linear controller CL have a similar pattern to those by manual control. However, the percentage of time that u{t) spends in high and low ends are much smaller compared to manual control. This shows some non-linear characteristic of manual control. CL needs the level setpoint jv, which was set to 30 %. From Figure 7.10, ayr of 40 % seems a safe choice, but 30 % was selected to give CL an advantage in controlling high X0S- As a result, X0 spends some time below zero. In practice, it is possible to add low level constraining logic similar to one for the 7.4. Simulation 126 .On J J L 100 200 300 400 500 600 700 800 900 1000 u(t) (GPM) Figure 7.12: Histograms of y(i) (left) and u(t) (right) by manual control QMOPC. However, CL was run as a pure linear controller to see the limitations of a linear scheme. CL needs its calculated control output base-loaded with the mean incoming flow fm as mentioned in Section 6.4.2. It was found during the simulations, that the chance of low level limit violation is sensitive to the value offm. This could be a problem in practice, which must be solved by some adaptive control logic adjusting the value of the base load, or the low level constraining logic. CL has very small cr[v(r)] with more swing in ii(f) than manual control or QMOPC. As was shown in Figure 7.9, cr[u(f)] cannot be made much smaller with higher <x[v(r)]. This is a limitation of linear controllers. Although <x[v(0] is very small, u{t) is changing all the time as shown in Figure 7.15. The effects of this kind of broke stock flow changes on paper quality and machine runnability are not known. D a 8000 7000 6000 6000 4000 3000 2000 1000 D level (%) 100 200 300 400 500 600 700 u(l) (GPM) 900 1000 Figure 7.13: Histograms ofy(0 (left) and u(t) (right) by CL 7.5. Summary 127 7.4.3 QMOPC 100 200 300 400 500 600 u(t) (GPM) 700 800 900 1000 Figure 7.14: Histograms of y(t) (left) and u(t) (right) by QMOPC The QMOPC shows a very different pattern of behaviour from either manual control or C L . For y{i), it is concentrated at a low range to maximize the tank capacity at the onsets of breaks. This can be observed from time series plots of Figure 7.16. It is difficult to keep y{t) this low without emptying the tank by either linear or manual control. To make the QMOPC robust for low level limit violation, the low level limit yL was set to 5 %. So y(f) is unlikely to go below zero even for very small incoming flow Fb(t). For u(t), it is concentrated at the high limit umax, which is considerably lower than the maximum value of both manual control and Q . Compared to CL, the QMOPC has much higher cr[v(r)], however still much lower than that of manual control. If the rate of change in u{f) has some influence on break frequency, the QMOPC should yield less frequent breaks than manual control. However, only way to test this hypothesis is mill trials. For practical matters, the QMOPC does not require the level setpoint, and also is not sensitive to the value of fm for low level violation due to the built-in logic. These characteristics of the QMOPC make it easy to implement. 7.5 Summary The incoming flow to the broke storage tank Fb(f) was estimated from industrial data from a 6 month period. The break and normal durations were obtained from Fb(f) and they follow the exponential distribution well. The optimal linear controller CL and hard constrained nonlinear 7.5. Summary 128 days Figure 7.15: Comparison of u{i) by the QMOPC and manual (upper) and the QMOPC and CL (bottom). controller, QMOPC were designed for the statistical characteristics of Fb(i). Their overflow prob-abilities were calculated from the integral equation approach of Chapter 5. The optimal linear controller produces very smooth u(t), but the maximum value of u(t) is approximately 50 % higher than that of the QMOPC. Also, selection of the level setpoint and base load could be difficult problems in practice. The QMOPC can control the level well with a smaller maximum u(t), but the rate of change (Var[v(x)]) is higher than C L , although much smaller than manual control. It has none of the implementation difficulties of CL and is recommended for field • i : , i 1 • tests. 7.5. Summary 129 31 32 33 34 35 36 37 38 39 40 days 31 32 33 34 35 36 37 38 39 40 days Figure 7.16: Comparison of y{i) by the QMOPC and manual (upper) and the QMOPC and CL (bottom). Chapter 8 Conclusions 8.1 Summary Averaging Level Control A design methodology for stochastic averaging level control is proposed. The problem is formu-lated to mirrimize an outlet flow roughness indicator r ( u ( t ) ) subject to the target level variance Var[y(0]- When r ( u ( t ) ) is a linear combination of Var[w(Y)] and Var[w(0], standard linear quadratic optimal control theory is applicable, which produces a state feedback controller. With a noise-free observer, an easily implementable scalar transfer function for the controller is obtained. The importance of input disturbance modelling and proper selection of the roughness indicator is high-lighted. • The random-walk model for the input disturbance, which is popular for regulatory control design, has a few shortcomings when applied for the averaging level control problem. Since Var[w(0] may not exist, it cannot be used in r ( u ( t ) ) . When r ( u ( i ) ) = V a r [ u ( t ) ] for the random-walk input, the optimal controller is a proportional+integral action (PI) controller. • A first-order low-pass process, which is the output of a first-order low-pass filter driven by white noise, is proposed as a more realistic model for the input disturbance. When r ( u ( t ) ) is a weighted sum of Var[w(0] and Var[w(0], the optimal controller is a lead-lag network {C L ) . The performance improvements by C L over the PI controller could be substantial if the input disturbance cut-off frequency is high relative to the closed-loop level control loop band width. 130 8.1. Summary 131 • When r{u{ij) = Var[u(t)], that is, it is desired to minimize the swing of u(t), and the input dis-turbance is a first-order low-pass process, the optimal controller is a proportional+derivative type. However, compared to a proportional only controller (P controller), the performance improvement is not substantial, and for practical purposes, the P controller would be a good choice. Since the purpose of the flow roughness indicator r(u(tj) is to express adverse effects of u(t) on the downstream process variable yiif), the case of r(u(t)) = var{y2(t)] is investigated. The downstream process is modelled as a first-order high-pass filter (cut-off frequency = a>p) driven by u(f). The optimal controller for the first-order low-pass input disturbance, is a second-order lead-lag network with one of its zeros at -a)p. The performance improvements over CL can be as much as 40 % reduction in Var[y2(0]- When iop/cod » 1 (o»</ is the cut-off frequency of the input disturbance), the performance improvement is small. Markov Jump Systems A continuous-time Markov jump system is a hybrid system with a continuous state x(t) and a dis-crete state £(r), where £(r) is a continuous-time finite-state Markov chain. A new theorem, which shows that the conditional probability distribution of x(t) given £(f) = £ is the same as the proba-bility distribution of x(f) sampled just before £(t) leaves state £, is presented. From this theorem, two practically important applications are derived: calculations of state probability distributions of jump systems, and calculation of the state moments of linear jump systems. The probability den-sity function of a filtered binary random signal is obtained analytically. Unlike approaches based on simulations, the new method does not have time-discretization approximations nor randomness due to computer generated random numbers. This makes jump Markov system control theory more applicable to practical problems. Broke Storage Tank Level Control Sheet breaks cause the biggest upsets in the paper machine operation. The broke storage tank is the first process to receive the upsets, and to manage the broke tank level properly goes a long way towards stabilizing the paper machine wet end. The incoming flow to the broke storage tank Fbif) can be modelled as a two-state continuous-time Markov chain. The spectrum oiFb{t) is obtained, 8.2. Future Research 132 which is a first-order low-pass process. It is difficult to find a good flow roughness indicator r(u(t)) because there is no identifiable single downstream process. However, a weighted sum of Var[w(0] and Var[zi(f)] is adopted as r(u(t)). Then, the .optimal controller is a lead-lag network CL-It was realized that since Fb(i) during breaks Fb2 is so much larger than normal flow FbX, it is not practical to design a controller that can keep the level below the high limit with probability 1. A practical nonlinear controller with hard constraints on the outlet flow u(t) and its rate of change it(t) is presented. The controller, which is called MOPC, minimizes the overflow probability and guarantees that the level will stay above the low limityL- The new method for calculating the state distribution of Markov jump systems makes the design process straightforward without relying on simulations. The actual incoming flow to a broke storage tank Fb(i) was estimated from the tank level and outlet flow measurements of a paper machine. It was confirmed that the Markov hypothesis for Fb(t) does hold for the industrial data. The MOPC is compared to the optimum linear controller and manual control with the industrial data. The MOPC controls the tank level better than both the linear controller and manual with smaller maximum outlet flow. It is recommended to conduct a mill trial for the MOPC. The optimal change strategy of u(f), where u(f) is changed from its initial value UQ to the final value «i in a given time T, is investigated. The performance index is the integrated square error of the downstream process variable. With a simple downstream process model, the problem is solved with the Maximum Principle's singular control. The closed form solution is obtained, which makes it possible to compare linear change strategy to the theoretical optimum. For process parameters for usual broke storage tank systems, it is shown that linear change strategy is not much different from the theoretical optimum. 8.2 Future Research Constrained Stochastic Control In Chapter 3, the problem was formulatedtQi'minimize some outlet flow roughness indicator r(t), such as Var[zi(0], subject to a given Var[y(*)]- To reflect the true nature of the averaging level control problem, the problem should be formulated to minimize r(u(t)) subject to yL < y(t) < yH-8.2. Future Research 133 However, there are a few obstacles for this approach. • To guarantee yL < y(t) < yH with probability 1 with finite u(t), the input disturbance must be modelled as a bounded stochastic process. This should not be a problem because all the actual input disturbances are bounded. However, random walk or Gaussian random process, which are unbounded, are popular choice due to their mathematical tractability. • The problem itself is very difficult as existing optimization techniques, such as quadratic programming, are not directly applicable. The present constrained control assumes step-like deterministic disturbances. In some process control applications, this assumption is quite different from reality and a performance improvement would be possible if a controller could be designed for a more realistic disturbance model. The problem is tough, and probably needs a theoretical breakthrough or two before effective solu-tions become available. Markov Jump System • The input disturbance is modelled as a piecewise constant random process. Inclusion of con-tinuous noise (Wiener process) in the disturbance term will make the model more useful for some applications. More sophisticated mathematics is necessary, but the extension should be straightforward. • The integral equation for the state distribution is solved with a simple discretization method. A more elaborate method that makes impossible to obtain accurate solutions of the integral equation with fewer grid points is desirable. This would shorten calculation times and lessen memory requirements for higher dimensional systems. • It would be helpful i f the error bounds of numerical solutions of the integral equation were known. Then the grid size could be selected according to the required accuracy. Broke Storage Tank Level Control Although the MOPC has been tested against the industrial data, a number of improvements may be necessary to make the MOPC accepted in the field. 8.2. Future Research 134 • In the simulation, the incoming flow is estimated with the least square method. This approach is not applicable in real time control because the future measurements cannot be used. A reliable method to estimate Fb(f) from the measurements available up to the present time must be devised. More challenging'would be a case where multiple paper machines share one or more broke storage tanks. Not only the controller but also the incoming flow estimator must be designed for a MIMO system. • It is assumed that the statistical characteristics of the input disturbance are constant, and the industrial data set shows good stationarity. However, non-stationarity has been reported in plastic film processes by Taylor and Duncan (2002), and an "avalanche of breaks", which suggests that the machine operation has bad spells, as mentioned in Lindstrom et al. (1994). Therefore, the MOPC may have to be equipped with some adaptive mechanism to change the maximum outlet flow umax or the ramp rate based on the situation. The above two problems do not seem insurmountable with the present technology, and field tests would provide a definite answer about the practicality and effectiveness of the MOPC. Bibliography Allison, B. J. and M. Khanbaghi (2001). Model predictive optimal averaging level control with state estimation. In: Proceedings of the IASTED International Conference: Control and Ap-plications. Banff, Alberta, Canada, pp; 88-93. o Astrom, Karl (1970). Introduction to Stochastic Control Theory. Academic Press. New York, NY. Batina, Ivo, Anton Stoorvogel and Siep Weiland (2002). Optimal control of linear, stochastic sys-tems with state and input constraints. In: Proc. 41st IEEE CDC. Las Vegas, NV. Bonhivers, J-C, M. Perrier and J. Paris (1999a). Management of broke ratio in an integrated pa-per mill. Post-Graduate Research Laboratory Report PGRLR725. Pulp and Paper Research Institute of Canada. Bonhivers, J-C, M. Perrier and J. Paris (19996). Model predictive control for integrated manage-ment of Whitewater tanks. Post-Graduate Research Laboratory Report PGRLR724. Pulp and Paper Research Institute of Canada. Bussiere, S., A. Roche and J. Paris (1987). S,tudy of broke handling and white water management using a dynamic simulation. Pulp & Paper Canada 88(11), 74-77. Campo, P.J. and M. Morari (1989). Model predictive optimal averaging level control. AICheE Journal 35(4), 579-591. Cheung, T. and W.L. Luyben (1979a). A further study of the pi level controller. ISA Transactions 18(2), 73-77. 135 Bibliography 136 Cheung, T. and W.L. Luyben (19796). Liquid-level control in single tanks and cascades of tanks with proportional-only and proportional-integral feedback controllers. Ind. Eng. Chem. Fun-dam. 18(1), 15-21. p . . u Costa, O. L. V., E.O. Assumpcao Filho, E.K. Boukas and R. R Marques (1999). Constrained quadratic state feedback control of discrete-time markovian junmp linear systems. Automat-ica 35, 617-626. Crisafulli, S. and R.D. Peirce (1999). Surge tank control in a cane raw sugar factory. Journal of Process Control 9, 33-39. Croteau, A.P. and A. A. Roche (1987). Study of broke handling and white water management using a dynamic simulation. Pulp & Paper Canada 88(11), 74-77. Cutler, C.R. (1982). Dynamic matrix control of unbalanced systems. ISA Transactions 21(1), 1-6. Dabros, Michal, Michel Perrier, Fraser Forbes, Martin Fairbanks and Paul Stuart (2003). Improving the broke recirculation strategy in a newsprint mill. In: Proc. PAPTAC 89th Annual Meeting. Montreal, PQ. Doss, J. E., T. W. Doub, J. J. Downs and E. F. Vogel (1983). New directions for process control in the eighties. In: AIChE Meeting. Washington, DC. Feller, William (1968). An Introduction to Probability Theory and Its Applications I, Third Edition. John Wiley & Sons. New York, NY. Feller, William (1971). An Introduction to Probability Theory and Its Applications II, Second Edi-tion. John Wiley & Sons. New York, NY. Fleming, W. H. and R. W. Rishel (1975). Deterministic and Stochastic Optimal Control. Springer-Verlag. New York, NY. Foley, M. W., K. E. Kwok and B. R. Copeland (2000). Polynomial lq solurtion to the averaging level control problem. In: 16th International Conference on CAD/CAM, Robotics and Facto-ries of the Future. Port of Spain, Trinidad, West Indies. Bibliography 137 Gess, J.M. and R.A. Kanitz (1996). Monitoring the stability of a paper machine. Tqppi Journal 79(4), 119-126. Grimmett, G. R. andD. R. Stirzaker (1992). Probabiliy and Random Processes, 2nd edition. Oxford University Press. Oxford, UK. Harris, T. J. and J. F. MacGregor (1987). Design of multivariable linear-quadratic controllers using transfer functions. AIChE Journal 33(9), 1481-1495. Ji, Yuandong and Howard Jay Chizeck (1990). Controllability, stabilizability, and continuous-time markovian jump linear quadratic control. IEEE Trans. Automat. Contr. 35(7), 777-788. Kao, Edward P. C. (1997). An Introduction tOiStochastic Processes. Duxbury Press. Belmont, CA. Kelly, J. D. (1998). Tuning digital pi controllers for minimal variance in manipulated input moves applied to unbalanced systems with delay. The Canadian Journal of Chemical Engineering 76, 967-975. Khanbaghi, M., R. Malhame and Perrier (2002). Optimal white water and broke recirculation poli-cies in paper mills via jump linear quadratic control. IEEE Tans. Control Systems Technology 10(4), 578-588. Khanbaghi, M. , R. Malhame, M . Perrier and A. Roche (1997). A statistical model of paper breaks in an integrated tmp-newsprint mill. Journal of Pulp and Paper Science 23(6), 282-288. Khanbaghi, Maryam (1998). Un Formalisme de Systemes A Sauts Pour La Recirculation Optimale Des Casses Dans Une Machine A Papier. PhD thesis. Ecole Polytechnique De Montreal. Kropholler, H., Wiseman N. and Harris E. (1994). Theory and practice in modelling paper machine white water systems. Appita 47(4), 301-304. Lindstrom, R., W.H. Manfield, A.F. Tkacz and J. Mardon (1994). Coping with an avalanche of breaks. Appita 47(2), 163-172. Liu, Zikuan and Fengsheng Tu (1999). Single sample pathe-based sensitivity analysis of Markov processes using uniformization. IEEE Trans. Automat. Contr. 44(4), 872-875. Bibliography 138 Loewen, Philip D. (2003). M403 lecture notes. University of British Columbia Mathematical De-partment. Luyben, W.L. and P. S. Buckley (1977). A proportional-lag level controller. Instrumentation Tech-nology 18(1), 65-68. • Mariton, Michel (1990). Jump Linear Systems in Automatic Control. Marcel Dekker, Inc.. New York, NY. Marlin, Thomas E. (1995). PROCESS CONTROL. Chap. 18 Level and Inventory Control. McGraw-Hill, Inc.. New York, N.Y McDonald, K.A. and TJ. McAvoy (1986). Optimal averaging level control. AICheE Journal 32(1), 75-86. Morari, M. and E. Zafiriou (1989). Robust Process Control. Prentice-Hall, Inc. Englewood Cliffs, N.J. Newton, G. C , L. A. Gould and J. F. Kaiser (1957). Analytical Design of Linear Feedback Con-trols. John Wiley & Sons, Inc. New York, NY. Ogawa, Shiro (1995). PI controller tuning for robust performance. In: Proc. 4th IEEE Conference on Control Applications. Albaby, N.Y. Ogawa, Shiro, Bruce Allison, Guy Dumont and Michael Davies (2002). A new approach to aver-aging level control. In: Proc. Control Systems 2002. Stockholm, Sweden. Orccotoma, J.A., D. Stiee, J. Paris and M. Perrier (1997a). Dynamics of fines distribution in a white-water nework. Pulp e> Paper Canada 98(9), 77-80. Orccotoma, J.A., J. Paris, M. Perrier and A. A. Roche (19976). Dynamics of white water networks during web breaks. Tappi Journal 80(12), 101-110. Pontryagin, L. S., V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mishchenko (1962). The Mathe-matical Theory of Optimal Processes^. Interscience Publishers. New York, NY. Bibliography 139 Rice, John A. (1988). Mathematical Statistics and Data Analysis. Wadsworth & Brooks/Cole Ad-vanced Books and Software. Pacific Grove, California. Shinskey, F. G. (1967). Process Control Systems. McGraw-Hill, Inc. New York, NY. Shinskey, F. G. (198 8). Process Control Systems, Third Edition. McGraw-Hill, Inc. New York, NY. Shinskey, F. G. (1994). Feedback Controllers for the Process Industries. McGraw-Hill, Inc. New York, NY. Shunta, J. P. and W. Fehervari (1976). Nonlinear control of liquid level. Instrumentation Technol-ogy 23, 43^18. Sidhu, Manpreet S. (2003). Master thesis, multivariable averaging level control. University of British Columbia. Stewart, W. J. (1994). Introduction to the Numerical Solution of Markov Chains. Princeton Uni-verstiy Press. Princeton, New Jersey. Sworder, David D. (1969). Feedback control of a class of linear systems with jump parameters. IEEE Trans. Automat. Contr. 14(1), 9-14. Sworder, David D. (1976). Control of systems subject to sudden change in character. Proc. IEEE 64(8), 1219-1225. Sworder, David D. and Lin L. Choi (1976). Stationary cost densities for optimally controlled stochastic systems. IEEE Trans. Automat. Contr. 21(4), 492-499. Sworder, David D. and Ronald O. Rogers (19.83). An LQ-solution to a control problem associated with a solar thermal central receiver. IEEE Trans. Automat. Contr. 28(10), 971-978. Taylor, A. R. and S. R. Duncan (2002). A comparison of techniques for monitoring process faults. In: Proc. Control Systems 2002. Stockholm, Sweeden. Tijms, H. C. (1994). Stochastic Models An Algorithmic Approach. John Wiley & Sons. Chichester, West Sussex, England. Bibliography 140 Wonham, W. M. (1970). Probabilistic Methods in Applied Mathematics, A. T. Bharucha-Reid ed.. Chap. Random Differential Equations in Control Theory, pp. 191-199. Vol. 2. Academic Press. New York, NY. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items