Identification of nonlinear control models using volterra-laguerre series

Linear model predictive control has been widely accepted in industry as an important tool for the operation of difficult interacting processes. Linear identification and control techniques are well developed and well understood. In the industry, it is rare to find a system that is truly linear. While for many systems linear modeling and control can approximate their performance in certain regions, there exist some systems whose nonlinearity is great enough that an approximate linear model and control scheme cannot yield the desired accuracy. In order to control these more complex nonlinear systems, significant research has been dedicated to extending model predictive control to nonlinear systems.

IDENTIFICATION OF NONLINEAR CONTROL MODELS USING VOLTERRA-LAGUERRE SERIES by Dale A. Smith A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah August 2010 Copyright © Dale A. Smith 2010 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of Dale A. Smith has been approved by the following supervisory committee members: Philip J. Smith , Chair 17 May 2010 Date Approved J.D. Seader , Member 17 May 2010 Date Approved Terry Ring , Member 17 May 2010 Date Approved Geoffrey D. Silcox , Member 17 May 2010 Date Approved Mark Minor , Member 17 May 2010 Date Approved and by JoAnn Lighty , Chair of the Department of Chemical Engineering and by Charles A. Wight, Dean of The Graduate School. ABSTRACT Linear model predictive control has been widely accepted in industry as an important tool for the operation of difficult interacting processes. Linear identification and control techniques are well developed and well understood. In the industry, it is rare to find a system that is truly linear. While for many systems linear modeling and control can approximate their performance in certain regions, there exist some systems whose nonlinearity is great enough that an approximate linear model and control scheme cannot yield the desired accuracy. In order to control these more complex nonlinear systems, significant research has been dedicated to extending model predictive control to nonlinear systems. The problem of implementing nonlinear model predictive control can be split into two main tasks: making the nonlinear model and calculating control inputs. The significant contributions of this dissertation are in the area of identification of nonlinear Volterra models from input-output data. Historically, the identification of Volterra models has been limited to lower order models because of the large amount of model parameters that need to be identified. By using the Laguerre polynomials, the number of model parameters can be greatly reduced, which limits the required input-output tests. The goal of this dissertation is to move nonlinear multivariable control closer to industrial application by addressing practical model identification questions. Results from three test cases are presented and discussed. The results have shown a decrease in parameters of as much as 99% without a significant loss in model fidelity. TABLE OF CONTENTS ABSTRACT ................................................................................................................................................. iii LIST OF FIGURES ...................................................................................................................................... vi NOMENCLATURE ..................................................................................................................................... ix ACKNOWLEDGEMENTS.......................................................................................................................... xi INTRODUCTION ......................................................................................................................................... 1 Chapters I BACKGROUND....................................................................................................................................... 5 Linear Multivariable Model Predictive Control ......................................................................................... 7 Volterra Series ......................................................................................................................................... 14 Volterra Nonlinear Model Predictive Control .......................................................................................... 16 Laguerre Polynomials ............................................................................................................................... 18 II LINEAR IDENTIFICATION USING VOLTERRA-LAGUERRE REDUCTIONS ............................ 21 III NONLINEAR IDENTIFICATION USING VOLTERRA-LAGUERRE REDUCTIONS .................. 32 Testing Procedure Issues .......................................................................................................................... 37 Input Sequence Form ............................................................................................................................... 38 Identification of SISO Second Order System ........................................................................................... 40 IV CASE STUDY I: ISOTHERMAL FREE RADICAL POLYMERIZATION CSTR ........................... 46 Linear Volterra-Laguerre .......................................................................................................................... 47 Second Order Volterra-Laguerre .............................................................................................................. 56 V CASE STUDY II: VAN DER VUSSE REACTION ............................................................................. 64 Linear Volterra-Laguerre .......................................................................................................................... 68 Second Order Volterra-Laguerre .............................................................................................................. 75 VI CASE STUDY III: IDENTIFICATION OF EXPERIMENTAL PH DATA ...................................... 83 Procedure ................................................................................................................................................. 85 Determination of Process Test Conditions ............................................................................................... 87 Process Step Testing ................................................................................................................................ 88 Data Conditioning ................................................................................................................................... 90 Mean Center and Scale the Input and Output Data .................................................................................. 92 Model Identification ................................................................................................................................ 93 v CONCLUSIONS AND RECOMMENDATIONS ....................................................................................... 97 Recommendations for Further Work ........................................................................................................ 99 APPENDIX EXTENSION OF EQUATIONS TO MULTI-INPUT SINGLE OUTPUT SYSTEMS ......101 REFERENCES ...........................................................................................................................................104 LIST OF FIGURES Figure Page 1. Block diagram of multivariable IMC ................................................................................................ 6 2. Volterra-Laguerre fits for identifications based on 1st order Volterra kernel .................................. 28 3. Laguerre approximation of Volterra system .................................................................................... 30 4. Linear Volterra-Laguerre approximation of process with 3 of 15 intervals deadtime ..................... 30 5. Typical input step sequence - 100 Steps ......................................................................................... 39 6. Volterra 2nd order process - R = 2 (2 linear coefficients and 3 2nd order coefficients.) ................... 42 7. Volterra 2nd order process - R=3 case (3 linear coefficients 6 2nd order coefficients) ..................... 42 8. Volterra 2nd order process - R=4 case (4 linear coefficients 10 2nd order coefficients) ................... 43 9. 2nd order Volterra process step test data with noise ......................................................................... 43 10. Volterra 2nd order process ID with noise - R=2 case ....................................................................... 44 11. Volterra 2nd order process ID with noise - R=3 case ....................................................................... 45 12. Volterra 2nd order process ID with noise - R=4 case ....................................................................... 45 13. Hardware configuration for polymerization CSTR ......................................................................... 49 14. Step response curves CSTR ............................................................................................................ 49 15. Improvement of linear Volterra-Laguerre ....................................................................................... 51 16. Standard deviation based on # of steps ........................................................................................... 52 17. Process data for 20 step identification data set with noise .............................................................. 53 18. Volterra-Laguerre linear improvement for process with noise ........................................................ 53 19. Volterra and Volterra-Laguerre linear predictions on process with noise ....................................... 54 20. Linear Volterra-Laguerre dependence on time scale factor ............................................................ 55 21. Volterra-Laguerre 2nd order error reduction .................................................................................... 57 22. Volterra-Laguerre 2nd order error reduction .................................................................................... 58 vii 23. Standard deviations of 2nd order identification ................................................................................ 59 24. 2nd order improvement for process with noise ................................................................................. 59 25. 2nd order standard deviation for identification of process with noise .............................................. 60 26. 2nd order predictions on 20 step data ............................................................................................... 61 27. Effect of time scale factor 2nd order objective function ................................................................... 62 28. Hardware configuration of Van Der Vusse isothermal CSTR reactor ............................................ 65 29. Step response curves for .5 normalized feed steps .......................................................................... 67 30. Process loci for Van Der Vusse reactor ........................................................................................... 68 31. Model comparison for steps on Van Der Vusse system .................................................................. 69 32. Improvement of linear Volterra-Laguerre ....................................................................................... 71 33. Standard deviation of linear Volterra-Laguerre IDs ........................................................................ 71 34. Typical 20 step Van Der Vusse ID set with noise ........................................................................... 72 35. Improvement of linear Volterra-Laguerre on process with noise .................................................... 72 36. Standard deviation of linear Volterra-Laguerre IDs with noise ...................................................... 73 37. Linear identification comparison on Van Der Vusse data with noise ............................................. 74 38. Time scale factor "a" effect on linear optimization ......................................................................... 74 39. Volterra-Laguerre 2nd order error reduction .................................................................................... 77 40. Volterra-Laguerre standard deviation of 2nd order ID ..................................................................... 77 41. Model prediction of noise free Van Der Vusse system ................................................................... 78 42. Comparison of steady state gain loci for process and discrete models ............................................ 78 43. Volterra-Laguerre 2nd order error reduction- process with noise .................................................... 80 44. Volterra-Laguerre standard deviation of 2nd order ID - process with noise .................................... 80 45. Typical 30 step model prediction of Van Der Vusse system with noise ......................................... 81 46. Effect of time scale factor "a" on 2nd order optimization ................................................................ 82 47. PH control experimental apparatus ................................................................................................. 83 48. PH step test raw data ....................................................................................................................... 91 49. Mean centered and detrended PH data ............................................................................................ 94 viii 50. PH process ID Volterra-Laguerre improvement ............................................................................. 96 51. PH identification fit with Laguerre order of four ............................................................................ 96 NOMENCLATURE a Time scale factor for Laguerre expansion d Disturbance variable value e Difference between desired trajectory and actual trajectory gi i'th order impulse coefficients m Time intervals to horizon of fading system p Time horizon of fading system u input (independent) variable value hi i'th order Volterra kernel N Order of the Volterra series r Order of the Laguerre approximation li Continuous Laguerre operators Li Discrete Laguerre operators t Time T Time step size y Control variable value y* Desired control variable value or control trajectory y ~ Estimated control variable value G ~ Impulse coefficient matrix H ~ Impulse coefficient matrix for prediction of past effects P General nonlinear operator d Unmeasured disturbance vector x u Vector of future manipulated variable moves  MV move penalty matrix  Delta time  Real Laguerre coefficients, Combined impulse coefficient matrix  Combined impulse coefficient matrix for past moves, Laguerre model matrix ref  Reference trajectory  CV weighting matrix ACKNOWLEDGEMENTS I would like to acknowledge and give thanks to Dr. Philip Smith for his never-ending support and mentoring in school and in life. I would also like to acknowledge the support of my long-suffering wife, Catherine. INTRODUCTION The petroleum and petrochemical industries have long recognized the value of good process control in solving many of the problems facing the daily operation of process plants. Processes are complex and interactive with many manipulated variables and competing control objectives. Advanced computer control is used to stabilize these processes' operation, thereby increasing efficiency while minimizing pollution and maximizing profit. Current industrial state of the art techniques utilize linear multivariable process control (MPC) in conjunction with unit optimizers to achieve these goals. Unfortunately, many of the most difficult processes are not linear. As a result, the attempt is made to apply the control over a very narrow operating region where the process can be treated as linear. Unit optimizers often move the process out of the region where the models were originally identified. These problems have spurred the need for improved control techniques that can be easily implemented in industry for complex nonlinear processes. Linear multivariable control techniques have been widely accepted in industry as a tool to control important constrained, interactive systems. The appeal of these control techniques throughout industry has stemmed from their ease of use and ability to handle interactive constrained systems. Linear models are easily identified with process step tests and linear multivariable control problems are easily solved using linear matrix math. There is no requirement to have even a simple fundamental model of the process. Thus, linear MPC has become a standard general purpose tool used throughout industry. There are several commercial multivariable controllers available (DMCplus - Aspentech, Profit Plus : RMPCT - Honeywell, and others). These controllers are internal model-based and 2 provide significant benefit over traditional single input single output (SISO) techniques. The tools for implementing these controllers are well developed and make implementation straight-forward. These controllers all use linear models to represent the process being controlled. The linear models for control are obtained from step tests and input/output (I/O) data from the process. These linear controllers work well for many processes over a limited range. Their best performance is around the point of process testing and identification. As the nonlinearity of the process increases, the range of acceptable control using these techniques decreases. In addition, due to process constraints, the final operating point of these controls is often at a different point than where testing was initially done. The result is significant mismatch between the model and the process at the true operating point of the controller. Model and process mismatch can also be caused by the use of on-line plant optimizers. Optimizers (depending on plant economics) often drive individual controllers to operate in regions different than the point of original process testing. Since linear models are identified at one point, this also may cause significant mismatch between the model and the process, ultimately contributing to poor control. Nonlinear controllers can fit nonlinear processes over a much broader range, thus solving many of the control problems associated with linear controllers. However, nonlinear controllers have traditionally had problems of their own. The mathematics involved in moving from linear to nonlinear models is significant. Unlike linear techniques, nonlinear process models are difficult to invert (an important step in implementation of control) and invertability is not usually guaranteed. Identification of nonlinear models from plant I/O data has also proven to be very difficult. Some success with nonlinear control has been obtained where first principle models (nonlinear differential equations) are available, but these cases are specific to individual pieces of equipment and thus are not attractive as a general tool (Gerstle et al., 1991, Zafiriou et al., 1991). What is needed in industry is a general form nonlinear controller that could extend the range of control of nonlinear processes significantly over that of linear MPC algorithms. In order to be 3 practical and general, the identification of the nonlinear model parameters would need to be obtained from process experimental testing and I/O data over a suitable operating range. This dissertation presents an extension of work done in the area of nonlinear control using Volterra models (Al-Baiyat et al., 1986, Doyle et al., 1995, Maner et al., 1994, and Pearson et al., 1992). The Volterra approach has the advantage of using a linear formulation with iterative corrections for higher order nonlinearities. It allows for the use of the linear matrix inversion techniques that are well developed and understood while accounting for process nonlinearities. A key step in the implementation of any model-based control system is process identification. As a general tool used in industry, process parameters for these models must be obtained from plant I/O data. In the past, the number of parameters necessary to describe nonlinear Volterra systems has been prohibitive. This dissertation also builds upon the work of Zheng and Zafirou (1995) where orthogonal Laguerre functions are used to reduce the number of coefficients required for nonlinear models to a feasible number. The significant contributions of the work in this dissertation are in the area of identification of industrial nonlinear Volterra models from input-output data. Until now, the application of nonlinear multivariable control using higher order Volterra models has been limited to cases where the Volterra model could be identified from fundamental physics-based models or on systems that could be tested ad infinitum to determine higher order Volterra coefficients. Using the approach described in this dissertation, nonlinear models end up in the form of second order Volterra models that can be used in conjunction with work done by researchers such as Doyle et al. (1995), Maner et al. (1994), and Pearson (1993). The result of combining the work done in this dissertation and the work done by the researchers above is a generic nonlinear multivariable controller identified from process step testing and input-output data. This project provides a methodology for the practical implementation of second order nonlinear multivariable control using the existing framework used by industry in linear multivariable control. 4 A significant benefit demonstrated in this work is the ability to identify models using I/O data for a broad range of nonlinear processes. Difficult process characteristics such as deadtime, inverse response, nonlinear dynamics, and changing gains including gain sign changes have been demonstrated. One goal of this dissertation is to move this nonlinear multivariable control technology toward industrial application by addressing practical model identification questions. Such questions include identification of second order models from I/O data, experimental test requirements for model identification, and nonlinear parameter estimation. Since most processes in industry are nonlinear to some degree, it is expected that this nonlinear controller will significantly improve control in many industrial applications. Model improvements eliminating up to 90% of the error from linear models have been demonstrated for systems tested in this work. In addition, models for processes with gain sign changes have been demonstrated. CHAPTER I BACKGROUND Model predictive control (MPC) is essentially defined as any control algorithm that uses models to predict the future behavior of a process and generates control moves based on that prediction to control the future behavior of the process. Linear model predictive control has been widely used in industry for some time. The appeal of these control techniques on industrial processes has stemmed from their ease of use and ability to handle interactive constrained systems. Linear models are simple to identify and linear multivariable control problems are easily solved using linear matrix math. The first industrial applications incorporated two separate but similar formulations. Dynamic matrix control (DMC) was presented by Cutler and Ramaker (1980) and Prett and Gillette (1980). Model algorithmic control (MAC) was first presented by Richalet et al. (1978). Garcia and Morari (1982) first pointed out the similarities of DMC and MAC approaches. Garcia and Morari also grouped these methods into a formulation they called internal model control (IMC). A block diagram of the IMC formulation is shown in Figure 1. These approaches have found widespread industrial acceptance (Cutler and Yocum, 1991; Gerstle and Hokanson, 1991; Grosdidier et al., 1992; Richalet et al., 1978) and are the basis for industrial multivariable control packages available today (DMCplus - Aspentech, Profit Plus : RMPCT - Honeywell, and others). It can be seen from the block diagram in Figure 1 that with model predictive control, the model is in parallel with the process. The controller is essentially a feedforward controller, with only the mismatch between the model and the process in the feedback loop. In addition, 6 Gc Gp G ~ + - ++ -+ d ysp y controller process model u Figure 1 Block diagram of multivariable IMC convolution models are used and the result is a controller with several appealing characteristics, including the following:  No parametric form of the process model is assumed. Therefore, little prior process knowledge is needed before implementing control. In addition, control algorithm application is independent of the process order. Thus, high order and complex processes where little is known about the physical laws governing the system can be modeled and controlled using the same procedures as those used for low order processes.  Linear process models for these controllers can be obtained from time series analysis of I/O process data that do not require a significant fundamental modeling.  The control approach handles nonminimum phase process behaviors such as dead time and inverse responses transparently.  The control takes care of the interactive nature of multiple input multiple output (MIMO) systems optimally.  Constraints on independent and dependent variables can be included in the control objectives.  The control performance is often very robust. It has been shown that model gain variations of greater than four can be tolerated with only small changes in performance 7 (Richalet et al., 1978). Other researchers have discussed robustness on a more general basis ( Morari and Zafiriou, 1989; Muske and Rawlings, 1993). Because of the many advantages of this control approach, research has been done on extending this formulation to the nonlinear realm. One promising approach for this extension of the IMC formulation uses Volterra series. Volterra series are power series first studied by the mathematician Vito Volterra in the 1880s (Schetzen, 1980). The Volterra series approach essentially extends the convolution integral approach currently used in many industrial linear MPC packages to include higher order models. In recent years, there has been substantial research on the use of these higher order Volterra series in control (Doyle et al., 1995; Maner, 1993; Maner et al., 1994; Pearson, 1992) Linear Multivariable Model Predictive Control Linear MPC will be reviewed here, as the proposed nonlinear MPC is an extension of and has its basis in work done previously in linear MPC. There are several good review articles on linear multivariable control (Garcia et al., 1989; Marchetti et al., 1983; Ricker, 1985; and Rouhani and Mehra, 1982). Below is a brief review of the MAC formulation theory. Linear continuous processes have been modeled using equations of the following form: y t  g  ut  d    1 0 1 ( ) (1.1) This equation is the linear convolution integral and also represents the continuous first order Volterra kernel (Schetzen, 1980). If the system is causal, time invariant, and has fading memory, the equation can be approximated in the discrete time domain by the finite sum:        m i i y t g u t i 1 (1.2) 8 The gi 's represent the impulse response coefficients of the process. yt is the predicted output at time t and ut  i is the input to the plant at time t  i. The series has been truncated at m terms because the system has fading memory. This means that effects caused by changes in independent variables at some point in time fade to approximately zero and thus, the process can be approximated with a finite series. The time horizon can be defined mathematically as mT  p where p is the time horizon of the process. This represents the time for input effects to fade to approximately zero. For systems that do not fade, such as integrating processes, this truncated formulation will not work. Equation (1.2) is the fundamental prediction equation for the MAC control technique discussed earlier, and also represents the first order kernel of the Volterra series. In matrix form, equation (1.2) can be rewritten as: y G~u m (1.3) where, G ~ is the first order impulse coefficient matrix model approximating the real process G. m y is the model predicted output over one time horizon, and uis the vector of inputs over one time horizon. For the SISO case, the overall prediction equation can be written as: y Gu  Hu  d m past ~ ~ (1.4) which states that the predicted future process output over one time horizon m y is calculated based on the effect of future moves, u G ~ term, plus the effect of past moves of both manipulated and measured disturbance variables, past u H ~ term, plus an unmeasured disturbance vector d. The G ~ and H ~ are dynamic impulse matrices built from impulse coefficients as: 9 1 2 1 3 2 1 1 2 1 0 0 0 0 0 0 m m m g g g g g g g g g g                 G           (1.5) 2 3 4 2 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m m m m m m m g g g g g g g g g g                      H             (1.6) In the SISO case, 􀝑 and 􀢛􀢖􀢇􀢙􀢚 are column vectors containing one time horizon 􀝌 of future and past independent variable positions, respectively. The vector 􀢊 is a column vector of expected future disturbances. In most circumstances, 􀢊 is a constant vector containing the most recent mismatch between the model and actual plant output. This assumes that the unmeasured disturbance will stay constant and the best estimate of the disturbance is the current value. It is through the unmeasured disturbance vector that the feedback action of the controller is implemented. In some cases, a filter on this feedback term called a "bias filter" is added to reduce noise and help stabilize the controller. In order to achieve stable control, it is sometimes desirable to have the process follow some trajectory y * which is usually defined as the filtered setpoint or deviation that would return to the setpoint sp y . If you assume a first order time constant for the closed loop reference trajectory, ref  , the setpoint trajectory can be calculated as: 10 sp y e y k e y k ref T ref T                       1 ( ) 1 1 *( 1) 1 1 * (1.7) where sp y is the desired setpoint value and T is the sample time. In this equation, ref  becomes one of the primary tuning parameters of the controller. Without ref  , or some other way to dampen the control action, the controller would be too aggressive and very sensitive to model mismatch. Other methods for dampening the control action are shown by Ricker (1985). One of these methods places a penalty on manipulated variable movement. This technique will be discussed later. With y * defined, an error vector can be represented as: err  y Hu d past * ~ (1.8) Here, err represents the error between the desired trajectory and the actual trajectory if no further control moves are implemented. In unconstrained control, the final step is to calculate the optimal future control moves opt u that will minimize the error. The least squares solution to this problem is: u GTG GT err opt ~ ~ 1 ~  (1.9) For the multivariable case, the equations are expanded as follows:            l l c cm m G G G G ~ ~ ~ ~ 1 11 1      (1.10) 11         ~ ~ ~ ~ H H H H 11 1 1      I c cI (1.11) where ml is the number of manipulated variables, c is the number of control variables, and I is the number of independent variables (manipulated plus measured disturbance). The dynamic impulse coefficient matrix,  , is used for future move calculations and  is the dynamic matrix used to calculate control variable trajectories based on past independent variable moves. The MIMO control move equation then becomes: err  y  u  d past * (1.12) u    err T  T opt    1 (1.13) where  is a diagonal weighting matrix.  then becomes another tuning parameter that is used to scale control variables based on units and relative importance of error. This is of particular value in the case of an over-constrained system. In that case, it may not be possible to obtain the desired set points of all control variables (CVs). The relative magnitude of CV's error will be set by the weighting factor. In the under-constrained case, equation (1.13) has an infinite number of solutions. In practice, another term is often included which will drive the manipulated variables to ideal resting values or minimize their movement. In this case, the problem to be solved is: err err u u  u T T  opt    2 max 1 (1.14) 12 where  is a diagonal matrix containing penalties for movement of the manipulated variables from their nominal values. The least squares solution of this problem is: u    I err T  T opt     1 (1.15) Techniques for adding constraints to MIMO systems have been described by Hanson (1986) and Ricker (1985). If at each sampling time we solve equation (1.13) subject to the constraints: L U u  u  u (1.16) L U y  y  y (1.17) where U u and L u are the manipulated variable upper and lower constraints, respectively. If we define a new vector: L u*  u  u (1.18) we can rewrite the independent variable constraints as: u*  0 (1.19) U L u*  u  u (1.20) In the same fashion, we can write the control variable constraints as: 13  u  y  u  u d L L past  *   (1.21) u  y  u  u d U L past  *   (1.22) These constraints are in the form: Cu*  r (1.23) where               I C (1.24) and                      y u u d y u u d u u r U L past L L past U L     (1.25) Using these transformations, Maner (1993) has shown the new optimization equation as:  * *  2 * 1 * max a u u Bu u T  T (1.26) where a y u d uTB L T past T  *     (1.27) B   T    (1.28) 14 which has the unconstrained least squares solution: u*  B1a (1.29) Other simplifying steps are also used to simplify the calculation and increase the controllers robustness. Ricker (1985) describes a technique called blocking, which reduces the number of control moves calculated while also reducing the size of the matrices in the calculations. In this discussion, the basic formulation for the constrained MAC controller has been developed. There are some simple but important differences in the MAC controller and the DMC control strategies. First, in industrial applications of the DMC controller, the desired trajectory of the controller, y * , is always the current setpoint of the controller. There is no first order filter applied to the setpoint. The controller is suppressed and stabilized with a diagonal MV weighting factor matrix using a move suppression factor 2. Second, the DMC controller uses step coefficients rather than impulse coefficients in the dynamic matrix model of the plant. As a consequence, the control moves are calculated as change from the current position rather than absolute position. Volterra Series Volterra series are power series first studied by the mathematician Vito Volterra (Schetzen, 1980). Volterra series include time invariant series that take the form: 15   1    1  1 2  1 2   1   2  1 2 y t   g  u t  d    g  , u t  u t  d d                             3 1 2 3 1 2 3 1 2 3 g  , , u t  u t  u t  d d d       n n n n   g  1   u t  1 u t  d 1d   , ,        (1.30) If the system is causal with fading memory, the equation (1.30) can be discretized and shown as: y t  y t   y t   y t    1 2 3 ( ) (1.31) where        m i y t g i u t i 1 1 1 ( ) (1.32)            m i m j y t g i j u t i u t j 1 1 2 2 ( ) ( , ) (1.33)                 m i m j m k y t g i j k u t i u t j u t k 1 1 1 3 3 ( ) ( , , ) (1.34)             m j m j m j n l n n n n y t g j j u t l 1 1 1 1 1 1 2 ( )  ( , ) (1.35) In this form, y (t) i is the i'th order Volterra kernel. Thus, ( ) 1 y t is the first order Volterra kernel and is equivalent to the linear impulse response. i g is the i'th order operator which contains m i coefficients. Therefore, if there is a one-by-one process model with a time horizon of 60 coefficients and a second order Volterra model is used, there are 1860 impulse coefficients required. For a third order Volterra model, 217,860 coefficients would be required. The main deterrent to using Volterra models for modeling nonlinear processes in industrial applications is 16 the large number of coefficients required. Nonetheless, there has been significant work done in developing a nonlinear control technique using these models. Volterra Nonlinear Model Predictive Control There has been work done by Doyle et al. (1994) and Maner et al. (1996), using Volterra operators as a means of separating a nonlinear control problem into a linear control calculation and an auxiliary loop, which augments a nonlinear correction term. The major advantage to this approach is that the need for a nonlinear inverse is eliminated. This removes a major stumbling block in nonlinear control because inversion of the nonlinear model is very difficult and even impossible in many cases. In addition, the nonlinear correction is a straightforward extension of the linear model predictive control currently being used in industry. The controller is constructed based on the linear inverse with nonlinear corrections that do not require inversion of the nonlinear portion of the model. The approach is outlined below. In a process that can be represented by a nonlinear I/O correlation: ~y  Pu (1.36) where P is a nonlinear operator that operates on the plant inputs u. If the operator can be broken into a linear and nonlinear portion: P  L  N (1.37) then P  L(I  L1N) (1.38) and 17 P 1  (I  L1N) 1 L1 (1.39) This separation assumes that the linear inverse L-1 and the inverse of (I+L-1N) exist. Both the Volterra and Volterra-Laguerre models discussed earlier fit into the category of the overall nonlinear model being represented by a linear plus nonlinear operator. As stated earlier, the linear operator of the Volterra model is the same form used in some of the commercial linear controllers. If we substitute equation (1.36) into (1.31) we get: n P  P1  P2  P3    P (1.40) where the n'th order Volterra or Volterra-Laguerre expansion is represented as Pn. This now fits the form of equation (1.37) with P1=L and N represents the sum of all nonlinear operators. If y represents the actual output and y* is the desired trajectory, the error can be defined as: e  y*  y (1.41) In this case, the error e contains the effects caused by the combination of model mismatch and unmeasured disturbances. Combining equation (1.41) and equation (1.36) gives: y  P[u]  e (1.42) The desired control action can be determined by the solution of: y y u min  *  (1.43) 18 which can be shown as:  1 y*  e opt u P (1.44) if   y*  e (1.45) Substituting equations (1.44) and (1.45) into equation (1.38) the overall control equation for the nonlinear case becomes:   1  1 1 1 3 1 1 2 1 1 1 u  I  P P  P P   P P  P opt  n (1.46) Once again, this solution assumes that the linear inverse and the inverse of the composite term can be obtained. Thus far, it has been assumed that the matrices P, L, and N are square. However, in a typical multivariable problem, the matrices are rarely square. Doyle et al. (1994) showed that generalized inverses can be used instead of the actual inverse, yielding the same result. Nonsquare systems are therefore not a problem. Equation (1.46) represents the optimum control for a square system with no manipulated variable or control variable constraints. In the nonsquare system, all inverses are replaced with pseudo inverses and the result is the least-squares optimum of the control problem. Laguerre Polynomials There has been work done with Laguerre functions (Dumont and Fu, 1993; and Zheng and Zafiriou, 1995) in an effort to reduce the number of coefficients required to obtain a Volterra representation. Laguerre functions are orthonormal functions that can be used to approximate 19 Volterra series. The advantage of the Laguerre functions is that they are orthogonal and all coefficients of the discrete series are independent of all other coefficients. As a result, a selection process can be used to eliminate coefficients with small contributions to the model without the need to recalculate the remaining coefficients. Thus, the number of coefficients required to approximate a nonlinear Volterra system can be greatly reduced. Laguerre series are obtained by orthonormalizing the power series (Lee, 1960): at n e  at for n  0,1,2 The resulting form is: at o l (t)  2ae (1.47) l t   2a2at  1e  at 1 (1.48) l t   2 a2a 2 t 2  4 at  1e  at 2 (1.49)                           2 2 1 1 2! 2 ! 1 2 1 ! 2 ! 2 2 n n n n n n n t n t n n a n t n a n l t a a (1.50) This form is orthonormal over the domain [0, ) . The time scale factor "a" is a tunable parameter that should be set near the dominate root of the system (Wahlberg, 1991; Zheng and Zafiriou, 1995). The time scale factor will be discussed in more detail later. A sampled data system can be expressed as: 20         i 1 i i f k  l k (1.51) where  i are the real Laguerre coefficients calculated as:         k 0 i i  f k l k (1.52) In practice, f k  can be approximated using a finite series of Laguerre functions. In the case of the Volterra-Laguerre approximation, the Laguerre equations are used to approximate the Volterra kernel. The problem in identifying these are that the Volterra kernel is not known and identifying the Volterra kernel for higher order systems is not practical because there are too many coefficients to be identified. If it were practical to identify the Volterra kernel directly, Laguerre reductions would not be necessary. Thus, the approach used here is to use Laguerre functions and input-output data to approximate the Volterra kernel for the system. CHAPTER II LINEAR IDENTIFICATION USING VOLTERRA-LAGUERRE REDUCTIONS The Identification step is usually the most difficult step in implementing linear MPC. Identification is even more difficult in nonlinear processes. As stated previously, one of the main advantages of using the Volterra-Laguerre model is that the number of coefficients required to describe the system is greatly reduced from the Volterra model alone. Work has been done on the identification problem for the Volterra model by Nowak et al. (1994), Pearson et al. (1992 and 1993), and Zheng et al. (1994). Because of the number of coefficients required, none of the papers suggest that it is practical to attempt identifying high order Volterra models from I/O data and all of the papers reviewed limit their approach to second order Volterra models. For Volterra-Laguerre models, there has been work done by Dumont et al. (1991 and 1993), Zervos et al. (1988), and Zheng et al. (1995). Zheng demonstrated the advantage of reduction of coefficients in identification and discussed control relevant identification techniques. The advantage of using Laguerre series comes from the reduction of parameters necessary to represent the system. An additional advantage is the theoretical ease with which the orthogonal coefficients can be determined. This discussion will be limited to Laguerre series which are complete in L2[0,) space. It is also assumed that the process being modeled, the impulse coefficient model g , is L2 stable, i.e.: 22       0 g 2 t dt (2.1) This is the case for fading memory, self regulating systems. Based on this, the Riesz-Fisher theorem (Dumont et al., 1991) states that the process (in this case linear Volterra model) can be approximated as:       1 ˆ i i i g  l k (2.2) and for any positive,  there is an n such that:                0 2 1 g t  l k dt  n i i i (2.3) where 's are the Laguerre coefficients, gˆ is the approximated Volterra model, and  represents the acceptable integrated error. Dumont et al. (1991) and Zheng et al. (1995) proposed two ways to determine the coefficients . The first is called the correlation method and is a natural consequence of orthogonality. In this case, the i'th Laguerre coefficient can be calculated as: gkTl kT n i n k i  1 1 ˆ (2.4) This equation assumes that the process input is Gaussian white noise. In the technique being proposed, the process being represented is the first order Volterra model so an input of Gaussian noise has little meaning. Also assumed in this technique is that there are no 23 unmeasured disturbances or any other correlated or uncorrelated noise, and that the time scale factor in the Laguerre equations is already known. The Laguerre time scale factor will be discussed later. This approach makes little sense for approximation of the Volterra kernel. The second approach presented by Dumont (1991) is the least-squares approach. In this approach, the Volterra kernel is represented as:     W u Tu gˆ T  (2.5) where the first term on the right hand side represents the model, the second term represents a similar Laguerre expansion for the unmodeled dynamics, and the third represents the effect of random noise on the process. Since there is no structure known about the unmodeled dynamics or the random noise, the best that can be done is assume that these will be zero in the mean. Because of this, Dumont goes on to determine that the expectation of the second and third terms of equation (2.5) go to zero. In the above equation, the parameters are defined as:  ( ), , ( ) 1 l kt l kT     n (2.6)  ( ), , ( ) 1 l kT l kT n T  u   (2.7)  T  (T ), , (nT ) (2.8)  (T ), , (nT ) u u Tu      (2.9) gˆ T gˆ (T ),, gˆ (nT ) (2.10) The least squares estimate of  becomes:   g   1 ˆ (2.11) 24 In the previous section, it was assumed that the Volterra kernel is already known. This is not the case in industrial processes. In the Volterra-Laguerre form, the above technique is used to approximate the Volterra-Laguerre kernel from I/O data. Zheng et al. (1995) presented the correlation method approach to estimate the Volterra-Laguerre kernel. Zheng did not discuss the requirements of the identification based on I/O data in his paper. In addition, Zheng assumed the time scale factor in the Laguerre equations had already been determined, making the problem linear in parameter identification. The input excitation signal to the process has been assumed to be Gaussian white noise to this point. However, within industrial processing plants, normally Gaussian random noise inputs are not practical and not allowed. The time scale factor and the input step requirements will be discussed in more detail later. The first step in the evaluation of Laguerre reductions is to test Laguerre functional approximations for Volterra kernels of linear systems. g   T (2.12) g is a vector with the Volterra first order model for one time horizon into the future. Equation (2.12) shows the approximated Volterra kernel as opposed to the actual Volterra kernel. In the industrial process identification problem, the Volterra kernel is not known, and identifying it directly is impractical. Much of the derivation in the previous section assumed the true Volterra kernel was available. From this point on, g will represent the Volterra kernel identified from process I/O data based on the Volterra-Laguerre form.  ( ), , ( ) 0 l kt l kT    r (2.13)  T  (T ), , ( pT ) (2.14) 25 In this equation, p is the number of coefficients in one time horizon as defined by the Volterra model and r is the number of Laguerre terms used in the approximation. In expanded form:                                      l pT l pT l pT l pT l T l T l T l T l l l l r r r         0 1 2 0 1 2 0 1 2 0 0 0 0 (2.15) In Equation (2.15),  is a transformation matrix that maps the Laguerre coefficients to Volterra coefficients. The result is a reduction in the number of coefficients required to represent the system from p to r where p is the number of coefficients in the time horizon of the Volterra model and r is the order of the Laguerre approximation. Recognizing the first order Volterra model as:        p i i y t g u t i 1 (2.16) The prediction equation for one point becomes:     k y(k ) u (2.17) where u uk  uk  uk p k  ,  1 , ,  Equation (2.17) gives the prediction at the k'th interval. What is needed for optimal control is the control move calculation for one full time horizon into the future. In order to keep 26 the Volterra-Laguerre form, there must be a matrix of moves U. In addition, the future prediction of the process value must be broken into two predictions. One prediction is based on past moves and the other prediction is based on future moves as done in the linear MPC approach discussed in Chapter 1. The total prediction for one time horizon in the future becomes the sum of the predictions based on past independent variable moves and moves the for one time horizon in the immediate future. The independent variable move matrices and process predictions take the form:                     0 0 0 0 0 0 0 0 0 2 0 2 1        u u u u u u U p p p past (2.18)     past past Y U (2.19)                  1 2 1 1 2 1 1 0 0 0 0 0 0 u u u u u u u u U p p fut         (2.20)    fut fut Y U (2.21) To complete the identification of the Laguerre coefficients directly from I/O data it is necessary to group terms.   past past Y U (2.22) 27 let    past A U (2.23) which leads to: Y  A past (2.24) Given an input-output data set (u, y), the Laguerre coefficients  can be approximated by the least squares fit. If U is sufficient for the identification step, this should provide reasonably accurate predictions for any input U in the prediction equations. At this point, it is assumed that the time scale factor for the Laguerre equations is known. Figure 2 shows example plots of linear predictions based on a first order plus deadtime process model. The comparison is between an impulse coefficient model, Volterra-Laguerre where the Laguerre coefficients were approximated from the Volterra kernel, and direct calculation of the Laguerre coefficients from I/O data. The input for the plots has a step up at interval two and a step down at interval seven. The final plot has random input. As can be seen, a relatively few number of Laguerre terms produce reasonably good predictions. The exponential time scale factor used in the Laguerre equations was determined from a nonlinear least squares optimization. When the time scale factor is not optimal, the Laguerre model still converges to the actual, but more Laguerre terms are required to obtain an equivalent fit. The next step in this identification process is to use equation (2.24) with large amounts of I/O data to generate an accurate set of Laguerre coefficients. To accomplish the identification, equation (2.24) is used with the entire I/O data set. Ypast is the vector containing all of the output data. A is an nxr matrix, where n is the number of points in the data set and r is the order of the 28 Figure 2 Volterra-Laguerre fits for identifications based on 1st order Volterra kernel 0 5 10 15 -1 0 1 2 3 4 5 0 5 10 15 -1 0 1 2 3 4 5 0 5 10 15 -1 0 1 2 3 4 5 0 5 10 15 -1 0 1 2 3 4 5 0 5 10 15 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 0 5 10 15 -0.5 0 0.5 1 1.5 2 2.5 3 3.5  Impulse Coefficient Model * - Laguerre Coefficients from Volterra Kernel o - Laguerre Coefficients from I/O Data 7 Laguerre Terms 7 Laguerre Terms Random Input 4 Laguerre Terms 6 Laguerre Terms 5 Laguerre Terms Random Input 2 Laguerre Terms 29 Laguerre approximation. Note that A is dependent on Upast. It is reasonable to expect that the quality of the least squares fit is also dependent on the form of Upast. For this work, the input signal will be pseudo-random binary sequences (PRBS), based on a paper by Pearson et al. (1993). To illustrate this identification, normally distributed random noise was imposed on a 15 coefficient first order Volterra model to simulate the process. The process had no deadtime. A new random set of input data was generated and the identified Volterra-Laguerre model was compared to the simulated process. In Figure 3, the solid line represents the process, and the "+" marks illustrate the Volterra-Laguerre approximation. As can be seen, if the process has no dead time and first order dynamics, it does not take many Laguerre coefficients to represent the system. In Figure 3, the simulated process was changed to have 3 deadtime intervals out of 15. Again, the I/O data used for the identification had normally distributed random noise imposed upon it. As a result, more Laguerre terms are required to adequately approximate the process. The number of coefficients required to achieve reasonable results is still significantly less than that of the Volterra model. This illustrates the advantage and practicality of using Laguerre equations. The ability to approximate deadtime is based on the Laguerre equations similarity to Padé approximations for deadtime. Figure 3 and Figure 4 illustrate that the Laguerre approximation provides good results for linear systems. The ability to model deadtime is one of the most important requirements for models used in MPC applications in chemical industries. Processes without deadtime can often be controlled with simple feedforward schemes such as lead-lag and interactive moves. In addition, deadtime is a major contributor in the failure of simple PID feedback control. The fundamentals of data regression dictate that if fewer coefficients are used in the model, then something is sacrificed in the regression model results. With the Volterra model, there is no assumption about the dynamics or smoothness of the process. With the Laguerre approximation, there is some smoothing effect. This smoothing 30 0 50 100 150 200 -1 0 1 2 3 4 5 2 Laguerre Coefficients Squared Error = 2.583 Figure 3 Laguerre approximation of Volterra system 0 50 100 150 200 0 1 2 3 4 5 2 Laguerre Coefficients Squared Error = 91.63 0 50 100 150 200 -1 0 1 2 3 4 5 3 Laguerre Coefficients Squared Error = 20.39 0 50 100 150 200 -1 0 1 2 3 4 5 6 4 Laguerre Coefficients Squared Error = 7.916 0 50 100 150 200 -1 0 1 2 3 4 5 6 Laguerre Coefficients Squared Error = 3.773 Figure 4 Linear Volterra-Laguerre approximation of process with 3 of 15 intervals deadtime 31 works well if the process is relatively smooth or the randomness is caused by noise in the system but, if the process is repeatable and consistently erratic, Volterra-Laguerre equations may not be the optimal approach. CHAPTER III NONLINEAR IDENTIFICATION USING VOLTERRA-LAGUERRE REDUCTIONS The Identification step in implementing linear MPC is typically the most difficult step. Identification of nonlinear processes is even more difficult. This difficult identification is the defining factor as to why nonlinear MPC has not been widely implemented throughout industry. The number of step tests required for regression-based nonlinear models has been prohibitive. As stated earlier, one of the main advantages of using the Volterra-Laguerre model is that the number of coefficients required to describe the system is greatly reduced from the Volterra model alone. Although the results will not be discussed here, there has been work done on the identification problem for the Volterra model by Pearson et al. (1992 and 1993), Nowak et al. (1993 and 1994), and Zheng et al. (1994). None of these papers suggest that it is practical to attempt identifying high order Volterra models from I/O data because of the number of coefficients required. The work done in these papers has been limited to obtaining Volterra coefficients from first principle models for individual systems. There has been work done by Dumont et al. (1991 and 1993), Zervos et al. (1988), and Zheng et al. (1995) for Volterra-Laguerre models. Once again, these works use first principle models to obtain Volterra-Laguerre coefficients. The advantage of using the orthonormal Laguerre series comes from the reduction of parameters necessary to represent the system. An additional advantage is the theoretical ease with which the orthogonal coefficients can be determined. This work is limited to Laguerre series which are complete in L2 [0,) space. Laguerre series are used here because of their similarity to 33 Padé approximations of deadtime. This characteristic is important since deadtime is usually present and significant in chemical processes. It is also assumed that the process being modeled is L2 stable, i.e.:       0 y 2 t dt (3.1) This is the case for fading memory, self regulating systems. If this is the case, the Riesz-Fisher theorem (Dumont et al., 1991) states that the process can be modeled as:       i 1 i i y  l k (3.2) and for any positive  there is an n such that:                0 2 1 y t  l k dt  n i i i (3.3) where, i are the Laguerre coefficients and  represents the acceptable integrated error. Dumont et al. (1991) and Zheng et al. (1995) have proposed two ways to determine the coefficients . The first is called the correlation method and is a natural consequence of orthogonality. In this case, the i'th Laguerre gain can be calculated as: ykTl kT n i n k i   1 1 ˆ (3.4) This equation assumes that the input is Gaussian random noise. Also assumed is that all disturbances are measured disturbances and there is no other correlated or uncorrelated noise. 34 The second approach presented by Dumont (1991) is the least-squares approach. In this approach, y is represented as: Y      W u Tu T  (3.5) where the first term on the right-hand side represents the model, the second term represents a similar Laguerre expansion for the unmodeled dynamics, and the third represents the effect of random noise on the process. In the above equation, the parameters are defined as:  ( ), , ( ) 1 l kt l kT     n (3.6)  ( ), , ( ) 1 l kT l kT n T  u   (3.7)  T  (T ), , (nT ) (3.8)  (T ), , (nT ) u u Tu      (3.9) YT y(T ), , y(nT ) (3.10) The least squares estimate of  becomes:   Y   1 ˆ (3.11) Dumont continues on to demonstrate that if the expectation is taken, the second and third terms of equation (3.5) go to zero if the input and disturbances have zero mean. In the previous sections, the process Y was approximated using Laguerre polynomials. In fading memory systems, the finite Laguerre series can be used to approximate the Volterra kernels i g as shown in equations (1.31) through (1.34), the advantage again being that 35 the i'th order Volterra kernel i g has on the order of m i coefficients where m is the time horizon and i is the order of the Volterra kernel. The Laguerre approximation has approximately r i coefficients where r is the order of the Laguerre approximation. In practical cases, r will be much smaller than m. Zheng (1995) reported a reduction from 23,750 coefficients using the Volterra model to 20 coefficients using the Volterra-Laguerre model in one test problem. Utilizing this substitution, the Volterra-Laguerre expansion becomes:                 N n r j r j n l m i n n j l l n l l y k j j l i u k i 1 1 1 0 1 1 ~   ,, (3.12) where n is the order of the Volterra series and r and m are as previously defined. The same technique could be used to approximate the Volterra kernels. Zheng et al. (1995) presents the correlation method approach to estimate the Volterra kernel. It has been assumed that the input signal for this derivation is Gaussian random noise. Other authors, such as Pearson et al. (1993) suggest that a pure white noise input signal may not be required. If the system is causal with fading memory, the first two terms in Volterra approximation can be represented as: y t y t  y t  1 2 ( )   (3.13) where        p i y t g i u t i 0 1 1 ( ) (3.14)           p i p j y t g i j u t i u t j 0 0 2 2 ( ) ( , ) (3.15) 36 In this form, y (t) i is the i'th order Volterra kernel. Thus, ( ) 1 y t is the first order Volterra kernel and is equivalent to the impulse response form of the linear multivariable controller. i g is the i'th order operator which contains pi coefficients. The Laguerre approximation for the second order Volterra kernel is:       2 0 0 2 1 2 1 2 1 2 1 2 1  ,    j r j r j j j j g   l l    (3.16) defining              r rr r D        1 21 22 11 12 1    (3.17)          r LT l0 l1  l  (3.18) then       2 1 2 1 2 g  ,  LT  DL  (3.19) Equation (3.19) can be substituted into equation (3.15) to produce:                 p i p i y t u t i u t i LT t i DL t i 0 0 2 1 2 1 2 1 2 ( ) (3.20) If 37 U T u(t) u(t  T )  u(t  pT ) (3.21) and recall that                                      l pT l pT l pT l pT l T l T l T l T l l l l r r r         0 1 2 0 1 2 0 1 2 0 0 0 0 (3.22) then y (t)  U T D TU 2 (3.23) Equation (3.23) is the second order prediction for one point in time. The overall prediction at this time period is the sum of the linear and second order contributions as shown in equation (3.13). In this context, equation (3.23) gives the nonlinear correction for the linear model. Testing Procedure Issues In doing identification testing to obtain I/O data for nonlinear systems, there are several issues to be considered. The first concern is the input excitation required to adequately generate control variable response over the expected operating conditions and ranges. Input step size, step frequency, number of steps, duration, and technique all need to be considered. The number of steps required will be discussed later in this dissertation. 38 Input Sequence Form In theory, the ideal step sequence has been shown to be completely random both in frequency and magnitude (Zheng et al., 1995). Zheng showed that if the input signal is Gaussian random noise, for a long enough time, it is expected that all of the relevant frequencies of the process are excited and the regression models are optimized. Unfortunately, this is not practical in industrial applications. Plant personnel and management will not allow for random inputs, neither in size or duration, on live process variables. It could produce dangerous conditions, equipment damage, off-spec products, or destabilize whole process units. Other authors, such as Pearson et al. (1993) have suggested that a pure white noise input signal may not be required. Pearson suggests that elliptically distributed inputs may also give good results based on autocorrelation matrix analysis. These results give credence to existing step test techniques used in industry. The most commonly used input technique in industry for linear systems identification is pseudo-random binary sequences (PRBS). PRBS does a good job of exciting the process for linear systems but does not adequately excite higher order Volterra terms (Nowak et al., 1994). In linear identification, the step size is constant and the magnitude is chosen based on the minimum step size needed to excite the control variables in the process. This technique identifies only a linear process model around the testing point. It is assumed that the linear model represents the behavior of the process over the whole expected operating range of control. Nowak recommends, based on analysis of the correlation matrix, that for nonlinear models, pseudo-random multilevel sequences (PRMS) be used. For the work in this dissertation, an extension of the PRMS approach is used. Instead of using a fixed number of levels for input amplitude, input amplitude is varied in a pseudo-random fashion as well. A typical input sequence is shown in Figure 5. The amplitude of the steps are restricted between a minimum step size which is required to see movement in the control variables and a maximum step size that is considered safe. The steps are randomly varied around an average step size. The random moves 39 Figure 5 Typical input step sequence - 100 Steps follow a uniform distribution. The frequency is also randomly varied around an average step time. The average step time is set based on the dominant time constant of the process response and also has uniform distribution. This is necessary to insure adequate data in both the low and high frequency ranges. In industry, it is commonly considered unacceptable to randomly move process variables. This approach is more appealing to industry because the random nature of the moves is limited based on prior knowledge of the process. This also has the advantage of a pseudo-random amplitude component which is important for the nonlinear Volterra terms. This gives a uniform distribution of amplitudes in the testing range. This is new and an extension of the traditional PRMS which uses a fixed number of amplitudes in the step testing. 0 500 1000 1500 2000 2500 3000 3500 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Time Normalized Step 40 Identification of SISO Second Order System In the truncated Volterra approximation shown in equation (1.33), the identification of the second order Volterra-Laguerre approximation from I/O data consists of two steps. The first step requires performing a linear identification to get ( ) 1 y t . The residual of ( ) ( ) 1 y t  y t can then be used as the process data in the second order identification. Given equation (3.23), a solution must be found for the coefficient matrix D. Note that D has dimensions rxr and is symmetric; therefore, there are ( ) 2 1 r 2  r independent coefficients in the system. The identification problem can be posed as: ( ) ( ) ( ) min 2 2     NS t p y t TU T D TU (3.24) where NS is the number of I/O data points. Assuming there is an exact solution for ( ) 2 y t , the equation becomes: ( ) ( ) ( ) 2 y t   TU T D  TU (3.25) Multiplying both sides by ( TU ) and solving for D yields:       1 2 1 ( )( ) ( ) ( ) ( )( )   D  TU TU T TU y t TU T TU TU T (3.26) where the term  1 ( )( )  TU TU T is the inverse. This is the solution for only one data point. For the complete time horizon of the second order system, it becomes necessary to perform a nonlinear optimization. In doing so, advantage can be taken of D being symmetric. The total 41 number of independent coefficients become ( ) 2 1 r 2  r . Thus, equation (3.26) is used for the optimization with the objective function being equation (3.24). The first test of these equations is to identify simulated I/O data generated from a simulation using a second order Volterra model. The Volterra process was created using 15 linear coefficients and 120 second order Volterra coefficients. The process was essentially a first order plus deadtime linear response with second order dynamics added. The results of the identification are shown in Figure 6 thru Figure 8. The optimization was calculated using the nonlinear least squares algorithm in Matlab's (Mathworks, 1999) optimization toolbox. The errors are the least squares error between the process simulation and the approximation from the linear and second order Volterra-Laguerre models. The process was identified several times with increasing Laguerre orders until improvement diminished. In this case, a Laguerre order of 3 was completely adequate for modeling the process. This provided a reduction of coefficients from 135 Volterra coefficients to 9 Laguerre coefficients in the Volterra-Laguerre model. This illustrates the effectiveness of the Volterra-Laguerre model in representing a well-behaved second order Volterra system. The implications of this are that the work done using second order Volterra models to represent industrial processes can be used to set expectations of the Volterra- Laguerre models. With that said, this is the most optimistic case for the Volterra-Laguerre modeling. The data used for identification had no noise and the system was a well-behaved and well-defined Volterra model. To evaluate the robustness of the identification with noise present, noise was added to the second order Volterra process step data and the process model identification was performed again. The same process data as the previous test was used but normalized random noise, with a standard deviation equal to the average process (CV) movement, was added to the simulated I/O data before the identification. Figure 9 shows the process data with and without the noise. 42 Figure 6 Volterra 2nd order process - R = 2 (2 linear coefficients and 3 2nd order coefficients.) Figure 7 Volterra 2nd order process - R=3 case (3 linear coefficients 6 2nd order coefficients) 0 20 40 60 80 100 120 140 160 180 200 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 Time Data 1st Order 2nd Order LS Error 1st = 41.27 LS Error 2nd = 31.74 0 20 40 60 80 100 120 140 160 180 200 -0.4 -0.2 0 0.2 0.4 0.6 LS Error 1st = 41.2659 LS Error 2nd = 0.0031 Time Data 1st Order 2nd Order 43 Figure 8 Volterra 2nd order process - R=4 case (4 linear coefficients 10 2nd order coefficients) Figure 9 2nd order Volterra process step test data with noise 0 20 40 60 80 100 120 140 160 180 200 -0.4 -0.2 0 0.2 0.4 0.6 LS Error 1st = 41.2381 LS Error 2nd = 0.030297 Time Data 1st Order 2nd Order 0 20 40 60 80 100 120 140 160 180 200 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Time Process w Noise Actual Process 44 In the evaluation of model performance, the squared error between the Identified Volterra- Laguerre model and the process data without noise was used. The purpose of this is to illustrate the stability and noise rejection characteristics of the Volterra-Laguerre models. For this test, 75 steps were used in the process step test. The effect of the number of steps in the test will be discussed later. Figure 10 through Figure 12 show the predictions from the identified Volterra- Laguerre models and the process data without noise. As with the identification without noise, a Laguerre order of 3 produces a good fit with an error reduction over the linear model error of 94%. Figure 10 Volterra 2nd order process ID with noise - R=2 case 0 20 40 60 80 100 120 140 160 180 200 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 LS Error 1st = 8.8141 LS Error 2nd = 6.2614 Time Data LS 2nd Order 45 Figure 11 Volterra 2nd order process ID with noise - R=3 case Figure 12 Volterra 2nd order process ID with noise - R=4 case 0 20 40 60 80 100 120 140 160 180 200 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 LS Error 1st = 8.6235 LS Error 2nd = 0.54706 Time Data LS 2nd Order 0 20 40 60 80 100 120 140 160 180 200 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 LS Error 1st = 8.5617 LS Error 2nd = 0.70158 Time Data LS 2nd Order CHAPTER IV CASE STUDY I: ISOTHERMAL FREE RADICAL POLYMERIZATION CSTR The following example consists of an isothermal free radical polymerization of methyl methacrylate using azo-bis-isobutyronitrile as initiator and toluene as solvent (Congalidis et al., 1989; Doyle III et al., 1994; Maner et al., 1996). The manipulated variable is initiator flow rate and the control variable is number average molecular weight (NAMW). The first principles differential equations for the process model (Daoutidis et al., 1990; Maner et al., 1996) were obtained based on the assumptions of:  Perfect mixing  Constant heat capacity  Quasi steady state  No gel effect  Monomer inlet flow is constant   V F C C k k C P dt dC m m p f m o m in m ) (      (4.1) V F C FC k C dt dC I I I I I I in     (4.2)   V FD k k P k C P dt dD o T T o f m o o C d m  0.5  2   (4.3) 47   V M k k C P FD dt dD m p fm m o 1    1 (4.4) where 0.5 2 *         Td Tc I I o k k P f k C (4.5) The process model parameters and steady-state operating conditions are given in Table 1 and Table 2 as given in Maner et al. (1995). In Maner's work, Volterra models were obtained by performing Carleman linearization (Rough, 1981) of the first principle models. For this work, the same first principle models are used to generate I/O data that will then be used for identification. Figure 13 shows the hardware configuration of the CSTR reactor. Open loop step response curves are shown in Figure 14 for a .5 normalized step up and a step down of equal magnitude. As can be seen, the process is not symmetrical. This is a common nonlinearity observed in industry. This form of nonlinearity is seen in high purity distillation, combustion, process heaters, chemical reactors (reformers, polymers, alkylation, dewaxing, etc.), to name a few. Linear Volterra-Laguerre One of the results of this research is the insight gained into the value of using Laguerre reductions for identification of linear models. The most difficult step during the implementation of linear MPC on industrial processes is obtaining good models. In processes with many independent variables, the number of step tests required to obtain first order Volterra models can be significant. Thus, step testing often takes long periods of time, during which there are measurable disturbances to plant operations. If the process is noisy, more time and more steps become necessary. As the noise increases, the number of process step tests must also increase. 48 Table 1 Isothermal free radical polymerization parameters Parameters For Isothermal Polymerization CSTR Tc k Equilibrium Constant 103281 X 1010 m3/(kmol h) Td k Equilibrium Constant 10.930 X 1011 m3/(kmol h) I k Equilibrium Constant 1.0225 X 10- 1 m3/(kmol h) p k Equilibrium Constant 2.4952 X 106 m3/(kmol h) fm k Equilibrium Constant 2.4522 X 103 m3/(kmol h) f * Constant 0.58 F Monomer Inlet Flow 1.00 m3/h V Volume of Tank 0.1 m3 I in C Initiator Inlet Concentration 8.0 Kmol/m3 m M Monomer MW 100.12 kg/kmol min C Monomer Inlet Concentration 6.0 kmol/m3 Table 2 Polymerization CSTR steady state operating conditions Steady-State Operating Conditions For Isothermal Polymerization CSTR x1 = Cm = 5.506774 kmol/m3 x2 = CI = 0.132906 kmol/m3 x3 = Do = 0.0019752 kmol/m3 x4 = D1 = 49.38182 kg/m3 u = FI = 0.016783 m3/h y = 25,000.5 kg/kmol 49 Figure 13 Hardware configuration for polymerization CSTR Figure 14 Step response curves CSTR -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0 0.5 1 1.5 2 2.5 3 Normalized NAMW Time (TSS) 50 One of the advantages of using Volterra-Laguerre reductions is the reduced number of step tests required for identifying linear models. Since Laguerre equations are used to approximate the Volterra kernel in the Volterra-Laguerre model, it would be expected that Volterra models would always produce better results than the Volterra-Laguerre models in the identification. Later in this section, this is shown to be true for pure Volterra systems without noise. Figure 15 shows a comparison of the linear Volterra and the linear Volterra-Laguerre identification improvements for the CSTR process as a function of the number of steps in the tests and the Laguerre order. Figure 15 is based on the average of several step tests for both the Volterra and Volterra-Laguerre equations. The frame of reference for these plots is the linear Volterra model. The first order Volterra identification is the current standard used throughout industry. A value of 0% improvement means that the linear Volterra-Laguerre performed as well as the first order Volterra model for the system. A value of 100% indicates that the model matched the process data perfectly and 100% of the modeling error was eliminated. Thus, the 51 Figure 15 Improvement of linear Volterra-Laguerre units of identification performance are indicated as percent error reduction between the linear Volterra model and the process data. Figure 16 illustrates that as the number of steps in the step test increases, the advantage of using linear Volterra-Laguerre identification over the standard first order Volterra identification diminishes. Figure 16 also provides an indication of the number of steps required for the linear Volterra identification. For this case, it appears that 30 to 40 steps are required. This would be more steps than are usually performed on industrial applications for linear control models. The characteristics of the process itself are one of the factors that affect the number of steps required in the process identification test. The fact that this case represents a nonlinear process may increase the number of steps required. Nevertheless, with less than 40 steps and a Volterra-Laguerre order of 4 or higher, the Volterra-Laguerre still performs better than the straight Volterra identification. The second comparison for linear system identification concerns what happens when the process has noise. To illustrate this, normally distributed random noise was imposed on the -10.00 -5.00 0.00 5.00 10.00 15.00 20.00 0 10 20 30 40 50 60 AVE % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 52 Figure 16 Standard deviation based on # of steps process data and the system was identified again using the noisy signal. Figure 17 illustrates a typical 20 step response with the actual process signal and the same process signal with added noise. To test this system, the standard deviation of the noise was set as 1 times the standard deviation of the process movement for the data set. Because the real process was known, both the Volterra and Volterra-Laguerre identified models were compared to the actual process without noise. Figure 18 illustrates the relative improvement of the Volterra-Laguerre prediction as a function of Laguerre order and the number of steps used in the identification. Figure 19 illustrates predictions for a typical 20 step system identified with noise. It is clear from the results shown in Figure 18 that the Volterra-Laguerre model produces an improved fit over the Volterra model. The use of the Laguerre equations to predict the Volterra kernel acts as a process filter. 0.00 5.00 10.00 15.00 20.00 25.00 30.00 0 10 20 30 40 50 60 Standard Deviation of % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 53 Figure 17 Process data for 20 step identification data set with noise Figure 18 Volterra-Laguerre linear improvement for process with noise 0 1 2 3 4 5 6 7 8 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 Time (TSS) Normalized Concentration of B Process w Noise Actual Process 0.00 10.00 20.00 30.00 40.00 50.00 60.00 0 20 40 60 80 100 120 AVE % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 54 Figure 19 Volterra and Volterra-Laguerre linear predictions on process with noise As expected, as the number of steps increases, the Volterra-Laguerre model approaches the performance of Volterra. As stated previously, if the process (without noise) was not "smooth," the Volterra-Laguerre approximation would not work as well. It is expected that the same results could be obtained from a straight Volterra identification if a filter were used on the data before identification. The two significant advantages illustrated for the Volterra-Laguerre identification approach are the following:  The number of coefficients required to identify the system using Laguerre reductions is always fewer than what is required for the straight Volterra form.  The Volterra-Laguerre equations naturally smooth the data. Therefore, model quality is improved while the number of step tests required is reduced. As a result, for processes with long times to steady state, or many independent variables, the Volterra-Laguerre identification would be superior in industrial application. 0 50 100 150 200 250 300 350 400 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time (TSS) Normalized NAMW Linear Volterra Error 1st = 3.4448 Volterra- Laguerre Error = 3.0233 Data Linear Volterra-Laguerre Voterra 55 As discussed earlier, the Laguerre time scale factor should also be included in the parameter identification. The exponential time scale factor makes the identification well-behaved and nonlinear. Figure 20 illustrates the effect of the Laguerre time scale factor on the linear Volterra-Laguerre identification. This plot shows the importance of the time scale factor in the Identification. The Volterra-Laguerre identification can take either less or more CPU time than the Volterra alone because of the nonlinearity. The fact that there are fewer parameters being identified reduces the CPU requirements but the nonlinear time factor increases identification time. In this CSTR process identification, the CPU time requirements ranged from half the CPU time to 50% more CPU time for the linear Volterra-Laguerre identification than the straight first order Volterra identification. Because the development was done in Matlab (Mathworks Inc.), this may not be representative of compiled code. In any case, the identification is a one-time process and the CPU time would not prove to be a problem for industrial applications. Figure 20 Linear Volterra-Laguerre dependence on time scale factor 0 2 4 6 8 10 12 14 16 18 20 29 30 31 32 33 34 35 Time Scale a Squared Error R=2 R=3 R=4 R=5 R=6 R=7 56 Second Order Volterra-Laguerre The main thrust of this research is to identify second order Volterra-Laguerre coefficients from I/O data and quantify the potential improvement of these models over the linear models being used in industry today. A comparison of Volterra-Laguerre dynamic models obtained from I/O data to the first principle models to illustrate overall quality has also been performed. Both linear and nonlinear Volterra-Laguerre identifications are compared. Identifications were completed for several (10-15) step data sets to evaluate robustness and sensitivity of the identification procedure. Some of the topics to be discussed include the following:  the number of steps required for the linear and nonlinear identification.  the minimum order of the Laguerre approximation to obtain reasonable results.  the effect and selection of the Laguerre Exponential time scale factor.  the effect of process noise on control models. A key indication of the practicality of the identification procedure is the minimum number of independent variable steps required to obtain acceptable nonlinear models. The relationship between the number of steps and relative error reduction for this system is shown in Figure 21 and Figure 22. These figures show the average error reduction over linear Volterra models based on the identification of ten different input-output data sets for each corresponding number of steps. Once again, a value of 0% represents no improvement over linear Volterra models, a value of 100% represents a second order Volterra-Laguerre model that matches the process perfectly, and there is no error. There are six lines plotted corresponding to the Laguerre order used in the approximation. As can be seen in these figures, the improvement over linear models ranges from 65% to about 75% at 30 steps. Figure 23 illustrates that the standard deviation of the error reduction varies with the number of steps. There is a breaking point at 30- 40 steps. This is approximately the same number of steps required to obtain the Volterra linear models as discussed in the previous section. Thus, if a single set of step tests with 40 steps were 57 performed, error reductions of 75%  6% could be expected. This is significant in that 30-40 steps would be an acceptable number of steps for nonlinear models on industrial processes that provide noteworthy benefit over linear models. The Volterra-Laguerre identification also requires a nonlinear optimization. The question of convergability is present in this approach as well. Though there may be a solution, the problem has multiple local optima and may also be stiff or otherwise ill conditioned. There is no guarantee of convergence or repeated convergence to the same solution in the process of identification. Figure 24 and Figure 25 illustrate that even with these problems, significant improvement is repeatable and can be consistently obtained. All of the identification cases completed for this project demonstrated identification convergability. Thus, the conclusion drawn is that the global optima for the identification may or may not be the one achieved, but the improvement for whatever optima is obtained is significant over linear models. Figure 21 Volterra-Laguerre 2nd order error reduction -80.00 -60.00 -40.00 -20.00 0.00 20.00 40.00 60.00 80.00 100.00 0 50 100 150 200 AVE % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 58 Figure 22 Volterra-Laguerre 2nd order error reduction -80.00 -60.00 -40.00 -20.00 0.00 20.00 40.00 60.00 80.00 100.00 0 10 20 30 40 50 AVE % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 59 Figure 23 Standard deviations of 2nd order identification Figure 24 2nd order improvement for process with noise 0.00 20.00 40.00 60.00 80.00 100.00 120.00 140.00 0 50 100 150 200 STDEV % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 -20.00 -10.00 0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 0 20 40 60 80 100 120 140 160 180 200 AVE % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 60 Figure 25 2nd order standard deviation for identification of process with noise In the work done for this dissertation, the mean time between steps used was one half time to steady state (TSS). Thus, the total testing time based on this manipulated variable and this control variable and a 30 step test would be 15*TSS. Industry currently uses 10 to 20 steps to identify linear models depending on the process, process noise, linearity, and other factors. There is always a tradeoff between the number of steps used in the identification and the quality of the models. Figure 26 illustrates model predictions for a typical 20 step prediction set where the input-output data has not been used in the identification. This prediction is based on models identified from a different 20 step test. In this example, the first order Volterra-Laguerre model and impulse coefficient model (first order Volterra) perform comparably while the second order Volterra-Laguerre model does a much better job of fitting the actual process data. The number of steps required in the model is always an issue. The order of the desired model also affects the number of steps required. In step tests with fewer steps a low order 0.00 5.00 10.00 15.00 20.00 25.00 30.00 35.00 0 50 100 150 200 STDEV % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 61 Figure 26 2nd order predictions on 20 step data Volterra-Laguerre model performs as well as a higher order model. This is because with the limited data set, only identification of a smaller coefficient set is supported. In addition, as the number of Laguerre coefficients in the model increase, the model begins capturing noise effects. The extreme of this effect is the straight Volterra model. When there are many steps in the step test and as the Volterra-Laguerre order increases, the quality of the model improves. Figure 27 illustrates the objective function for a second order identification as a function of the Laguerre time scale factor and Laguerre order. Figure 27 also indicates that as the Laguerre order is increased, the relative magnitude of the improvement decreases. In this case, a Laguerre approximation with 4 coefficients would appear to be adequate. These results support an approach of identifying input-output data with continually increasing orders until the improvement caused by increasing the order becomes lower than some threshold. 0 2 4 6 8 10 12 14 -0.5 0 0.5 1 Time (TSS) Normalized NAMW Data Linear Volterra-Laguerre 2nd Order Volterra-Laguerre Impulse Coefficients 62 Figure 27 Effect of time scale factor 2nd order objective function As in the linear case, the selection of the Laguerre time scale factor is important. A plot depicting the effect of the time scale factor on the model error is shown in Figure 27. Since the identification requires a nonlinear optimization, the question of convergability is present here as well. Even though there may be a solution, the problem has multiple local optima and may again be stiff or otherwise ill conditioned. There is no guarantee of convergence or repeated convergence to the same solution in the process of identification. The work done in this dissertation is not meant to prove unconditional convergence criteria, but to illustrate practical convergability, convergence rates, and repeatability. Statistical information on the identification convergability is given in Figure 25. Figure 27 illustrates the shape of the objective function for the nonlinear fit of the time scale factor. A gradient-based optimization (nonlinear least squares) was performed to determine parameters for the model. It is apparent, by looking at the shape of 0 5 10 15 20 25 30 35 40 6 8 10 12 14 16 18 Time Scale a Squared Error R=2 R=3 R=4 R=5 R=6 R=7 63 the error function, that the solution of the optimization may not be the global optimum for each Laguerre order. However, the differences between the local optimum and the global optimum are not significant. CHAPTER V CASE STUDY II: VAN DER VUSSE REACTION The second case is the Van Der Vusse reaction: A k1 B k2C 2Ak3D The reaction is carried out in an isothermal CSTR as shown in Figure 28. This system has been studied by several researchers including Doyle III (1995), Maner et al. (1996), and Van Der Vusse (1964) and is considered a benchmark problem for nonlinear control. The system exhibits both nonminimum phase behavior, and nonlinear gain and dynamic effects. The problem addressed is to control the concentration of B by manipulating the inlet flow rate. The mass balance for the reacting components is given by:   V F C C k C k C dt dC Af A A A A a     2  1 3 (5.1) A B B B C V k C k C F dt dC    1 2 (5.2) B y  C (5.3) The process model parameters and steady-state operating conditions are given in Table 3 and Table 4 as given in Doyle et al. (1995). In Maner's work, Volterra models were obtained by 65 Figure 28 Hardware configuration of Van Der Vusse isothermal CSTR reactor Table 3 Van Der Vusse reaction parameters Parameters For Isothermal Van Der Vusse Reaction 1 k = 50 hr-1 2 k = 100 hr-1 3 k = 10 L/(mol hr) F = 1.12 L/hr V = 1 L Table 4 Van Der Vusse CSTR steady state operating conditions Steady-State Operating Conditions Van Der Vusse Reaction CAo = 3.0149 mol/L CBo = 1.12 mol/L Fo/V = 34.5939 L/hr y = 1.12 mol/L 66 using Carleman linearizations (Rough, 1981) of the first principle models. For this work, the same first principle models are used to generate I/O data that will then be used for identification. Open loop step response curves are shown in Figure 29 for a .5 normalized step up and step down of equal magnitude. As can be seen, the process is not symmetrical. The process has both changing gains and dynamic characteristics. Thus, the dominant root for the step up and the step down appear to be different. This is a common nonlinearity found in industry. In addition, this process exhibits nonminimum phase behavior with the inverse response. Inverse response is usually difficult to model and controller performance is affected dramatically if inverse response is not modeled properly. One of the advantages of using Volterra form nonlinear models is that changing gain, inverse response, and dynamic nonlinearities can be captured. The steady state gain loci for the Van Der Vusse reaction are shown in Figure 30. One of the characteristics that makes this system interesting is that the process has a maximum value with obtainable lower values on each side. The process gain is the tangent of this curve at any given operating point and thus changes sign. The implications of the gain sign change are that if the process is operating near the maximum, a linear controller will show a zero gain between the feed flow rate and the concentration of B. If the process is identified on the left side of the maxima a positive gain will be identified. Later, if the process is pushed to operate on the other side, the process gain will change signs and the controls will become unstable. This result is extremely important in industry. This behavior is common in many reacting systems, including combustion systems. In many of these systems, controlling close to the zero gain point is the optimum process operating point. It represents optimal stoichiometry for reactants in a reacting system. The ability to capture gains that not only change values but also change signs provides the ability to control many systems over a range that has proven to be uncontrollable with standard linear model predictive control. This result helps fill a hole in existing control technology for reacting systems. Currently, operators are 67 Figure 29 Step response curves for .5 normalized feed steps 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 Time (hrs) Normalized B Concentration 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 Time (hrs) Normalized B Concentration 68 Figure 30 Process loci for Van Der Vusse reactor forced to operate significantly away from the stoichiometric optimum to maintain stable control. The danger is that if they exceed the stoichiometric ratio with current controls, the process goes immediately unstable. This is a strong justification for using the higher order Volterra-Laguerre approach for reacting systems. Linear Volterra-Laguerre As in the previous case, the Volterra-Laguerre provided similar advantages for linear identification of the Van Der Vusse system. The main advantages were found in the reduced number of steps required to identify the linear models and the filtering characteristics of the Laguerre reductions. The Impulse and the Volterra- Laguerre systems were both able to capture the inverse response behavior of the system. Figure 31 illustrates a comparison of a typical step response for the process linear impulse coefficient model and Volterra-Laguerre model. As can be seen, neither of the linear models captures the nonlinear characteristics of the process. In Figure 31, the process was 69 Figure 31 Model comparison for steps on Van Der Vusse system 70 modeled with fifty coefficients in the impulse model and five coefficients in the Volterra- Laguerre model. Notice that the steady state gain of the second order Volterra-Laguerre model overshoots the process data in both directions. This is a result of the multiple step sizes used in the identification. The identification does the best fit possible for an array of step sizes. Thus, unlike the linear binary stepped identification, there is no guarantee that for a given step size there will be symmetrical errors around the step. Figure 32 illustrates the relative improvement for a noise free identification. As expected, the first order Volterra-Laguerre approaches zero improvement over the impulse coefficient identification as the number of steps increases. Figure 33 depicts the standard deviation of the linear identification based on the number of steps used in the step test. Figure 32 and Figure 33 show it is clear that two Laguerre coefficients are not adequate to describe the system. Even without process noise, it appears that between 10 and 20 steps are required to obtain an adequate linear model. When noise is added to the identification data, the situation changes. Figure 34 illustrates a typical 20 step identification data set with noise. Again, the noise is normally distributed random noise with a standard deviation equal to the average movement of the system without noise. Figure 35 illustrates the improvement of the Volterra-Laguerre over the impulse coefficient model. Because of the noisy signal, the Volterra- Laguerre identification demonstrates significant improvement at fewer steps and approaches that of the impulse model as the number of steps in the identification increases. Once again, this demonstrates the advantage of using the Volterra-Laguerre form even in linear identification. Figure 36 depicts the variation of identifications based on the number of steps in the step test. In this case, a 30 step test appears to be adequate. Again, it is important to note that there is a relationship between the number of required steps and the magnitude of the noise. For the Van Der Vusse process with no noise, 10 to 20 steps seemed to be adequate. In this case, three Laguerre coefficients appear to be adequate. The proposed identification approach would require successive identifications increasing the number of Laguerre coefficients until the relative 71 Figure 32 Improvement of linear Volterra-Laguerre Figure 33 Standard deviation of linear Volterra-Laguerre IDs -20.00 -15.00 -10.00 -5.00 0.00 5.00 0 10 20 30 40 50 60 Average % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 0 10 20 30 40 50 60 Standard Deviation of % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 72 Figure 34 Typical 20 step Van Der Vusse ID set with noise Figure 35 Improvement of linear Volterra-Laguerre on process with noise 0 2 4 6 8 10 12 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Time (TSS) Normalized Concentration of B Process w Noise Actual Process -10.00 -5.00 0.00 5.00 10.00 15.00 20.00 25.00 30.00 35.00 40.00 45.00 0 20 40 60 80 100 120 Average % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 73 Figure 36 Standard deviation of linear Volterra-Laguerre IDs with noise improvement diminished. In this process, there is significant improvement increasing from two to three Laguerre coefficients, but almost none in moving from three to four Laguerre coefficients. Using the prescribed approach, the identification would stop at four Laguerre coefficients. Figure 37 illustrates a typical plot of a linear prediction based on models obtained from a 30 step test using four Laguerre coefficients and fifty impulse coefficients. The models are compared to the process data without noise. Again, the impulse coefficient model clearly degenerates to a greater extent by the process noise than the Volterra-Laguerre model. In this case, no filtering was performed on the data. It is expected that the linear impulse model could be improved significantly by filtering the data. This step is eliminated in the Volterra-Laguerre identification. As in the previous case, the overall optimization of parameters is highly dependent on the Laguerre time scale factor. Figure 38 illustrates the effect of time scale selection on the objective function. The time scale factor was again included as one of the parameters in the identification. 0.00 5.00 10.00 15.00 20.00 25.00 30.00 0 20 40 60 80 100 120 Standard Deviation of % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 74 Figure 37 Linear identification comparison on Van Der Vusse data with noise Figure 38 Time scale factor "a" effect on linear optimization 0 5 10 15 20 25 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Time (TSS) Normalized Concentration of B Data Linear Volterra-Laguerre 2nd Order Volterra-Laguerre Impulse Coefficients 0 2 4 6 8 10 12 14 16 18 20 12 12.5 13 13.5 14 14.5 15 15.5 Time Scale "a" Squared Error R=2 R=3 R=4 R=5 R=6 R=7 75 The time scale factor is very important if you want to use a low order Laguerre approximation for your model. Second Order Volterra-Laguerre Once again, the main thrust of this research is to demonstrate the feasibility and practicality of dynamic nonlinear identification based on process step testing and I/O data sets. A significant area of contribution to the science is in dealing with the practical limitations and requirements of this identification. The Van Der Vusse system has been used extensively to illustrate nonlinear process characteristics. The plot in Figure 30 illustrated the process gain loci for the Van Der Vusse system. One interesting and important aspect of this system is that the process gain actually changes sign. At a reactant flow rate of 60, the process gain is approximately zero. This characteristic would make control of this process with a standard linear controller unstable across the range where the gain changes sign regardless of how mildly the controller was tuned. In many processes, the optimal place to operate the process is at this maximum value. This makes control locally uncontrollable. This process characteristic is common in many industrial applications where reactions are occurring, the most prevalent being combustion systems. If the manipulated variable were fuel rate and the control variable were process temperature, the plot of increasing fuel would show the same temperature characteristics. Many polymeric reactions also exhibit the same behavior. A second form of nonlinearity also exists in many processes. Not only does the gain change, but the dynamic response changes as well. The Volterra-Laguerre model does a nonlinear mapping of all of the dynamic coefficients and thus captures nonlinearities in the dynamics of the process. This nonlinear dynamic mapping provides advantage over a simpler gain scheduling model predictive control where model gains are changed but no dynamic changes are included. Figure 31 illustrates both the dynamic and the gain nonlinear effects for a simple step up and step down of the process and the models. 76 This process was identified using the same nonlinear identification technique as in the previous case. The approach uses the linear identification from the previous section and a second order correction matrix to yield second order nonlinearities. Figure 39 depicts the average improvement for ten identifications at each number of steps for the Van Der Vusse system without process noise. In this case, 20 to 30 steps appear to be adequate. Depending on the number of steps, a Laguerre order of three or four appears to be sufficient. These results demonstrate that 80% to 90% of the linear model error can be eliminated using the Volterra- Laguerre second order model. An indication of the variability of the identification is illustrated in Figure 40. Thirty to forty steps produces a fairly stable response. This is more steps than required in the linear identification but still reasonable. Even when fewer steps are done in the testing, the magnitude of the improvement outweighs the testing standard deviation in the data sets. It appears that regardless of the number of steps performed, the second order Volterra-Laguerre predicts the process significantly better than a linear model. Figure 41 illustrates a typical prediction plot where twenty feed rate steps are used in the identification. Figure 42 shows the steady state gain loci for the process, linear model, and second order Volterra-Laguerre. The feed flow rate and the concentration of B are scaled and normalized around the initial starting point for the step test. For this plot of gain loci, the mean scaled step size was 0.5. This plot demonstrates that the Volterra-Laguerre model performs well over the tested range and that the change in sign of the gain is captured. The Volterra-Laguerre model breaks down when moving out of the tested range. This verifies the well-known fact that extrapolation outside the testing range can be dangerous. Therefore, when an increased range of operation is expected, the process should be tested over the entire range. This differs from linear identification. In linear identification, when the process is nonlinear, it is necessary to test at the expected operating point. For processes with significant nonlinearity, process gain is very dependent on the operating point. In addition, the size of the 77 Figure 39 Volterra-Laguerre 2nd order error reduction Figure 40 Volterra-Laguerre standard deviation of 2nd order ID 30.00 40.00 50.00 60.00 70.00 80.00 90.00 100.00 0 10 20 30 40 50 60 Aveerage % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 0.00 5.00 10.00 15.00 20.00 25.00 30.00 0 10 20 30 40 50 60 Standard Deviation of % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 78 Figure 41 Model prediction of noise free Van Der Vusse system Figure 42 Comparison of steady state gain loci for process and discrete models 0 2 4 6 8 10 12 -1 -0.5 0 0.5 Time (TSS) Normalized Concentration of B Data Linear Volterra-Laguerre 2nd Order Volterra-Laguerre Impulse Coefficients -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Feed Flowrate (F) Concentration of B Process Linear Model Volterra-Laguerre 2nd Order Mean Step Size Test Range 79 steps in the step tests on these processes can significantly affect the gain of the model. One advantage of the second order Volterra-Laguerre model is an increased range of acceptable predictions over the linear model. A second advantage is the simplicity of stepping over the expected operating range of the process. The size of the steps in the step test should be limited by the constraints of process operation, not based on fitting a linear model to a nonlinear process. This actually makes defining the step tests much easier in the sense that the process nonlinearity does not require minimizing the step size around the current operating point. Making larger steps helps with both the linear and is required for the nonlinear portions of the model. Any errors in the linear model caused by the larger step sizes are picked up in the second order portion of the model. Figure 43 illustrates the results of the second order Volterra-Laguerre identification for the Van Der Vusse process with noise added. The process data used are the same as the noisy data used in the linear identification. Several aspects of these results are interesting. First, as with the previous cases, there is significant improvement in the second order Volterra-Laguerre model's ability to predict process behavior over the impulse model. Second, as before, a Laguerre order of two is inadequate but still better than the linear model. Third, a continual increase of Laguerre order does not necessarily improve model quality. A Laguerre order of three does better than an order of six or seven. This is true because as the number of coefficients in the model increases, the identification is able to capture unwanted noise effects. Therefore in noisy systems, increasing Laguerre order R can actually produce worse models than lower Laguerre orders. Again, the proposed procedure will alleviate this problem. In this case, the Laguerre order in the identification would be increased until no significant improvement was shown in the model, thus resulting in a Laguerre order of four. The standard deviation of the identification based on the number of steps in the step test for the noisy data is plotted in Figure 44. The standard deviation suggests that 40 steps would be optimal but, even fewer steps would provide significant improvement. A typical prediction plot of the noisy system is shown in Figure 45. 80 Figure 43 Volterra-Laguerre 2nd order error reduction- process with noise Figure 44 Volterra-Laguerre standard deviation of 2nd order ID - process with noise 50.00 60.00 70.00 80.00 90.00 0 20 40 60 80 100 120 Average % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 0.00 5.00 10.00 15.00 20.00 25.00 0 20 40 60 80 100 120 Standard Deviation of % Error Reduction # of Steps R=2 R=3 R=4 R=5 R=6 R=7 81 Figure 45 Typical 30 step model prediction of Van Der Vusse system with noise The plot displays the linear impulse and Volterra-Laguerre models and nonlinear Volterra- Laguerre prediction, as well as the actual process data. In this plot, the process data are not the same data used in the identification and thus demonstrates the ability of the models to predict data sets not used in the identification. Finally, Figure 46 illustrates the effect of the time scale factor on the identification objective function. Once again, it is clear that the time scale factor is a critical parameter in the identification process. If lower Laguerre orders are expected to be used, using the optimal time scale factor becomes critical. 0 5 10 15 20 25 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Time (TSS) Normalized Concentration of B Data Linear Volterra-Laguerre 2nd Order Volterra-Laguerre Impulse Coefficients 82 Figure 46 Effect of time scale factor "a" on 2nd order optimization 0 10 20 30 40 50 60 70 80 90 100 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 Time Scale a Squared Error R=2 R=3 R=4 R=5 R=6 CHAPTER VI CASE STUDY III: IDENTIFICATION OF EXPERIMENTAL PH DATA The final identification was performed on experimentally obtained PH data from a PH control experimental apparatus. This chapter will illustrate the procedure for identifying the nonlinear Volterra-Laguerre models from input-output data. Thus, the goal of this chapter is to show both the methodology of the identification and the results of an experimentally tested system. The process consists of one feed tank containing a base solution, dissolved concentrated sodium hydroxide in water, and another feed tank with acidic solution containing hydrochloric acid in water. A hand valve is used to set the acid flow to a general value. Base flow is adjusted as the manipulated variable of the control system. A schematic of the experimental apparatus is shown in Figure 47. Figure 47 PH control experimental apparatus Acid Solution Base Solution Waste Solution PH Meter 84 The PH control problem exists in a broad range of industries. It has proven to be a nagging problem for control engineers. This system has several characteristics that make it interesting. First, PH control is typically considered off-limits for linear controls because of the extreme process nonlinearity. Second, the logarithmic system, process drift, and other properties make PH control notoriously difficult. Unlike the previous two cases in this work, this data was obtained from the process as input-output data through a distributed control system (DCS). No attempt to do fundamental physics-based modeling was done. Because a fundamental model does not exist, this experimental system illustrates typical problems encountered during industrial process operation. These problems include the following: The noise in the system is not random noise and is not well defined. There is drift in the base feed solution caused by adsorption of CO2 from the surrounding air. In addition, flow rates of the base and acid solutions vary with the level in the tanks. Thus, the system is slightly time variant. There is drift in the PH sensor (a combination of three sensors) also causing variation in time. The valve used as the manipulated variable has nonlinear characteristics and exhibits a small amount of hysteresis "stickiness," depending on step direction and magnitude. The feed solutions are prepared off-line in batches and have slightly different characteristics. The feed system consists of two tanks and centrifugal pumps, with a control valve connecting the two, and a hand control valve on each line. There was strong coupling between the acid and base flows caused by line pressure drop effects causing additional nonlinearities in the system. These problems are typical of what would be found in an industrial setting on an active production process. In this case, there was no fundamental dynamic model. The process consisted of multiple pieces of equipment, each with their own set of dynamics, nonlinearities, and peculiarities. The total process was a combination of mechanical, chemical, environmental, and electronic effects. 85 The objective of this case is not to show that the Volterra-Laguerre model is "ideal" for PH systems but to demonstrate that in systems that are not well defined, the Volterra-Laguerre nonlinear identification can provide significant improvement over linear systems. Second, the Volterra-Laguerre system can be identified from easily obtainable input-output data sets in these systems. In documenting this case, the "identification algorithm" will also be defined and illustrated. The identification algorithm represents a method of how to perform a nonlinear identification using Volterra-Laguerre systems. It includes all of the steps required from step testing to final model validation. The Volterra-Laguerre characteristics have been illustrated extensively in the previous sections of this dissertation. This chapter will present the procedure and results of what would be a typical industrial identification. Procedure The first component in any process control identification and modeling project should be to study the process and its characteristics. Some key elements include the following:  Identify as many peculiarities in the process as possible that may affect the process response.  Identify any known process disturbances, measured or unmeasured.  Check for instrumentation problems and characteristics (any problems in instrumentation will result in a degradation of the final models).  Determine the expected response based on the physics of the process.  Learn the objectives of the process operation as well as process constraints.  Determine with the help of the process engineers and operators the goals of the control system. 86 In the PH control experiment, the process analysis revealed several characteristics that affect process dynamics, including the following:  The PH control is inherently nonlinear because PH is a logarithmic scale. Thus, linear incremental changes in base flow rate will produce logarithmic, not linear, responses in the PH.  The process noise in the system is not pure white noise and is not well defined.  There is drift in the base feed solution caused by adsorption of CO2 from the surrounding air.  The feed solutions are prepared off-line in batches and have slightly different characteristics.  Flow rates of the base and acid solutions vary with the level in the tanks. Thus, the system is slightly time variant.  There is drift in the PH sensor (a combination of three sensors) also causing variation in time. The valve used as the manipulated variable has nonlinear characteristics and exhibits a small amount of hysteresis "stickiness," depending on step direction and magnitude.  The feed system uses two tanks and centrifugal pumps with a control valve connecting the two, and a hand control valve on each line. There is strong coupling between the acid and base flows caused by line pressure drop effects that cause additional nonlinearities in the system. The second component in the identification process is to evaluate lower level control loops in the process. Aspects of this include the following:  Ensure control loops are configured properly and are consistent with the overall objectives of the control strategy. Every control loop should have a function. The two primary functions for regulatory controls are to Linearize a portion of the process and to 87 reject process disturbances that show up in measured variables quickly. An example of this is steam flow control. The flow controller will remove any nonlinearity due to valve trim type and reject high frequency disturbances caused by fluctuating steam pressure.  Regulatory control loops should be tuned so the process is as responsive as possible without being cyclic or approaching instability. Regulatory loop tuning will have an effect on the dynamics in the resulting models. Tuning should be done before identification of the process. It is important to recognize that any changes in process equipment, regulatory control loop configuration or tuning, or basic operating parameters (feedstock, catalyst, etc.) can significantly change the process dynamics and process gains. Thus, any changes should be made before step testing of the process. Any changes made after the identification will show up as model mismatch between the identified process model and the actual process response. In this PH control case, there were no regulatory controls. The dynamics of this process are relatively fast (~30 seconds) so there were no high frequency disturbances to be rejected. Thus, in this case, regulatory configuration and tuning were not an issue. Determination of Process Test Conditions The third step in the identification process is to define the operating conditions where testing is to be done and to establish the range of step sizes. Since the Volterra-Laguerre model is a regression-based model, it is necessary to define the mean operating point for the test as well. Ideally, this point will be close to what will be the desired process setpoint. When the process is operating at exactly the mean testing point, the linear portion of the Volterra-Laguerre model represents the process without the need for nonlinear corrections. In reality, the actual location of operation after the controls are implemented may change. Improved control usually makes it possible to operate the process near process constraints. In most processes, the optimal operating point is at an operational constraint. This operational shift can cause significant model error in 88 the linear portion of the model if the process is significantly nonlinear. It is through the nonlinear correction that the process model accounts for nonlinearity. As discussed in the previous chapter, it is important to include as much of the expected operational area as possible in the nonlinear identification test. As a result, the step size magnitude should be as broad as possible. The mean step size and the maximum and minimum limits on the manipulated variable should be selected carefully. Extrapolation of the model predictions outside the tested area often lead to significant model mismatch and thus poor process control. The PH system testing range was limited by the PH sensor accuracy. The sensor was inaccurate at PH values above 10 and below 4. The mean valve position of 11% was decided upon after preliminary testing. The mean step size of 1.5% was selected because it resulted in a PH data variation of about 5 PH units or 5 orders of magnitude in concentration. This range is broad enough to illustrate significant nonlinearities. The mean time between steps was set at 30 seconds based on the dominant time constant in the system response. This was determined via preliminary step tests as well. Process Step Testing The next component in the identification process is the actual step testing and data collection. As discussed in Chapter 3, the step frequency and magnitude is based on pseudo-random multilevel sequences (PRMS). There are several items that should be checked before the step testing begins, including the following:  The process should be at the predetermined mean operating point for the test.  The regulatory controls should be in the configuration and state they will be in when the controls will be turned on. For example, if a regulatory loop will be in manual mode when the model predictive control is implemented, it should be in manual mode during the step test. There are many cases where the fallback position of the control system 89 (when the MPC is off) will have additional regulatory control loops activated to help maintain control of the process. These control loops should be off during step testing.  Unmeasured disturbances should be minimized.  Data should be collected for all variables that might become a manipulated, disturbance, or control variable in the model. In many cases, there are several variables (or combination of variables) that could become important for the controls after the step testing. Often, the plant operators have never moved the process enough to obtain dynamic cause and effect data. Additional control points for the process often become apparent during the test or during the ensuing analysis.  The data collection frequency must be high enough to capture the highest frequency dynamics of the process response to be modeled.  The distributed control system must have a PRMS test sequence testing program implemented that allows for the automatic execution of the manipulated variable moves. This is critical. Manual stepping of the process does not usually work well. It is difficult for even the best operator to implement random steps. The operator has been trained to make multiple moves on a process to deal with process coupling. Dealing with coupling between process variables is one of the main features of the MPC. The MPC must have models that represent interactions between all independent variables and dependent variables in order to address coupling in the process. When steps are correlated based on process response, it becomes impossible to identify the interactions and the resulting models can be seriously compromised. Once the above items have been completed, the step test can begin. During the step test, the process should be continuously monitored with special attention given to events that will compromise the data being collected. These events include the following: 90  Manipulated and control variable saturation which cause errors in both process gain and process dynamics of the identified model.  Any unmeasured disturbances that might skew the data.  Any operator intervention on unmonitored or hand controlled variables. Examples of this include hand valves, pump bypass valves, pump motors, etc. In short, any event that could affect the process data that is not recorded in the data set itself should be recorded in a testing log so that problem sections of the resulting data set can be removed before identification. The minimum step test duration is based on the number of steps required to obtain an acceptable model, the mean frequency of the steps and the number of independent variables being tested. An optimal test would include enough moves so that data not used in the identification could be used to evaluate the identified models prediction capabilities. In practice, it is often necessary to do preliminary identification of the process during the step testing to determine if enough data have been collected. Guidelines for the number of steps required for the nonlinear identification was discussed in previous chapters. In the PH case, data were collected at a sample rate of once per second. The dominant time constant for the process appeared to be approximately 30 seconds. The one-second sample rate allowed enough resolution for higher frequency disturbances to be observed. A graph showing raw data from a step test of this system is shown in Figure 48. Data Conditioning Once the process step testing is complete and the data is collected, the data needs to be preconditioned before the identification is performed. The elements of preconditioning include the following: 91 Figure 48 PH step test raw data  Removing data slices for periods where unmeasured disturbances or other process events affecting data quality occurred and detrending the data. These removed periods should be recorded in the testing log as described in the previous section and as was done in the pH case.  One of the characteristics of MPC is that it assumes the process being controlled is time invariant. Thus, if the manipulated variable is set at the same position at two different times, the control variable will return to the same value. In practice, this is often not exactly true. An example of this time variance can be seen by examining the raw process data in Figure 48. It is clear that there are sections of the process data where the PH is drifting significantly while the valve position is being varied around the mean testing point. The drift can be caused by any one of many unmeasured effects. Examples of these effects include instrumentation drift and environmental changes. One of the identified causes, in the PH case, is CO2 adsorption from the air into the base solution. This causes changes in the feed base solution PH, resulting in product 5 7 9 11 13 15 17 19 21 23 25 12:00 12:30 13:00 13:30 14:00 14:30 15:00 15:30 16:00 16:30 Time Valve Position % Open 0 1 2 3 4 5 6 7 8 9 10 PH VP PH 92 PH drift. Another source of drift noted in the PH process was due to changing levels in the feed tanks. Centrifugal pumps are used to force feed solutions through the system. The feed system is essentially a batch process where feed solutions are prepared and loaded in batches. The levels of the feed tanks and the feed rates of the acid and base drift with time. In a SISO system, drift is usually easy to identify and thus easily removed. As the number of independent variables being stepped increases, it becomes more difficult to identify the drift. One approach used to solve this problem stems from the separation of the dynamics. Drift is usually a lower frequency process than the process dynamics of the system to be modeled. One effective technique used to remove this low frequency drift is to process the data through a high pass filter. As the frequency of the dynamics of the drift approach the frequency of dynamics of the process, this technique fails. In some situations, drift is a