- Email: [email protected]
Benefits of factorized RBF-based NMPC
Computers and Chemical Engineering 26 (2002) 1185– 1199 www.elsevier.com/locate/compchemeng
Benefits of factorized RBF-based NMPC Sharad Bhartiya, James R. Whiteley * School of Chemical Engineering, Oklahoma State Uni6ersity, 423 Engineering North, Stillwater, OK 74078 -5021, USA Received 13 June 2001; received in revised form 6 September 2001; accepted 25 February 2002
Abstract The computation and application benefits of factorized RBF-based nonlinear model predictive control (NMPC) are presented. The NMPC algorithm derives its computational efficiency by factorizing the radial basis function (RBF) model response. A brief description of the factorized RBF-based NMPC algorithm is provided. Theoretical computation benefits are quantified for both SISO and MIMO formulations. Computation and application benefits of the algorithm are documented for a 4 ×4 MIMO subset of the Eastman Challenge problem. Results confirm computation requirements are reduced by more than an order of magnitude relative to application of traditional MPC using a non-factorized RBF model. The expected control performance benefits from using a nonlinear process model are also achieved. © 2002 Elsevier Science Ltd. All rights reserved. Keywords: Nonlinear MPC; RBF model; Multivariable control
1. Introduction The current generation of commercially available model predictive control (MPC) packages relies on linear process models to provide future predictions. Excellent performance is realized as long as the plant operates close to the conditions used to create the linear approximation of the process. The goal for the next generation of control software is to provide similar capability across the whole range of possible plant operating conditions. The notion of using a nonlinear model within the MPC paradigm has led to an active interest in the development of nonlinear model predictive control (NMPC) methods. A common approach to the nonlinear optimization required in NMPC is based on successive linearization of nonlinear models. Garcia (1984) proposed a nonlinear QDMC algorithm, a simple extension of DMC/ QDMC based on online successive linearization of a mechanistic nonlinear model. Nonlinear MPC using closed-loop state estimation by an extended Kalman filter has been proposed by Lee and Ricker (1994). Gattu and Zafiriou (1995) augmented the system states with stochastic states to account for modeling errors * Corresponding author. Fax: + 1-405-744-6338. E-mail address: [email protected] (J.R. Whiteley).
and disturbances. Banerjee, Arkun, Ogunnaike, and Pearson (1997) describe a method of state estimation for nonlinear systems that are subject to multiple operation regimes and make transitions between them. The nonlinear process is approximated by a linear parameter varying system which consists of local linear models. Krishnan and Kosanovich (1998) also present a multiple model based MPC scheme. The linear time invariant models are computed offline along a pre-defined reference trajectory of a batch process. The NMPC method addressed by this paper utilizes an artificial neural network (ANN) for the process model. The nonlinear modeling capability of ANNs is well documented (Cybenko, 1987; Hartman, Keeler, & Kowalski, 1990). Multi-layer perceptron and radial basis function (RBF) networks are the most extensively utilized ANNs in control and identification applications. Hussain (1999) provides a summary of a number of control applications reported in the literature. The current work employs an RBF network to model the process. Development of process models using RBF networks has been extensively explored by Seborg and co-workers. Pottman and Seborg (1992) applied an RBF based predictive control strategy to an experimental pH neutralization process. Bomberger, Seborg, and Ogunnaike (1998) studied identification of a multivariable copolymerization model with emphasis on design
0098-1354/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 2 ) 0 0 0 2 9 - 7
1186
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
of input sequence. Gurumoorthy and Kosanovich (1998) incorporated process knowledge via regularization theory to improve the prediction capability of RBF networks. Also, online adaptation of RBF networks has been studied. Chen, Billings, and Grant (1992) used a hybrid clustering algorithm along with recursive least squares to update the network centers and weights, respectively. Sanner and Slotine (1992) applied model update rules for RBF networks in an adaptive control framework. Because of the computational cost associated with nonlinear optimization, there are practical incentives to integrate a nonlinear model in an MPC scheme in the most efficient manner possible. RBF networks provide good modeling capability but are not ‘‘compact’’ from a traditional serial computing standpoint. RBF models for nontrivial processes require a large number of parameters and involve many more mathematical operations than time series and other common types of models. The factorized RBF-based NMPC method (Bhartiya & Whiteley, 2001) was developed to reduce the computational penalty associated with repeated RBF model predictions during nonlinear optimization searches. The remainder of this paper is organized as follows. The essential points for the factorized RBF-based NMPC algorithm are presented in Section 2. The computation benefits for both SISO and MIMO formulations of the algorithm are discussed in Section 3. For illustrative purposes, the algorithm is applied to a 4× 4 MIMO subset of the Eastman Challenge problem (Downs & Vogel, 1993). An overview of the Eastman process and development of the NMPC controller are described in Section 4. Benefits associated with the NMPC controller are presented in Section 5 followed by conclusions in Section 6.
2. Factorized RBF-based nonlinear MPC algorithm Feedforward RBF networks approximate process dynamics by using past process information as network inputs. Thus, yˆk = F(yk − 1, …, yk − Ny, uk − 1, …, uk − Nu )
(1)
where a time-delay of unity is assumed between the model output, yˆk, and the previous process inputs, uk − 1. Ny and Nu are integers defined by the order of the model. Scalars uk − 1, …, uk − Nu represent the sequence of inputs used by the model. Function F depends on the network architecture and the type of activation function employed by the nodes. Model orders, Ny and Nu, may be determined by methods discussed by He and Asada (1993) (based on Lipschitz numbers), and Kennel, Brown, and Abarbanel (1992) (using false nearest neighbors). Bomberger and Seborg (1998) demonstrate
use of these methods in conjunction with an RBF network. Bhartiya and Whiteley (2001) proposed a factorized RBF based NMPC algorithm that exploits the factorizability property of the Gaussian function yielding a numerically efficient algorithm. Model identification is achieved via training of an RBF network using Gaussian nodes to provide sequential predictions over a prediction horizon of length p. The NMPC controller is directly parameterized in terms of the RBF network parameters. This allows analytical expressions for the gradient and Hessian of the objective function.
2.1. Model formulation The RBF network consists of an input layer, a hidden layer, and an output layer. In the input layer, unweighted inputs, x, are directly transmitted to the hidden layer nodes. Each hidden layer node consists of a RBF. In the current work, the Gaussian function is employed,
gj (x) =exp −
x− tj 2 |2
(2)
where tj is the center of the jth hidden node. It is assumed that each node is of fixed width |. The output layer performs a weighted summation of hidden node outputs to give the network output, m
yˆi (x) = % wij gj (x)
(3)
j=1
where m represents the number of hidden nodes and wij is the weight between the ith output node and the jth hidden node. The node centers, tj, are fixed in an unsupervised training step by partitioning the input data space using a clustering algorithm (e.g. k-means). Node width, |, accounts for the desired degree of overlap between the clusters and is fixed using heuristics (e.g. mean distance between cluster centers) (Hush & Horne, 1993; Moody & Darken, 1989). Finally, the network weights are obtained by linear regression to approximate the input–output mapping provided by the data. The number of nodes, m, directly contributes to the model complexity (number of free parameters). Henrique, Lima, and Seborg (2000) discuss node pruning techniques based on orthogonal least squares which yield parsimonious models. Often, the training data consists of noisy measurements. This may cause the minimization of prediction error to be an ill-posed problem. Regularization theory can be used to address this problem by incorporating smoothness constraints or a priori knowledge (Poggio & Girosi, 1990). The RBF model in the proposed NMPC scheme provides non-iterative sequential predictions over a prediction horizon of length p by avoiding dependency of future model predictions on previous model predic-
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
tions. Future predictions, yˆ k, are related to past measurements yi and inputs ui, delayed by p samples and input moves, Dui, as follows, yˆ k/k − p = F(yk − p, …, yk − p + 1 − Ny, uk − p − 1, …, uk − p + 1 − Nu, Duk − p, …, Duk − 1)
(4)
where the control move at the kth instant, Duk = uk − uk − 1
(5)
The delayed input/output measurements, viz. yk − p, …, yk − p + 1 − Ny and uk − p − 1, …, uk − p − Nu, respectively, provide a reference to the state of the system p samples in the past. Based on the model in Eq. (3), the following form of the input vector is chosen for a single-input, single-output system, xk = [yk − p …yk − p + 1 − Ny uk − p − 1…uk − p + 1 − Nu Duk − p…Duk − 1]
(6)
Fig. 1 illustrates the input vector with elements located on a timeline. Measurable disturbances can be accounted by augmenting the input vector x to include the disturbance variables. The RBF prediction of the process output at instant k is then obtained by Eq. (3). The form is general and can accommodate both multiple inputs and outputs.
1187
2.2. NMPC controller formulation The algorithm computes the manipulated variable profile over a control horizon by optimizing an objective function defined over the prediction horizon, subject to constraints. Only the first move is implemented and the procedure is repeated at every sampling instant. The following objective function is used: p
c−1
i=1
i=0
= % Yi (rk + i − yˆ k + i )2 + % \i (Duk + i )2
(7)
Variables p and c represent the prediction and control horizons, respectively. Yi and \i denote the error penalty and move suppression factors at the ith instant. The MPC control law can be stated as, arg(min ) Du , Du ,…,Du k k+1 k + c − 1 such that ymin 5 yˆ k + i 5 ymax (8) Dumin 5 Duk + i 5 Dumax umin 5 uk + i 5 umax Future model predictions, yˆ k + i, depend on past control moves and the future control move variables, Duk + i (see Eq. (4)). The future control moves, Duk, …, Duk + c − 1, represent the decision variables for the optimization problem in Eq. (8). For calculation purposes, it is desirable to express yˆ k + i such that the unknown decision variables appear explicitly in the objective function. The key idea in the factorized RBF
Fig. 1. Timelines showing inputs to the p-step control model. For this example, p = 4, Ny =3, and Nu =2. Predicted outputs are generated from known information only, previous model predictions for y are not used in model input.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1188
model lies in expressing the model prediction, yˆ k + i, as an inner product of two vectors. The unknown decision variables of the nonlinear program (Eq. (7)) are contained in one vector and all other known past quantities, including the network weights in the other. Thus, the RBF output can be rearranged as follows: m1
yˆ k + i = % wj exp(past +future) j=1
Æ w1 exp(past) ÇT Æexp(future) Ç Ã Ã Ã Ã =Ã Ã Ã Ã Ã Ã Ã Ã Èwm 1 exp(past) É Èexp(future) É
(9)
yˆ k + i =yˆ p,k +Tiyˆ f,k + i
(10)
Subscripts p and f refer to the fact that the corresponding factors contain all known (past) and unknown (future) terms, respectively. Thus, only yˆ f,k + i needs to be computed during every function call by the optimization algorithm. For a SISO problem, the factors for the predictions, yˆ k + 1, …, yˆ k + p, based on c future moves, Duk, …, Duk + c − 1, take the following form:
yˆ p,k + i =
w exp −
1 |2
Ny − 1
% (yk − p − j + i 1
j=0
Nu − 1
−tyj + 1)2 + % (uk − p − j + i 1 −tuj )2
n
j=1 p−1
+ % (Duk − p + j + i 1 − tj + 1)2
(11)
j=0
The symbol 1 represents a column vector of size m with unity elements. The operator ‘‘ ’’ is used to denote the Hadamard product (element by element multiplication). Note that when the index, i, equals p, the final summation drops out since j varies from 0 to − 1. The factor yˆ f,k + i that contains the future moves the optimization algorithm must calculate is:
yˆ f,k + i = exp −
i−1
i= 1, …, c yˆ f,k + i = exp −
n
1 2 % (Du 1 − tDu , p − j) | 2 j=0 k+i−j−1 i−1
n
(12a)
1 2 % (Duk + i − j − 11 − tDu , p − j) | 2 j=i−c
i=c+ 1, …, p
the gradient and Hessian of the objective function. The availability of analytic expressions significantly reduces the computational requirements associated with the algorithm. In addition, the separation of the decision variables in the model prediction ensures that only the unknown parts of the objective function, and the gradient and Hessian required by the sequential quadratic programming (SQP) algorithm are recalculated during optimization.
3. Computational efficiency To quantify the computational efficiency of the factorized approach, we first consider the number of floating point operations (flops) required when the SQP solver calls the function to produce the future p model predictions yˆ k + i required by the controller objective function (Eq. (7)). The total number of flops depends on, (a) the formulation of the RBF model, factorized or non-factorized, (b) the number of nodes in the RBF network, m, (c) the size of the control problem (the number of control outputs, No, and control inputs, Ni ), (d) the orders of the model, Ny and Nu, and (e) the length of the prediction and control horizons, p and c, respectively. Counting addition, subtraction, multiplication, division, and exponentiation as single floating point operations, the total number of flops required to evaluate vectors yˆ p,k + i and yˆ f,k + i, can be determined directly from Eqs. (11) and (12), respectively. For the SISO case, it is straightforward to show that the number of flops, Cp,SISO, required to generate yˆ p,k + i is, Cp,SISO = pm(3Ny + 3Nu + 3p −2)
Likewise, the number of flops, Cf,SISO, required to generate yˆ f,k + i is, Cf,SISO = mc(c+ 2)+ m(p −c)(2c + 3)
The exponential function in Eqs. (11) and (12) implies a term by term application to each element of the column vector in the square brackets. Column vector tj is constructed by using the jth element of all m centers. Thus, any future prediction, yˆ k + i, can be computed using Eqs. (10)–(12). Extension to the MIMO case is straightforward (Bhartiya & Whiteley, 2001). The factorized form provides analytic expressions for
(14)
Analogously, for the Ni × No MIMO problem, assuming each input and output have the same associated orders, Ny and Nu, respectively, the following results are obtained, Cp,MIMO = pm(3+3NoNy + 3Ni Nu + 3Ni p− 5Ni )
(12b)
(13)
(15)
and Cf,MIMO = mc[Ni (c− 1)+ 3]+ m(p−c)(2Ni c+ 3) (16) For every function call when using the factorized model, calculation of the p future model predictions using Eq. (10) requires, Cfactorized = Cp + Cf + p(2m −1)
(17)
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1189
However, the model prediction function is called numerous times at every control step. Thus, if n total function calls are made by the SQP solver for computation of the control inputs, the factorized approach will need, Ctotal,factorized = Cp + nCf + np(2m − 1)
(19)
total flops (since Cp is only calculated once), whereas the non-factorized model predictions will require, Ctotal,non-factorized = n(Cp + Cf − 3mp) Fig. 2. Computational speedup associated with the factorized RBF model increases with number of model predictions required at each time step by the NMPC optimization algorithm. Plot prepared with m = 230, p = 20, c= 5, Ni = 4, No = 4, Nu = 2 and Ny = 1.
flops. The reduction in computational expense associated with the factorized RBF model is difficult to appreciate from inspection of Eqs. (13)–(20). A sensitivity analysis using typical variable values was performed to characterize the benefits of model factorization in terms of a speedup factor defined as, Speedup factor=
Fig. 3. Contour plot of speedup factor as a function of the prediction and control horizons. Plotted for m = 230, n = 100, Ni = 4, No =4, Nu = 2 and Ny =1.
Fig. 4. Contour plot of speedup factor as a function of the number of MIMO inputs and outputs. Plotted for m= 230, n =100, Ni =4, No = 4, Nu =2 and Ny = 1.
floating point operations. The last term in Eq. (17) is associated with the operations performed to calculate the inner product in Eq. (10) (m products followed by (m− 1) additions). For the non-factorized model, the number of flops needed for calculation of p model predictions is,
Cnon-factorized = (Cp +Cf −3mp)
(18)
(20)
Ctotal,non-factorized Ctotal,factorized
(21)
A speedup factor of 10 indicates an order of magnitude fewer computations are required to produce model predictions using a factorized RBF model. The sensitivity analysis reveals that the speedup factor is independent of the number of nodes, m. The speedup factor is dependent only on n, p, c, Ni, No, Nu, and Ny and not on the size of the RBF model. The speedup factor increases in a first-order manner with the total number of model predictions or function calls, n, required to solve the MPC control law, Eq. (8). Fig. 2 was generated for a 4×4 MIMO system with prediction and control horizons of 20 and 5, respectively, input and output model orders of 2 and 1, respectively, and an RBF model with 230 hidden nodes. The trends and speedup values shown in Fig. 2 (as well as Figs. 3 and 4 to follow) are representative for a wide range of variable values typical for industrial MPC applications. As indicated in Fig. 2, a five-fold improvement is realized after only 13 function calls. As the number of function calls increases, the speedup factor approaches an asymptotic limit of approximately seven. The number of function calls depends on the efficiency of the nonlinear optimization search. A contour plot of the speedup factor as a function of the lengths of the prediction and control horizons, p and c, respectively, is presented in Fig. 3. The variable values used to generate the plot are the same as those for Fig. 2 with n fixed at 100. Typical practice is to implement MPC with p/c ratios ranging from 3/1 to 5/1. As revealed in Fig. 3, the speedup factor is nearly constant for a fixed p/c ratio and varies from 6 to 7 for typical MPC applications. The size of the MIMO control problem impacts the speedup factor through Ni (number of inputs or MVs)
1190
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
and No (number of outputs or CVs). As seen in Fig. 4, all combinations with between 3– 15 inputs and 3–20 outputs produce speedup factors between 6.8 and 7.8. Fig. 4 was generated in the same manner as Figs. 2 and 3 with n=100 and a p/c ratio of 4 with p =20. As one would expect, the speedup factor increases with increasing complexity or size of the MIMO problem. A review of all results generated by the different sensitivity analyses indicates a speedup factor of approximately seven can be expected for a typical MIMO control problem. However, there is an additional computational benefit associated with the factorized RBF model. Analytical expressions can be provided for the gradient and Hessian of the objective function. Significant computational savings are realized by elimination of the need for numerical approximations. With cNi decision variables in the optimization program in Eq. (8), a single numerical approximation of the gradient requires cNi additional function calls. Likewise, a single numerical approximation of the Hessian requires (cNi )2 calls. The actual computational savings are magnified based on the need for repeated evaluations of the gradient and Hessian to satisfy convergence criteria for both the numerical approximations (inner computation loop) and overall optimization search (outer computation loop).
4. Application of RBF-based NMPC to the Eastman problem The Eastman challenge problem (Downs & Vogel, 1993) was employed to characterize the computational and application benefits of the factorized RBF-based NMPC algorithm. The process (Fig. 5) involves producing two liquid products, G and H, from four gaseous reactants, A, C, D, and E.
4.1. Plantwide control strategy A variety of decentralized and MPC control strategies have been proposed for the Eastman process. Ricker (1996) discusses the key issues involved with design of the plantwide control system with references to many of the previously published designs. The current work is not intended to provide another comprehensive solution to the Eastman problem. Rather, the control system described below was developed strictly to examine performance characteristics of the proposed NMPC algorithm on a meaningful scale. The SISO strategy previously reported by McAvoy and Ye (1994) was used for the base regulatory control system. The SISO loop pairings recommended by McAvoy and Ye are presented in Table 1. Two minor modifications were made to McAvoy and Ye’s scheme:
Fig. 5. Schematic of the Eastman challenge process.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199 Table 1 Input–output pairing determined by McAvoy and Ye (1994) Controlled variable
Manipulated variable
Reactor temperature Reactor pressure Reactor level Separator level Stripper temperature Stripper level Production rate G/H composition in product E composition in product Inert B composition in purge stream Recycle flow
Reactor coolant temperature setpoint A feed rate setpoint E feed rate setpoint Underflow rate setpoint Steam flow setpoint Sripper underflow rate setpoint A&C feed rate setpoint D/E ratio setpoint
Compressor power
Stripper temperature setpoint Purge flow rate setpoint Condenser cooling water temperature setpoint Recycle valve
Table 2 Controlled and manipulated variables for the RBF based MPC controller MPC controlled variables (CVs)
MPC manipulated variables (MVs)
Reactor pressure Composition of B in purge Production rate A/C mole ratio in reactor feed
Reactor temperature setpoint Purge rate setpoint A&C feed rate setpoint A feed rate setpoint
1. The setpoint of the condenser cooling water return is used to control the condenser temperature (Gain=5.0%/°C; integral time=50 min). 2. The flow controller for the A&C feed stream (stream 4) was tuned more aggressively (Gain= 0.3%/kscmh; integral time=0.09 min). Controlled variables were selected to meet the following servo specifications suggested by Downs and Vogel: 1. production rate step change (−15%), 2. product mix step change (50 G/50H to 40 G/60H on mass basis), 3. reactor operating pressure step change (− 60 kPa), and 4. purge gas composition of component B step change ( +2 mol%). The selected controlled variables include the production rate, reactor pressure, and purge gas composition. Since the PI controller in the McAvoy and Ye scheme performs adequately for control of product quality (within the tolerance of 9 5 mol% G as suggested by Downs and Vogel), product quality was not included in the NMPC scheme. A fourth controlled variable, the A/C mole ratio in the reactor feed, was selected to address reaction stoichiometry disturbances. Disturbance IDV(1) is charac-
1191
terized by a step change in the A/C feed ratio in stream 4. This disturbance upsets the reaction stoichiometry causing the gaseous reactants to accumulate in the reactor–recycle loop with an associated increase in reactor pressure. Controlling the A/C feed ratio provides an effective means of stabilizing the reactor pressure. The NMPC scheme was implemented in a supervisory mode, where the manipulated variables were setpoints for the lower level PI controllers. The manipulated variables were selected based on the steady-state gain matrix presented by McAvoy and Ye (1994) along with the following observations. 1. A&C feed stream consists of reactant C and the bulk of reactant A, both of which are needed to form products, G and H. Thus, this stream directly affects the production rate. It is also the carrier of inert B and hence directly influences the purge gas composition of component B. 2. A decrease in reaction temperature quenches the reaction. In this case, the unreacted reactants, which are essentially noncondensibles, accumulate in the recycle loop causing the reactor pressure to rise. On the other hand, an increase in reaction temperature increases the rate of formation and subsequent vaporization of products. The vaporized products are then condensed leading to a drop in the recycle loop pressure. 3. Stream A has a much smaller throughput compared to A&C stream and supplies the remainder of component A needed for the reaction. Thus, this stream can be effectively used for control of the A/C ratio in the reactor feed. 4. Purge rate is a logical choice for control of component B in purge. It also affects the pressure by avoiding accumulation of noncondensibles. Thus, A feed rate, A&C feed rate, reactor temperature, and purge rate were selected as the set of manipulated variables for the NMPC controller. A summary of the resulting 4× 4 subset selected for NMPC is provided in Table 2. PI controller tuning parameter values used are as in McAvoy and Ye except for the NMPC subset and the two modifications noted previously. Note that use of a square system for the factorized RBF-based NMPC system is not required. The demonstration system is square by consequence rather than intent. The selected 4×4 system exhibits adequate complexity for MPC demonstration purposes.
4.2. De6elopment of the RBF model A sample interval of 3 min was used. Based on characteristics of the Eastman process, the prediction and control horizons were chosen as 1 h (p= 20) and 15 min (c =5), respectively.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1192
Table 3 Operating region represented in network training Manipulated variable
Maximum value
Minimum value
Reactor temperature (°C) Purge rate (kscmh) A&C feed rate (kscmh) A feed (kscmh)
116 0.05 7.5 0.05
128 0.85 11 1.05
Training data was generated by perturbing the NMPC manipulated variables around the base case. Table 3 gives the range of the perturbations. The magnitude of the steps was selected randomly within the ranges specified in Table 3. The duration of each step was selected randomly between 45 min and 2 h. A total of 35 runs with simulated process run times varying between 16 and 175 h were made. Each simulation was terminated when process conditions approached known operating constraints for the process (e.g. high pressure shutdown). Disturbances were turned off during each run. Before generating network input patterns, the data were normalized to zero mean and unity standard deviation. Measurements included in each network input pattern, xk are listed in Table 4. In addition to the control inputs, the input vector was augmented to include the key disturbance variable measurements of D feed, E feed and the total reactor feed rates. A total of 10 000 network input patterns were selected by uniformly choosing from the entire set of approximately 40 000 available patterns. The RBF network model was developed using 5000 input patterns.
The model was validated on the remaining 5000 patterns. Training and test set error statistics are provided in Table 5. The node centers and width were determined using methods discussed in Section 2.1. The model employed 230 hidden nodes with uniform widths, |, of 8.6. The quality of the RBF model is illustrated by the following experiment. At time zero, the A&C feed rate was stepped down from its base value of 9.3477 m3 h − 1 to 9.0 m3 h − 1 for a period of 2 h. The A&C feed rate was then decreased in steps of 0.5 m3 h − 1 at hourly intervals to 7.5 m3 h − 1. After 5 h, the A&C feed rate was stepped to 6.0 m3 h − 1. The final step lies outside the operating range used to generate the model (see Table 3). During this experiment, the A feed rate, purge rate and reactor temperature were maintained at their base values. A comparison of the actual and predicted behavior is shown in Fig. 6. As expected, good agreement is observed until after 5 h when the model begins to operate outside the range of the training data. While the absolute value of the model predictions degrade during this period, the model generally continues to correctly predict long term trends.
4.3. NMPC controller tuning Table 6 shows the error penalty and move suppression factors that were used. These values were established empirically using engineering judgment in the same manner as generated for industrial MPC applications. Hard constraints were implemented only on the manipulated variables. Table 3 documents the constraints used.
Table 4 Elements of RBF input pattern vector xk Delayed CVs Delayed MV Past input moves Disturbance inputs
Pressure,k−p A feed,k−p−1 DA feed,k−p, …, DA feed,k−1 D feed rate,k
% B in purge,k−p A&C feed,k−p−1 DC feed,k−p, …, DC feed,k−1 E feed rate,k
Product rate,k−p A/C mole ratio,k−p Purge rate,k−p−1 Temperature,k−p−1 DPurge,k−p, …, DPurge,k− DTemp.,k−p, …, DTemp.,k−1 1 Reactor feed rate,k
Table 5 Training and test set error statistics Reactor pressure Product flow rate (m3 Component B in purge (kPa) h−1) (mol%) Training set
Test set
Mean error −0.0024 Standard deviation of 8.3387 error Mean error 0.0042 Standard deviation of 8.7749 error
(A/C) mole ratio in reactor feed
−0.0001 0.1186
0.0000 0.1038
0.0000 0.0337
0.0006 0.1207
−0.0039 0.1089
0.0008 0.0331
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1193
Fig. 6. Comparison of RBF model predictions with plant measurements for step changes in A&C feed rate. All other variables are maintained at their base values. After 5 h, the A&C feed rate is brought to 6 m3 h − 1, which lies outside the lower limit of training data of 7.5 m3 h − 1. Table 6 Weights used in the RBF based MPC simulations Controlled variable
Error penalty
Manipulated variable
Move suppression
Reactor pressure Product flow rate Component B in purge A/C mole ratio in reactor feed
0.16 0.42 0.28 0.35
A feed setpoint C feed setpoint Purge rate Reactor temperature
0.7 1.0 0.8 1.4
5. Computation and application benefits
5.1. Computational benefits
In the simulation results reported below, setpoint changes or disturbances were implemented after 2 h of operation at the base case conditions. Further, all flow and pressure measurements were filtered using a CUSUM filter (Rhinehart, 1992). The MATLAB implementation of the Eastman process developed by Professor N.L. Ricker was used for all simulations. The MATLAB function constr was used to perform the constrained nonlinear minimization via sequential quadratic programming with active set method for constraint handling.
The speedup factor as described in Section 3 was determined as an estimate of the expected computational savings attributable only to model prediction efficiencies. For the 4-input, 4-output problem modeled using a 230 node network, with p= 20, c= 5, Ny =1, Nu = 2, and assuming n= 100, the calculated speedup factor is 7.05. Table 7 documents the time needed to compute the control relevant instructions for the factorized and nonfactorized NMPC controllers for the pressure setpoint change described below. In the non-factorized ap-
1194
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
Table 7 Comparison of computation time needed for implementation of a setpoint change of −60 kPa in reactor pressure (Fig. 8) (simulation for 1000 samples or 50 h) Factorized RBF Non-factorized RBF based NMPC based NMPC (gradients evaluated numerically) Real-time needed for computation (h)
1.82
38.31
proach, the gradients and Hessian were computed numerically. The factorized approach used analytical expressions for the gradient and numerical approximations for the Hessian (limitation of MATLAB function constr). The results presented in Table 7 were calculated using the tic-toc command in MATLAB and are typical for a 50 h simulation using a 550 MHz processor. As indicated by the results in Table 7, the factorized formulation required 21 times less computation effort than the non-factorized approach. The additional speedup beyond the predicted factor of 7 is attributable to the use of analytic expressions for the gradient. For the factorized formulation, the number of model prediction function evaluations needed for optimization convergence at each control step ranged from 3 to 15. For these values of n, the predicted speedup factor would be 3 – 5 rather than 7.
5.2. Application benefits Results for four of the test cases suggested by Downs and Vogel are presented to characterize performance of the NMPC algorithm. Fig. 7 shows simulation results for a step decrease in production rate by −15% at t= 2 h. The controller reduces the flowrates of A and A&C streams. The net decrease in the gaseous reactants causes the reactor pressure to initially drop approximately 100 kPa. The reaction rate is modified by a nominal 3 °C decrease in reactor temperature to meet the new production requirement and restore the reactor pressure. Fig. 8 shows results for the setpoint change in reactor pressure (− 60 kPa). The long-term response required an increase in the reactor temperature. The NMPC controller initially satisfied the new pressure setpoint by increasing the purge rate. As the reactor temperature slowly increased, the purge rate was returned to the original steady-state value and composition. The nonlinear process interactions were handled smoothly with minimal impact on production rate and quality. Results for a setpoint change in the composition of B (+ 2 mol%) in the purge are presented in Fig. 9. The desired response was achieved with negligible impact on production rate or quality through a combination of purge rate and reactor temperature adjustments. All three previous cases maintained operation in the region represented by the training data. Consequently,
Fig. 7. Product flowrate setpoint change ( − 15%).
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1195
Fig. 8. Reactor pressure setpoint change ( − 60 kPa).
the nonlinear characteristics of the process were well represented by the RBF model at all times with the result that the full benefits of MPC were realized. Fig. 10 shows the mismatch between plant measurements and RBF model predictions for the setpoint change in purge gas composition. As expected, process-model mismatch is essentially reduced to measurement noise. This result is typical for all three previously discussed cases. Disturbance rejection associated with IDV(1) (step change in A/C feed ratio in stream 4) and the corresponding process-model mismatch are shown in Figs. 11 and 12, respectively. The unmeasured step disturbance is first detected by a drop in a controlled vari-
able, the A/C molal ratio in the reactor feed (Fig. 13). Consequently, the A feed rate is increased to bring the A/C ratio in reactor feed back to setpoint. Since IDV(1) was a step disturbance, process-model mismatch ultimately reaches steady-state values (Fig. 12). Satisfactory control performance is realized since the NMPC algorithm employs plant-model mismatch as output feedback.
5.3. Discussion In addition to the ratio or speedup factor, the absolute values of the computation times required for the factorized and non-factorized approaches are impor-
1196
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
Fig. 9. Purge B composition setpoint change ( +2 mol%).
Fig. 10. Process-model mismatch during implementation of + 2% setpoint change in composition of component B in purge stream.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1197
Fig. 11. Disturbance IDV(1) (step change in A/C feed ratio in A&C stream).
tant. Use of a non-factorized controller that requires 38 + h of computation time during 50 h of plant operation is unrealistic. The factorized algorithm completed the control calculations in 1.82 h or 3.6% of the total run time. In this case, the difference in computation costs has a major impact on the practicality of utilizing an RBF-based NMPC control algorithm. The performance results presented are expected based on the nonlinear modeling capabilities of RBF neural networks. Hence, an NMPC scheme that employs an RBF model would be expected to provide good control whenever the plant operated in regions used to train the model. Practical implementation of an RBF (or any empirical model) based MPC scheme
would require a watchdog to verify that operation lies within model development bounds. The problem/concept is the same as required for the linear MPC systems currently employed in industry. However, the level of concern is greater for an RBF based system due to the less certain (compared to a linear impulse or step response model) extrapolation characteristics. Another issue affecting practical implementation of an RBF based NMPC scheme is collection of the plant data necessary to develop the model. The amount of data required to develop any nonlinear model is orders of magnitude greater than required for a simple linear model. However, the time and investment required to develop the linear models used with today’s linear MPC
1198
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
Fig. 12. Process-model mismatch during occurrence of disturbance IDV(1). The unmeasured disturbance is not modeled by the RBF network, leading to poor predictions. However, due to the step nature of the disturbance, the process-model mismatches settle to near-steady values.
systems is already an impediment in many cases. One area of research that can potentially address this concern is model identification using closed-loop data. Falcioni and Seborg (1997) have studied this issue using neural network models. Bomberger et al. (1998) used historical data to identify RBF models for a copolymerization reactor. Potential improvement in performance using nonlinear control methods provides the incentive to find new ways to develop robust RBF models from existing closed-loop plant data.
6. Conclusions From a manufacturing standpoint, the development of faster and more powerful computers provides the opportunity to apply the concepts of MPC with process models of almost any type. Computationally intensive modeling methods such as RBF networks can now be realistically considered to achieve the benefits of nonlinear MPC. However, good engineering practice demands efficient utilization of all resources, including computation. Methods to improve the computational efficiency of new or improved control algorithms as addressed in
Fig. 13. The step disturbance IDV(1) is detected by the drop in the controlled variable, A/C ratio in reactor feed, from its setpoint. The variable is then brought back to its setpoint by increasing A feed rate.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
this paper are essential to accelerate potential acceptance and implementation in practice. Among the issues affecting potential application of RBF-based NMPC, model identification and maintenance remain the most prominent. As discussed previously, these issues represent hurdles for essentially all forms of NMPC. The highly distributed, non-compact form of RBF models is undesirable from a maintenance standpoint where simplicity and parsimony are highly valued. Ultimately, the nonlinear modeling capabilities of RBF networks provide the economic incentive to overcome this and all other impediments to industrial implementation. Future work on the factorized RBF-based NMPC technique will progress along two lines: identification and control. In the identification phase, development of RBF models using closed-loop data will be stressed. Also the broader question of input sequence design for nonlinear system identification will be explored. The control phase will expand the current version of NMPC with output feedback to optimal state estimation procedures.
Acknowledgements The authors express their gratitude to Professor N.L. Ricker and his group for use of their MATLAB simulation of the Eastman process.
References Banerjee, A., Arkun, Y., Ogunnaike, B., & Pearson, R. (1997). Estimation of nonlinear systems using linear multiple models. American Institute of Chemical Engineering Journal, 43, 1204 – 1226. Bhartiya, S., & Whiteley, J. R. (2001). Factorized approach to nonlinear MPC using a radial basis function model. American Institute of Chemical Engineering Journal, 47, 358–368. Bomberger, J. D., & Seborg, D. E. (1998). Determination of model order for NARX models directly from input –output data. Journal of Process Control, 8, 459–468. Bomberger, J. D., Seborg, D. E., & Ogunnaike, B. A. (1998). Radial basis function identification of a solution copolymerization reactor model (pp. 3182 – 3188). Philadelphia, PA: ACC. Chen, S., Billings, S. A., & Grant, P. M. (1992). Recursive hybrid algorithm for non-linear system identification using radial basis function networks. International Journal of Control, 55, 1051 – 1070.
1199
Cybenko, G. (1987). Approximation by superposition of a sigmoidal function. Mathematical Control Signal Systems, 2, 303. Downs, J. J., & Vogel, E. F. (1993). A plant-wide industrial process control problem. Computers and Chemical Engineering., 17, 245– 255. Falcioni, M., & Seborg, D. E. (1997). On the identification of neural network models during closed-loop operation. In European control conference. Brussels. Garcia, C. E. (1984). Quadratic/dynamic matrix control of nonlinear processes: An application to batch reaction process. In AIChE annual meeting. San Francisco, CA. Gattu, G., & Zafiriou, E. (1995). Observer based nonlinear quadratic dynamic matrix control for state space and input/output models. The Canadian Journal of Chemical Engineering, 73, 883 –895. Gurumoorthy, A., & Kosanovich, K. A. (1998). Improving the prediction capability of radial basis function networks. Industrial Engineering and Chemical Research, 37, 3956 – 3970. Hartman, E. J., Keeler, J. D., & Kowalski, J. M. (1990). Layered neural networks with Gaussian hidden units as universal approximations. Neural Computation, 2, 210 – 215. He, X., & Asada, H. (1993). A new method for identifying orders of input – output models for nonlinear dynamic systems (pp. 2520 – 2523). San Francisco, CA: ACC. Henrique, H., Lima, E. L., & Seborg, D. E. (2000). Model structure determination in neural network models. Chemical Engineering Science, 55, 5457 – 5469. Hush, D. R., & Horne, B. G. (1993). Progress in supervised neural networks; What’s new since Lippmann, IEEE Signal Processing Magazine, 10, 8 – 39. Hussain, M. A. (1999). Review of the applications of neural networks in chemical process control-simulation and online implementation. Artificial Intelligence in Engineering, 13, 55 – 68. Kennel, M. B., Brown, R., & Abarbanel, H. D. I. (1992). Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical Re6iew, 45, 3403 – 3411. Krishnan, A., & Kosanovich, K. (1998). Batch reactor control using a multiple model-based controller design. The Canadian Journal of Chemical Engineering, 76, 806 – 815. Lee, J. H., & Ricker, N. L. (1994). Extended Kalman filter based nonlinear model predictive control. Industrial Engineering Chemical Research, 33, 1530 – 1541. McAvoy, T. J., & Ye, N. (1994). Base control for the Tennessee Eastman problem. Computers and Chemical Engineering, 18, 383– 413. Moody, J., & Darken, C. J. (1989). Fast learning in networks of locally-tuned processing units. Neural Computation, 1, 281–294. Poggio, T., & Girosi, F. (1990). Networks for approximation and learning. Proceedings of the IEEE, 78, 1481 – 1497. Pottman, M., & Seborg, D. E. (1992). Identification of nonlinear processes using reciprocal muliquadric functions. Journal of Process Control, 2, 189 – 203. Rhinehart, R. R. (1992). A CUSUM type on-line filter. Process Control and Quality, 2, 169 – 176. Ricker, N. L. (1996). Decentralized control of the Tennessee Eastman challenge process. Journal of Processing Control, 6, 205 –221. Sanner, R. M., & Slotine, J. J. (1992). Gausian networks for direct adaptive control. IEEE Transactions on Neural Networks, 3, 837– 863.
Benefits of factorized RBF-based NMPC Sharad Bhartiya, James R. Whiteley * School of Chemical Engineering, Oklahoma State Uni6ersity, 423 Engineering North, Stillwater, OK 74078 -5021, USA Received 13 June 2001; received in revised form 6 September 2001; accepted 25 February 2002
Abstract The computation and application benefits of factorized RBF-based nonlinear model predictive control (NMPC) are presented. The NMPC algorithm derives its computational efficiency by factorizing the radial basis function (RBF) model response. A brief description of the factorized RBF-based NMPC algorithm is provided. Theoretical computation benefits are quantified for both SISO and MIMO formulations. Computation and application benefits of the algorithm are documented for a 4 ×4 MIMO subset of the Eastman Challenge problem. Results confirm computation requirements are reduced by more than an order of magnitude relative to application of traditional MPC using a non-factorized RBF model. The expected control performance benefits from using a nonlinear process model are also achieved. © 2002 Elsevier Science Ltd. All rights reserved. Keywords: Nonlinear MPC; RBF model; Multivariable control
1. Introduction The current generation of commercially available model predictive control (MPC) packages relies on linear process models to provide future predictions. Excellent performance is realized as long as the plant operates close to the conditions used to create the linear approximation of the process. The goal for the next generation of control software is to provide similar capability across the whole range of possible plant operating conditions. The notion of using a nonlinear model within the MPC paradigm has led to an active interest in the development of nonlinear model predictive control (NMPC) methods. A common approach to the nonlinear optimization required in NMPC is based on successive linearization of nonlinear models. Garcia (1984) proposed a nonlinear QDMC algorithm, a simple extension of DMC/ QDMC based on online successive linearization of a mechanistic nonlinear model. Nonlinear MPC using closed-loop state estimation by an extended Kalman filter has been proposed by Lee and Ricker (1994). Gattu and Zafiriou (1995) augmented the system states with stochastic states to account for modeling errors * Corresponding author. Fax: + 1-405-744-6338. E-mail address: [email protected] (J.R. Whiteley).
and disturbances. Banerjee, Arkun, Ogunnaike, and Pearson (1997) describe a method of state estimation for nonlinear systems that are subject to multiple operation regimes and make transitions between them. The nonlinear process is approximated by a linear parameter varying system which consists of local linear models. Krishnan and Kosanovich (1998) also present a multiple model based MPC scheme. The linear time invariant models are computed offline along a pre-defined reference trajectory of a batch process. The NMPC method addressed by this paper utilizes an artificial neural network (ANN) for the process model. The nonlinear modeling capability of ANNs is well documented (Cybenko, 1987; Hartman, Keeler, & Kowalski, 1990). Multi-layer perceptron and radial basis function (RBF) networks are the most extensively utilized ANNs in control and identification applications. Hussain (1999) provides a summary of a number of control applications reported in the literature. The current work employs an RBF network to model the process. Development of process models using RBF networks has been extensively explored by Seborg and co-workers. Pottman and Seborg (1992) applied an RBF based predictive control strategy to an experimental pH neutralization process. Bomberger, Seborg, and Ogunnaike (1998) studied identification of a multivariable copolymerization model with emphasis on design
0098-1354/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 2 ) 0 0 0 2 9 - 7
1186
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
of input sequence. Gurumoorthy and Kosanovich (1998) incorporated process knowledge via regularization theory to improve the prediction capability of RBF networks. Also, online adaptation of RBF networks has been studied. Chen, Billings, and Grant (1992) used a hybrid clustering algorithm along with recursive least squares to update the network centers and weights, respectively. Sanner and Slotine (1992) applied model update rules for RBF networks in an adaptive control framework. Because of the computational cost associated with nonlinear optimization, there are practical incentives to integrate a nonlinear model in an MPC scheme in the most efficient manner possible. RBF networks provide good modeling capability but are not ‘‘compact’’ from a traditional serial computing standpoint. RBF models for nontrivial processes require a large number of parameters and involve many more mathematical operations than time series and other common types of models. The factorized RBF-based NMPC method (Bhartiya & Whiteley, 2001) was developed to reduce the computational penalty associated with repeated RBF model predictions during nonlinear optimization searches. The remainder of this paper is organized as follows. The essential points for the factorized RBF-based NMPC algorithm are presented in Section 2. The computation benefits for both SISO and MIMO formulations of the algorithm are discussed in Section 3. For illustrative purposes, the algorithm is applied to a 4× 4 MIMO subset of the Eastman Challenge problem (Downs & Vogel, 1993). An overview of the Eastman process and development of the NMPC controller are described in Section 4. Benefits associated with the NMPC controller are presented in Section 5 followed by conclusions in Section 6.
2. Factorized RBF-based nonlinear MPC algorithm Feedforward RBF networks approximate process dynamics by using past process information as network inputs. Thus, yˆk = F(yk − 1, …, yk − Ny, uk − 1, …, uk − Nu )
(1)
where a time-delay of unity is assumed between the model output, yˆk, and the previous process inputs, uk − 1. Ny and Nu are integers defined by the order of the model. Scalars uk − 1, …, uk − Nu represent the sequence of inputs used by the model. Function F depends on the network architecture and the type of activation function employed by the nodes. Model orders, Ny and Nu, may be determined by methods discussed by He and Asada (1993) (based on Lipschitz numbers), and Kennel, Brown, and Abarbanel (1992) (using false nearest neighbors). Bomberger and Seborg (1998) demonstrate
use of these methods in conjunction with an RBF network. Bhartiya and Whiteley (2001) proposed a factorized RBF based NMPC algorithm that exploits the factorizability property of the Gaussian function yielding a numerically efficient algorithm. Model identification is achieved via training of an RBF network using Gaussian nodes to provide sequential predictions over a prediction horizon of length p. The NMPC controller is directly parameterized in terms of the RBF network parameters. This allows analytical expressions for the gradient and Hessian of the objective function.
2.1. Model formulation The RBF network consists of an input layer, a hidden layer, and an output layer. In the input layer, unweighted inputs, x, are directly transmitted to the hidden layer nodes. Each hidden layer node consists of a RBF. In the current work, the Gaussian function is employed,
gj (x) =exp −
x− tj 2 |2
(2)
where tj is the center of the jth hidden node. It is assumed that each node is of fixed width |. The output layer performs a weighted summation of hidden node outputs to give the network output, m
yˆi (x) = % wij gj (x)
(3)
j=1
where m represents the number of hidden nodes and wij is the weight between the ith output node and the jth hidden node. The node centers, tj, are fixed in an unsupervised training step by partitioning the input data space using a clustering algorithm (e.g. k-means). Node width, |, accounts for the desired degree of overlap between the clusters and is fixed using heuristics (e.g. mean distance between cluster centers) (Hush & Horne, 1993; Moody & Darken, 1989). Finally, the network weights are obtained by linear regression to approximate the input–output mapping provided by the data. The number of nodes, m, directly contributes to the model complexity (number of free parameters). Henrique, Lima, and Seborg (2000) discuss node pruning techniques based on orthogonal least squares which yield parsimonious models. Often, the training data consists of noisy measurements. This may cause the minimization of prediction error to be an ill-posed problem. Regularization theory can be used to address this problem by incorporating smoothness constraints or a priori knowledge (Poggio & Girosi, 1990). The RBF model in the proposed NMPC scheme provides non-iterative sequential predictions over a prediction horizon of length p by avoiding dependency of future model predictions on previous model predic-
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
tions. Future predictions, yˆ k, are related to past measurements yi and inputs ui, delayed by p samples and input moves, Dui, as follows, yˆ k/k − p = F(yk − p, …, yk − p + 1 − Ny, uk − p − 1, …, uk − p + 1 − Nu, Duk − p, …, Duk − 1)
(4)
where the control move at the kth instant, Duk = uk − uk − 1
(5)
The delayed input/output measurements, viz. yk − p, …, yk − p + 1 − Ny and uk − p − 1, …, uk − p − Nu, respectively, provide a reference to the state of the system p samples in the past. Based on the model in Eq. (3), the following form of the input vector is chosen for a single-input, single-output system, xk = [yk − p …yk − p + 1 − Ny uk − p − 1…uk − p + 1 − Nu Duk − p…Duk − 1]
(6)
Fig. 1 illustrates the input vector with elements located on a timeline. Measurable disturbances can be accounted by augmenting the input vector x to include the disturbance variables. The RBF prediction of the process output at instant k is then obtained by Eq. (3). The form is general and can accommodate both multiple inputs and outputs.
1187
2.2. NMPC controller formulation The algorithm computes the manipulated variable profile over a control horizon by optimizing an objective function defined over the prediction horizon, subject to constraints. Only the first move is implemented and the procedure is repeated at every sampling instant. The following objective function is used: p
c−1
i=1
i=0
= % Yi (rk + i − yˆ k + i )2 + % \i (Duk + i )2
(7)
Variables p and c represent the prediction and control horizons, respectively. Yi and \i denote the error penalty and move suppression factors at the ith instant. The MPC control law can be stated as, arg(min ) Du , Du ,…,Du k k+1 k + c − 1 such that ymin 5 yˆ k + i 5 ymax (8) Dumin 5 Duk + i 5 Dumax umin 5 uk + i 5 umax Future model predictions, yˆ k + i, depend on past control moves and the future control move variables, Duk + i (see Eq. (4)). The future control moves, Duk, …, Duk + c − 1, represent the decision variables for the optimization problem in Eq. (8). For calculation purposes, it is desirable to express yˆ k + i such that the unknown decision variables appear explicitly in the objective function. The key idea in the factorized RBF
Fig. 1. Timelines showing inputs to the p-step control model. For this example, p = 4, Ny =3, and Nu =2. Predicted outputs are generated from known information only, previous model predictions for y are not used in model input.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1188
model lies in expressing the model prediction, yˆ k + i, as an inner product of two vectors. The unknown decision variables of the nonlinear program (Eq. (7)) are contained in one vector and all other known past quantities, including the network weights in the other. Thus, the RBF output can be rearranged as follows: m1
yˆ k + i = % wj exp(past +future) j=1
Æ w1 exp(past) ÇT Æexp(future) Ç Ã Ã Ã Ã =Ã Ã Ã Ã Ã Ã Ã Ã Èwm 1 exp(past) É Èexp(future) É
(9)
yˆ k + i =yˆ p,k +Tiyˆ f,k + i
(10)
Subscripts p and f refer to the fact that the corresponding factors contain all known (past) and unknown (future) terms, respectively. Thus, only yˆ f,k + i needs to be computed during every function call by the optimization algorithm. For a SISO problem, the factors for the predictions, yˆ k + 1, …, yˆ k + p, based on c future moves, Duk, …, Duk + c − 1, take the following form:
yˆ p,k + i =
w exp −
1 |2
Ny − 1
% (yk − p − j + i 1
j=0
Nu − 1
−tyj + 1)2 + % (uk − p − j + i 1 −tuj )2
n
j=1 p−1
+ % (Duk − p + j + i 1 − tj + 1)2
(11)
j=0
The symbol 1 represents a column vector of size m with unity elements. The operator ‘‘ ’’ is used to denote the Hadamard product (element by element multiplication). Note that when the index, i, equals p, the final summation drops out since j varies from 0 to − 1. The factor yˆ f,k + i that contains the future moves the optimization algorithm must calculate is:
yˆ f,k + i = exp −
i−1
i= 1, …, c yˆ f,k + i = exp −
n
1 2 % (Du 1 − tDu , p − j) | 2 j=0 k+i−j−1 i−1
n
(12a)
1 2 % (Duk + i − j − 11 − tDu , p − j) | 2 j=i−c
i=c+ 1, …, p
the gradient and Hessian of the objective function. The availability of analytic expressions significantly reduces the computational requirements associated with the algorithm. In addition, the separation of the decision variables in the model prediction ensures that only the unknown parts of the objective function, and the gradient and Hessian required by the sequential quadratic programming (SQP) algorithm are recalculated during optimization.
3. Computational efficiency To quantify the computational efficiency of the factorized approach, we first consider the number of floating point operations (flops) required when the SQP solver calls the function to produce the future p model predictions yˆ k + i required by the controller objective function (Eq. (7)). The total number of flops depends on, (a) the formulation of the RBF model, factorized or non-factorized, (b) the number of nodes in the RBF network, m, (c) the size of the control problem (the number of control outputs, No, and control inputs, Ni ), (d) the orders of the model, Ny and Nu, and (e) the length of the prediction and control horizons, p and c, respectively. Counting addition, subtraction, multiplication, division, and exponentiation as single floating point operations, the total number of flops required to evaluate vectors yˆ p,k + i and yˆ f,k + i, can be determined directly from Eqs. (11) and (12), respectively. For the SISO case, it is straightforward to show that the number of flops, Cp,SISO, required to generate yˆ p,k + i is, Cp,SISO = pm(3Ny + 3Nu + 3p −2)
Likewise, the number of flops, Cf,SISO, required to generate yˆ f,k + i is, Cf,SISO = mc(c+ 2)+ m(p −c)(2c + 3)
The exponential function in Eqs. (11) and (12) implies a term by term application to each element of the column vector in the square brackets. Column vector tj is constructed by using the jth element of all m centers. Thus, any future prediction, yˆ k + i, can be computed using Eqs. (10)–(12). Extension to the MIMO case is straightforward (Bhartiya & Whiteley, 2001). The factorized form provides analytic expressions for
(14)
Analogously, for the Ni × No MIMO problem, assuming each input and output have the same associated orders, Ny and Nu, respectively, the following results are obtained, Cp,MIMO = pm(3+3NoNy + 3Ni Nu + 3Ni p− 5Ni )
(12b)
(13)
(15)
and Cf,MIMO = mc[Ni (c− 1)+ 3]+ m(p−c)(2Ni c+ 3) (16) For every function call when using the factorized model, calculation of the p future model predictions using Eq. (10) requires, Cfactorized = Cp + Cf + p(2m −1)
(17)
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1189
However, the model prediction function is called numerous times at every control step. Thus, if n total function calls are made by the SQP solver for computation of the control inputs, the factorized approach will need, Ctotal,factorized = Cp + nCf + np(2m − 1)
(19)
total flops (since Cp is only calculated once), whereas the non-factorized model predictions will require, Ctotal,non-factorized = n(Cp + Cf − 3mp) Fig. 2. Computational speedup associated with the factorized RBF model increases with number of model predictions required at each time step by the NMPC optimization algorithm. Plot prepared with m = 230, p = 20, c= 5, Ni = 4, No = 4, Nu = 2 and Ny = 1.
flops. The reduction in computational expense associated with the factorized RBF model is difficult to appreciate from inspection of Eqs. (13)–(20). A sensitivity analysis using typical variable values was performed to characterize the benefits of model factorization in terms of a speedup factor defined as, Speedup factor=
Fig. 3. Contour plot of speedup factor as a function of the prediction and control horizons. Plotted for m = 230, n = 100, Ni = 4, No =4, Nu = 2 and Ny =1.
Fig. 4. Contour plot of speedup factor as a function of the number of MIMO inputs and outputs. Plotted for m= 230, n =100, Ni =4, No = 4, Nu =2 and Ny = 1.
floating point operations. The last term in Eq. (17) is associated with the operations performed to calculate the inner product in Eq. (10) (m products followed by (m− 1) additions). For the non-factorized model, the number of flops needed for calculation of p model predictions is,
Cnon-factorized = (Cp +Cf −3mp)
(18)
(20)
Ctotal,non-factorized Ctotal,factorized
(21)
A speedup factor of 10 indicates an order of magnitude fewer computations are required to produce model predictions using a factorized RBF model. The sensitivity analysis reveals that the speedup factor is independent of the number of nodes, m. The speedup factor is dependent only on n, p, c, Ni, No, Nu, and Ny and not on the size of the RBF model. The speedup factor increases in a first-order manner with the total number of model predictions or function calls, n, required to solve the MPC control law, Eq. (8). Fig. 2 was generated for a 4×4 MIMO system with prediction and control horizons of 20 and 5, respectively, input and output model orders of 2 and 1, respectively, and an RBF model with 230 hidden nodes. The trends and speedup values shown in Fig. 2 (as well as Figs. 3 and 4 to follow) are representative for a wide range of variable values typical for industrial MPC applications. As indicated in Fig. 2, a five-fold improvement is realized after only 13 function calls. As the number of function calls increases, the speedup factor approaches an asymptotic limit of approximately seven. The number of function calls depends on the efficiency of the nonlinear optimization search. A contour plot of the speedup factor as a function of the lengths of the prediction and control horizons, p and c, respectively, is presented in Fig. 3. The variable values used to generate the plot are the same as those for Fig. 2 with n fixed at 100. Typical practice is to implement MPC with p/c ratios ranging from 3/1 to 5/1. As revealed in Fig. 3, the speedup factor is nearly constant for a fixed p/c ratio and varies from 6 to 7 for typical MPC applications. The size of the MIMO control problem impacts the speedup factor through Ni (number of inputs or MVs)
1190
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
and No (number of outputs or CVs). As seen in Fig. 4, all combinations with between 3– 15 inputs and 3–20 outputs produce speedup factors between 6.8 and 7.8. Fig. 4 was generated in the same manner as Figs. 2 and 3 with n=100 and a p/c ratio of 4 with p =20. As one would expect, the speedup factor increases with increasing complexity or size of the MIMO problem. A review of all results generated by the different sensitivity analyses indicates a speedup factor of approximately seven can be expected for a typical MIMO control problem. However, there is an additional computational benefit associated with the factorized RBF model. Analytical expressions can be provided for the gradient and Hessian of the objective function. Significant computational savings are realized by elimination of the need for numerical approximations. With cNi decision variables in the optimization program in Eq. (8), a single numerical approximation of the gradient requires cNi additional function calls. Likewise, a single numerical approximation of the Hessian requires (cNi )2 calls. The actual computational savings are magnified based on the need for repeated evaluations of the gradient and Hessian to satisfy convergence criteria for both the numerical approximations (inner computation loop) and overall optimization search (outer computation loop).
4. Application of RBF-based NMPC to the Eastman problem The Eastman challenge problem (Downs & Vogel, 1993) was employed to characterize the computational and application benefits of the factorized RBF-based NMPC algorithm. The process (Fig. 5) involves producing two liquid products, G and H, from four gaseous reactants, A, C, D, and E.
4.1. Plantwide control strategy A variety of decentralized and MPC control strategies have been proposed for the Eastman process. Ricker (1996) discusses the key issues involved with design of the plantwide control system with references to many of the previously published designs. The current work is not intended to provide another comprehensive solution to the Eastman problem. Rather, the control system described below was developed strictly to examine performance characteristics of the proposed NMPC algorithm on a meaningful scale. The SISO strategy previously reported by McAvoy and Ye (1994) was used for the base regulatory control system. The SISO loop pairings recommended by McAvoy and Ye are presented in Table 1. Two minor modifications were made to McAvoy and Ye’s scheme:
Fig. 5. Schematic of the Eastman challenge process.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199 Table 1 Input–output pairing determined by McAvoy and Ye (1994) Controlled variable
Manipulated variable
Reactor temperature Reactor pressure Reactor level Separator level Stripper temperature Stripper level Production rate G/H composition in product E composition in product Inert B composition in purge stream Recycle flow
Reactor coolant temperature setpoint A feed rate setpoint E feed rate setpoint Underflow rate setpoint Steam flow setpoint Sripper underflow rate setpoint A&C feed rate setpoint D/E ratio setpoint
Compressor power
Stripper temperature setpoint Purge flow rate setpoint Condenser cooling water temperature setpoint Recycle valve
Table 2 Controlled and manipulated variables for the RBF based MPC controller MPC controlled variables (CVs)
MPC manipulated variables (MVs)
Reactor pressure Composition of B in purge Production rate A/C mole ratio in reactor feed
Reactor temperature setpoint Purge rate setpoint A&C feed rate setpoint A feed rate setpoint
1. The setpoint of the condenser cooling water return is used to control the condenser temperature (Gain=5.0%/°C; integral time=50 min). 2. The flow controller for the A&C feed stream (stream 4) was tuned more aggressively (Gain= 0.3%/kscmh; integral time=0.09 min). Controlled variables were selected to meet the following servo specifications suggested by Downs and Vogel: 1. production rate step change (−15%), 2. product mix step change (50 G/50H to 40 G/60H on mass basis), 3. reactor operating pressure step change (− 60 kPa), and 4. purge gas composition of component B step change ( +2 mol%). The selected controlled variables include the production rate, reactor pressure, and purge gas composition. Since the PI controller in the McAvoy and Ye scheme performs adequately for control of product quality (within the tolerance of 9 5 mol% G as suggested by Downs and Vogel), product quality was not included in the NMPC scheme. A fourth controlled variable, the A/C mole ratio in the reactor feed, was selected to address reaction stoichiometry disturbances. Disturbance IDV(1) is charac-
1191
terized by a step change in the A/C feed ratio in stream 4. This disturbance upsets the reaction stoichiometry causing the gaseous reactants to accumulate in the reactor–recycle loop with an associated increase in reactor pressure. Controlling the A/C feed ratio provides an effective means of stabilizing the reactor pressure. The NMPC scheme was implemented in a supervisory mode, where the manipulated variables were setpoints for the lower level PI controllers. The manipulated variables were selected based on the steady-state gain matrix presented by McAvoy and Ye (1994) along with the following observations. 1. A&C feed stream consists of reactant C and the bulk of reactant A, both of which are needed to form products, G and H. Thus, this stream directly affects the production rate. It is also the carrier of inert B and hence directly influences the purge gas composition of component B. 2. A decrease in reaction temperature quenches the reaction. In this case, the unreacted reactants, which are essentially noncondensibles, accumulate in the recycle loop causing the reactor pressure to rise. On the other hand, an increase in reaction temperature increases the rate of formation and subsequent vaporization of products. The vaporized products are then condensed leading to a drop in the recycle loop pressure. 3. Stream A has a much smaller throughput compared to A&C stream and supplies the remainder of component A needed for the reaction. Thus, this stream can be effectively used for control of the A/C ratio in the reactor feed. 4. Purge rate is a logical choice for control of component B in purge. It also affects the pressure by avoiding accumulation of noncondensibles. Thus, A feed rate, A&C feed rate, reactor temperature, and purge rate were selected as the set of manipulated variables for the NMPC controller. A summary of the resulting 4× 4 subset selected for NMPC is provided in Table 2. PI controller tuning parameter values used are as in McAvoy and Ye except for the NMPC subset and the two modifications noted previously. Note that use of a square system for the factorized RBF-based NMPC system is not required. The demonstration system is square by consequence rather than intent. The selected 4×4 system exhibits adequate complexity for MPC demonstration purposes.
4.2. De6elopment of the RBF model A sample interval of 3 min was used. Based on characteristics of the Eastman process, the prediction and control horizons were chosen as 1 h (p= 20) and 15 min (c =5), respectively.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1192
Table 3 Operating region represented in network training Manipulated variable
Maximum value
Minimum value
Reactor temperature (°C) Purge rate (kscmh) A&C feed rate (kscmh) A feed (kscmh)
116 0.05 7.5 0.05
128 0.85 11 1.05
Training data was generated by perturbing the NMPC manipulated variables around the base case. Table 3 gives the range of the perturbations. The magnitude of the steps was selected randomly within the ranges specified in Table 3. The duration of each step was selected randomly between 45 min and 2 h. A total of 35 runs with simulated process run times varying between 16 and 175 h were made. Each simulation was terminated when process conditions approached known operating constraints for the process (e.g. high pressure shutdown). Disturbances were turned off during each run. Before generating network input patterns, the data were normalized to zero mean and unity standard deviation. Measurements included in each network input pattern, xk are listed in Table 4. In addition to the control inputs, the input vector was augmented to include the key disturbance variable measurements of D feed, E feed and the total reactor feed rates. A total of 10 000 network input patterns were selected by uniformly choosing from the entire set of approximately 40 000 available patterns. The RBF network model was developed using 5000 input patterns.
The model was validated on the remaining 5000 patterns. Training and test set error statistics are provided in Table 5. The node centers and width were determined using methods discussed in Section 2.1. The model employed 230 hidden nodes with uniform widths, |, of 8.6. The quality of the RBF model is illustrated by the following experiment. At time zero, the A&C feed rate was stepped down from its base value of 9.3477 m3 h − 1 to 9.0 m3 h − 1 for a period of 2 h. The A&C feed rate was then decreased in steps of 0.5 m3 h − 1 at hourly intervals to 7.5 m3 h − 1. After 5 h, the A&C feed rate was stepped to 6.0 m3 h − 1. The final step lies outside the operating range used to generate the model (see Table 3). During this experiment, the A feed rate, purge rate and reactor temperature were maintained at their base values. A comparison of the actual and predicted behavior is shown in Fig. 6. As expected, good agreement is observed until after 5 h when the model begins to operate outside the range of the training data. While the absolute value of the model predictions degrade during this period, the model generally continues to correctly predict long term trends.
4.3. NMPC controller tuning Table 6 shows the error penalty and move suppression factors that were used. These values were established empirically using engineering judgment in the same manner as generated for industrial MPC applications. Hard constraints were implemented only on the manipulated variables. Table 3 documents the constraints used.
Table 4 Elements of RBF input pattern vector xk Delayed CVs Delayed MV Past input moves Disturbance inputs
Pressure,k−p A feed,k−p−1 DA feed,k−p, …, DA feed,k−1 D feed rate,k
% B in purge,k−p A&C feed,k−p−1 DC feed,k−p, …, DC feed,k−1 E feed rate,k
Product rate,k−p A/C mole ratio,k−p Purge rate,k−p−1 Temperature,k−p−1 DPurge,k−p, …, DPurge,k− DTemp.,k−p, …, DTemp.,k−1 1 Reactor feed rate,k
Table 5 Training and test set error statistics Reactor pressure Product flow rate (m3 Component B in purge (kPa) h−1) (mol%) Training set
Test set
Mean error −0.0024 Standard deviation of 8.3387 error Mean error 0.0042 Standard deviation of 8.7749 error
(A/C) mole ratio in reactor feed
−0.0001 0.1186
0.0000 0.1038
0.0000 0.0337
0.0006 0.1207
−0.0039 0.1089
0.0008 0.0331
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1193
Fig. 6. Comparison of RBF model predictions with plant measurements for step changes in A&C feed rate. All other variables are maintained at their base values. After 5 h, the A&C feed rate is brought to 6 m3 h − 1, which lies outside the lower limit of training data of 7.5 m3 h − 1. Table 6 Weights used in the RBF based MPC simulations Controlled variable
Error penalty
Manipulated variable
Move suppression
Reactor pressure Product flow rate Component B in purge A/C mole ratio in reactor feed
0.16 0.42 0.28 0.35
A feed setpoint C feed setpoint Purge rate Reactor temperature
0.7 1.0 0.8 1.4
5. Computation and application benefits
5.1. Computational benefits
In the simulation results reported below, setpoint changes or disturbances were implemented after 2 h of operation at the base case conditions. Further, all flow and pressure measurements were filtered using a CUSUM filter (Rhinehart, 1992). The MATLAB implementation of the Eastman process developed by Professor N.L. Ricker was used for all simulations. The MATLAB function constr was used to perform the constrained nonlinear minimization via sequential quadratic programming with active set method for constraint handling.
The speedup factor as described in Section 3 was determined as an estimate of the expected computational savings attributable only to model prediction efficiencies. For the 4-input, 4-output problem modeled using a 230 node network, with p= 20, c= 5, Ny =1, Nu = 2, and assuming n= 100, the calculated speedup factor is 7.05. Table 7 documents the time needed to compute the control relevant instructions for the factorized and nonfactorized NMPC controllers for the pressure setpoint change described below. In the non-factorized ap-
1194
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
Table 7 Comparison of computation time needed for implementation of a setpoint change of −60 kPa in reactor pressure (Fig. 8) (simulation for 1000 samples or 50 h) Factorized RBF Non-factorized RBF based NMPC based NMPC (gradients evaluated numerically) Real-time needed for computation (h)
1.82
38.31
proach, the gradients and Hessian were computed numerically. The factorized approach used analytical expressions for the gradient and numerical approximations for the Hessian (limitation of MATLAB function constr). The results presented in Table 7 were calculated using the tic-toc command in MATLAB and are typical for a 50 h simulation using a 550 MHz processor. As indicated by the results in Table 7, the factorized formulation required 21 times less computation effort than the non-factorized approach. The additional speedup beyond the predicted factor of 7 is attributable to the use of analytic expressions for the gradient. For the factorized formulation, the number of model prediction function evaluations needed for optimization convergence at each control step ranged from 3 to 15. For these values of n, the predicted speedup factor would be 3 – 5 rather than 7.
5.2. Application benefits Results for four of the test cases suggested by Downs and Vogel are presented to characterize performance of the NMPC algorithm. Fig. 7 shows simulation results for a step decrease in production rate by −15% at t= 2 h. The controller reduces the flowrates of A and A&C streams. The net decrease in the gaseous reactants causes the reactor pressure to initially drop approximately 100 kPa. The reaction rate is modified by a nominal 3 °C decrease in reactor temperature to meet the new production requirement and restore the reactor pressure. Fig. 8 shows results for the setpoint change in reactor pressure (− 60 kPa). The long-term response required an increase in the reactor temperature. The NMPC controller initially satisfied the new pressure setpoint by increasing the purge rate. As the reactor temperature slowly increased, the purge rate was returned to the original steady-state value and composition. The nonlinear process interactions were handled smoothly with minimal impact on production rate and quality. Results for a setpoint change in the composition of B (+ 2 mol%) in the purge are presented in Fig. 9. The desired response was achieved with negligible impact on production rate or quality through a combination of purge rate and reactor temperature adjustments. All three previous cases maintained operation in the region represented by the training data. Consequently,
Fig. 7. Product flowrate setpoint change ( − 15%).
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1195
Fig. 8. Reactor pressure setpoint change ( − 60 kPa).
the nonlinear characteristics of the process were well represented by the RBF model at all times with the result that the full benefits of MPC were realized. Fig. 10 shows the mismatch between plant measurements and RBF model predictions for the setpoint change in purge gas composition. As expected, process-model mismatch is essentially reduced to measurement noise. This result is typical for all three previously discussed cases. Disturbance rejection associated with IDV(1) (step change in A/C feed ratio in stream 4) and the corresponding process-model mismatch are shown in Figs. 11 and 12, respectively. The unmeasured step disturbance is first detected by a drop in a controlled vari-
able, the A/C molal ratio in the reactor feed (Fig. 13). Consequently, the A feed rate is increased to bring the A/C ratio in reactor feed back to setpoint. Since IDV(1) was a step disturbance, process-model mismatch ultimately reaches steady-state values (Fig. 12). Satisfactory control performance is realized since the NMPC algorithm employs plant-model mismatch as output feedback.
5.3. Discussion In addition to the ratio or speedup factor, the absolute values of the computation times required for the factorized and non-factorized approaches are impor-
1196
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
Fig. 9. Purge B composition setpoint change ( +2 mol%).
Fig. 10. Process-model mismatch during implementation of + 2% setpoint change in composition of component B in purge stream.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
1197
Fig. 11. Disturbance IDV(1) (step change in A/C feed ratio in A&C stream).
tant. Use of a non-factorized controller that requires 38 + h of computation time during 50 h of plant operation is unrealistic. The factorized algorithm completed the control calculations in 1.82 h or 3.6% of the total run time. In this case, the difference in computation costs has a major impact on the practicality of utilizing an RBF-based NMPC control algorithm. The performance results presented are expected based on the nonlinear modeling capabilities of RBF neural networks. Hence, an NMPC scheme that employs an RBF model would be expected to provide good control whenever the plant operated in regions used to train the model. Practical implementation of an RBF (or any empirical model) based MPC scheme
would require a watchdog to verify that operation lies within model development bounds. The problem/concept is the same as required for the linear MPC systems currently employed in industry. However, the level of concern is greater for an RBF based system due to the less certain (compared to a linear impulse or step response model) extrapolation characteristics. Another issue affecting practical implementation of an RBF based NMPC scheme is collection of the plant data necessary to develop the model. The amount of data required to develop any nonlinear model is orders of magnitude greater than required for a simple linear model. However, the time and investment required to develop the linear models used with today’s linear MPC
1198
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
Fig. 12. Process-model mismatch during occurrence of disturbance IDV(1). The unmeasured disturbance is not modeled by the RBF network, leading to poor predictions. However, due to the step nature of the disturbance, the process-model mismatches settle to near-steady values.
systems is already an impediment in many cases. One area of research that can potentially address this concern is model identification using closed-loop data. Falcioni and Seborg (1997) have studied this issue using neural network models. Bomberger et al. (1998) used historical data to identify RBF models for a copolymerization reactor. Potential improvement in performance using nonlinear control methods provides the incentive to find new ways to develop robust RBF models from existing closed-loop plant data.
6. Conclusions From a manufacturing standpoint, the development of faster and more powerful computers provides the opportunity to apply the concepts of MPC with process models of almost any type. Computationally intensive modeling methods such as RBF networks can now be realistically considered to achieve the benefits of nonlinear MPC. However, good engineering practice demands efficient utilization of all resources, including computation. Methods to improve the computational efficiency of new or improved control algorithms as addressed in
Fig. 13. The step disturbance IDV(1) is detected by the drop in the controlled variable, A/C ratio in reactor feed, from its setpoint. The variable is then brought back to its setpoint by increasing A feed rate.
S. Bhartiya, J.R. Whiteley / Computers and Chemical Engineering 26 (2002) 1185–1199
this paper are essential to accelerate potential acceptance and implementation in practice. Among the issues affecting potential application of RBF-based NMPC, model identification and maintenance remain the most prominent. As discussed previously, these issues represent hurdles for essentially all forms of NMPC. The highly distributed, non-compact form of RBF models is undesirable from a maintenance standpoint where simplicity and parsimony are highly valued. Ultimately, the nonlinear modeling capabilities of RBF networks provide the economic incentive to overcome this and all other impediments to industrial implementation. Future work on the factorized RBF-based NMPC technique will progress along two lines: identification and control. In the identification phase, development of RBF models using closed-loop data will be stressed. Also the broader question of input sequence design for nonlinear system identification will be explored. The control phase will expand the current version of NMPC with output feedback to optimal state estimation procedures.
Acknowledgements The authors express their gratitude to Professor N.L. Ricker and his group for use of their MATLAB simulation of the Eastman process.
References Banerjee, A., Arkun, Y., Ogunnaike, B., & Pearson, R. (1997). Estimation of nonlinear systems using linear multiple models. American Institute of Chemical Engineering Journal, 43, 1204 – 1226. Bhartiya, S., & Whiteley, J. R. (2001). Factorized approach to nonlinear MPC using a radial basis function model. American Institute of Chemical Engineering Journal, 47, 358–368. Bomberger, J. D., & Seborg, D. E. (1998). Determination of model order for NARX models directly from input –output data. Journal of Process Control, 8, 459–468. Bomberger, J. D., Seborg, D. E., & Ogunnaike, B. A. (1998). Radial basis function identification of a solution copolymerization reactor model (pp. 3182 – 3188). Philadelphia, PA: ACC. Chen, S., Billings, S. A., & Grant, P. M. (1992). Recursive hybrid algorithm for non-linear system identification using radial basis function networks. International Journal of Control, 55, 1051 – 1070.
1199
Cybenko, G. (1987). Approximation by superposition of a sigmoidal function. Mathematical Control Signal Systems, 2, 303. Downs, J. J., & Vogel, E. F. (1993). A plant-wide industrial process control problem. Computers and Chemical Engineering., 17, 245– 255. Falcioni, M., & Seborg, D. E. (1997). On the identification of neural network models during closed-loop operation. In European control conference. Brussels. Garcia, C. E. (1984). Quadratic/dynamic matrix control of nonlinear processes: An application to batch reaction process. In AIChE annual meeting. San Francisco, CA. Gattu, G., & Zafiriou, E. (1995). Observer based nonlinear quadratic dynamic matrix control for state space and input/output models. The Canadian Journal of Chemical Engineering, 73, 883 –895. Gurumoorthy, A., & Kosanovich, K. A. (1998). Improving the prediction capability of radial basis function networks. Industrial Engineering and Chemical Research, 37, 3956 – 3970. Hartman, E. J., Keeler, J. D., & Kowalski, J. M. (1990). Layered neural networks with Gaussian hidden units as universal approximations. Neural Computation, 2, 210 – 215. He, X., & Asada, H. (1993). A new method for identifying orders of input – output models for nonlinear dynamic systems (pp. 2520 – 2523). San Francisco, CA: ACC. Henrique, H., Lima, E. L., & Seborg, D. E. (2000). Model structure determination in neural network models. Chemical Engineering Science, 55, 5457 – 5469. Hush, D. R., & Horne, B. G. (1993). Progress in supervised neural networks; What’s new since Lippmann, IEEE Signal Processing Magazine, 10, 8 – 39. Hussain, M. A. (1999). Review of the applications of neural networks in chemical process control-simulation and online implementation. Artificial Intelligence in Engineering, 13, 55 – 68. Kennel, M. B., Brown, R., & Abarbanel, H. D. I. (1992). Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical Re6iew, 45, 3403 – 3411. Krishnan, A., & Kosanovich, K. (1998). Batch reactor control using a multiple model-based controller design. The Canadian Journal of Chemical Engineering, 76, 806 – 815. Lee, J. H., & Ricker, N. L. (1994). Extended Kalman filter based nonlinear model predictive control. Industrial Engineering Chemical Research, 33, 1530 – 1541. McAvoy, T. J., & Ye, N. (1994). Base control for the Tennessee Eastman problem. Computers and Chemical Engineering, 18, 383– 413. Moody, J., & Darken, C. J. (1989). Fast learning in networks of locally-tuned processing units. Neural Computation, 1, 281–294. Poggio, T., & Girosi, F. (1990). Networks for approximation and learning. Proceedings of the IEEE, 78, 1481 – 1497. Pottman, M., & Seborg, D. E. (1992). Identification of nonlinear processes using reciprocal muliquadric functions. Journal of Process Control, 2, 189 – 203. Rhinehart, R. R. (1992). A CUSUM type on-line filter. Process Control and Quality, 2, 169 – 176. Ricker, N. L. (1996). Decentralized control of the Tennessee Eastman challenge process. Journal of Processing Control, 6, 205 –221. Sanner, R. M., & Slotine, J. J. (1992). Gausian networks for direct adaptive control. IEEE Transactions on Neural Networks, 3, 837– 863.