2010
Ind. Eng. Chem. Res. 2000, 39, 2010-2023
Robust Adaptive Predictive Control of Nonlinear Processes Using Nonlinear Moving Average System Models Yugender Chikkula Advanced Control and Optimization Division, Aspen Technology Inc., 9896 Bissonnet, Houston, Texas 77036-8220
Jay H. Lee* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47906
An adaptive predictive control algorithm is presented for nonlinear moving average systems with parametric uncertainty. The algorithm is developed in the stochastic optimal control framework in which the parameters are modeled as random processes and their probability distributions are recursively updated and used explicitly in the optimal control computation. The framework yields an open-loop optimal feedback control (OLOFC) algorithm in which openloop optimal input trajectories minimizing the expectation of a multistep quadratic loss function are computed repeatedly as feedback updates occur. The algorithm is shown to be robust with respect to parametric uncertainty, with features such as on-line parameter refinement (and/or adaptation) and “cautious control”. Some potential additions to incorporate an active-learning feature to an otherwise passive-learning OLOFC are considered. Numerical examples are provided to illustrate the merits of the proposed method. 1. Introduction The dynamic behavior for most chemical processes is inherently nonlinear. Despite this fact, the vast majority of the existing industrial controllers are designed based on linear control techniques. These linear controllers are typically detuned and retuned in order to account for the underlying process nonlinearities, often resulting in conservative control laws that require high maintenance. For some processes, linear control may be intrinsically incapable of maintaining a stable operation. In these situations, a controller that accounts for some of the process nonlinearity can expand the region of operation and improve the overall closed-loop performance.1,2 Despite several important advances made in the nonlinear control system synthesis area, the availability of nonlinear dynamic models has been the main limiting factor for a widespread industrial practice of nonlinearmodel-based control. When a fundamental model is unavailable and prohibitively expensive to develop, empirical modeling is a practical alternative.3-5 Primarily inspired by linear time-series models, a number of nonlinear model structures have been suggested that relate past inputs/outputs to present/future outputs. They range from the relatively simple block-oriented structures to the more general Volterra or NARMAX structures. These models, typically identified using limited amounts of input/output data gathered from plant tests, have both structural and parametric uncertainties about them and should be used with caution for control purposes. Structural uncertainty depends on the choice of model structure and operating range within which the model is used. Model structure selection and uncertainty * To whom all correspondence should be addressed. Phone: (765)494-4088. Fax: (765)494-0805. E-mail:
[email protected].
characterization for general nonlinear systems are very difficult issues, and solutions do not appear to be within reach any time soon. Parametric uncertainty is also important because the commonly used model structures contain a large number of parameters and persistent excitation conditions (conditions for parameter convergence) for most of them are either not understood or difficult to meet during a typical plant experiment. For instance, the pseudo random binary signal input, which is popular for linear system identification, turns out to be a poor choice for identifying a second-order Volterra model.6,7 The focus of this work is placed on the issue of parametric uncertainty in the context of empirical nonlinear-model-based control. Parametric uncertainty can be combated in two ways. First, one can account for it explicitly in the controller design or in the control input optimization by minimizing their effect on the closed-loop performance. This is the premise of robust control. Second, one can reduce the level of uncertainty and distribute it in a manner least damaging to the closed-loop control, by intelligently designing on- and off-line identification signals and also by fully utilizing available information, including on-line data as in adaptive control. Of course, it would be best if both aspects can be incorporated into a single method. In this paper, we first present a framework for robust model predictive control (MPC) of nonlinear systems represented by nonlinear moving average (NMA) model structures. The Hammerstein and Volterra structures are special cases of the NMA structure, and these structures have been shown to well approximate many industrially important processes.1,6,8 For NMA systems we derive a model-based predictive control scheme in a stochastic optimal control setting. In this approach the parameters are modeled as random variables with a given prior density. Then their posterior density is computed recursively and used in the control computa-
10.1021/ie990393e CCC: $19.00 © 2000 American Chemical Society Published on Web 04/20/2000
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2011
tion. Following the idea of model predictive control, an open-loop optimal control problem is solved at each sample time and the solution is implemented in a receding horizon fashion, resulting in the so-called “open-loop optimal feedback control (OLOFC)” policy. A prototypical algorithm is developed which consists of two components: (1) a Kalman filter used for the calculation of the posterior density of the parameter vector and (2) an on-line optimizer similar to that appearing in the conventional MPC but with an additional nonlinear input penalty term accounting for the parametric uncertainty. The input weighting scheme is different from that for the conventional MPC algorithm in that it is directional and time-varying; i.e., various nonlinear modes are penalized in accordance with the level of their uncertainty, as shown by the parameter covariance matrix. The control algorithm is robust with respect to parametric uncertainty, with features such as on-line parameter refinement (or adaptation) and “cautious control”. One potential disadvantage of the control algorithm is the lack of active probing, which can make the controller turn off itself when the uncertainty becomes too large. We discuss two different ways to incorporate the active probing feature into the OLOFC. The paper is organized as follows. In section 2 we provide a brief overview of the general control problem for uncertain systems. In section 3 we develop the proposed nonlinear-model-based control strategies, with a detailed discussion of the underlying theoretical as well as practical issues. In section 4 we present numerical examples to illustrate various aspects of the proposed algorithms. Finally, in section 5 we give some concluding remarks. 2. Review of Robust Control of Systems with Parametric Uncertainty The design of feedback systems in the presence of parametric uncertainty has been widely investigated over the last few decades. Depending on the uncertainty description considered, the proposed methods can be classified into two categories: the hard-bound approach and the probabilistic approach. In the hard-bound approach, deterministic bounds are assigned to the parameters. Proposed control formulations for this uncertainty description include (1) design of a fixed parameter controller with constraints given by robust stability/performance conditions,9,10 (2) on-line optimization of which the objective is to minimize the “worstcase” error, i.e., the largest possible error over a given set of models.11-13 The drawbacks of the former type methods are that the optimization to be solved is usually nonconvex and the conditions for robust stability and performance are often conservative. In addition, these methods lack flexibility and are difficult to generalizes to cases involving input constraints, nonlinear models, etc. The latter type methods, while being more flexible and general, require solving a computationally demanding min-max optimization problem on-line. Moreover, the worst-case minimization results in a conservative control. In the probabilistic approach, the parameters are modeled as random variables and an optimal control problem is formulated, such as the minimization of the expected future output error. The dual control14 and various suboptimal control strategies such as the OLOFC15 belong to this category. This approach can
potentially lead to algorithms that are computationally less demanding and less conservative compared to the deterministic approach. With a few exceptions,16-18 most of the research effort in this area has been directed at linear systems. On the other hand, we believe that parametric uncertainty is more relevant and important for nonlinear systems for the reasons mentioned earlier. The recent work of Genceli and Nikolaou16 deals with parametric uncertainty of a second-order Volterra model in the nonlinear MPC setting. They provide sufficient conditions for robust stability of nonlinear MPC for a particular parametric uncertainty description. The uncertainty description they consider, independent hard parameter bounds (the probabilistic equivalent of which would be a mutually independent truncated uniform distribution for all of the parameters), however, can be conservative. In addition, the potential availability of on-line dynamic information which can be used to tighten parameter bounds is not considered. On-line parameter adaptation and evolution of the uncertainty in a deterministic setting can be demanding from a computational standpoint. The probabilistic framework, on the other hand, facilitates computationally simpler implementations (e.g., Kalman filter, extended Kalman filter (EKF), etc.). Although with the probabilistic approach it is difficult to guarantee stability, it is conceivable that some constraints (e.g., contraction constraint) can be added to the optimization to ensure closed-loop stability for all parameters within some probability level set. 3. Problem Formulation 3.1. System Model. To keep the presentation simple, all of the developments are carried out for single-input/ single-output (SISO) systems. We will consider systems or operating domains of systems that are well described by the following SISO linear-in-parameter NMA structure. M
yk )
Ri,kfi,k-1 + dk ∑ i)1
(1)
where
fi,k ) fi(uk,uk-1, ..., uk-N+1)
(2)
yk is the measured output and Ri,k, i ) 1, ..., M, are the system parameters. fi,k-1, i ) 1, ..., M, represent the nonlinear modes and are functions of {uk-1, ..., uk-N}, the past input sequence. The sequence dk represents the unstructured uncertainty which includes the effects of structural mismatch, output disturbances, effects of sampling a continuous-time system, etc. Note that the Hammerstein and Volterra structures are special cases of the above NMA structure. Rewrite the above model to obtain
yk ) φTk θk
(3)
θTk ) [R1,k R2,k ... RM,k dk]
(4)
φTk ) [f1,k-1 f2,k-1 ... fM,k-1 1]
(5)
where
The parameter vector φ, which also includes the dis-
2012
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000
turbance term θk, is assumed to be partially known with uncertainty because of a limited off-line identification experiment or the time-varying nature of the process. To account for these uncertainties, we model the parameter vector θ as a random variable with a given prior distribution characterizing the level of confidence. Further assuming the parameter vector to be a Gauss-Markov process, we rewrite the above system description in the following state-space form:
θk+1 ) Φθk + vk
p
uk
ωk ∼ N(0,Qw); vk ∼ N(0,Qν) This means the parameters and the residual term are either constant or change as Wiener processes. In the absence of better knowledge, this description represents a good choice. The choice for Qν further characterizes the time variation of the parameters. Note that, if one chooses all of the elements of Qν corresponding to the system parameters (R1, ..., RM) to be zero, the parameters are modeled as just constants. In this case, the estimation will refine the parameter estimates but will not be able to follow parameter changes. The initial parameter vector θ0 is also modeled as a Gaussian random variable vector independent from wk and vk:
θ0 ∼ N(θˆ 0,P0) The mean and covariance of θ0 are assumed to be available, for example, from an off-line identification experiment. 3.2. Parameter Estimation. On the basis of the above stochastic parameter model, the parameter refinement/adaptation and the computation of the associated uncertainty can be conveniently formulated within the Kalman filter framework as follows:
θk|k ) θk-1|k-1 + Lk(yˆ k - φTk θk-1|k-1) Lk ) Pk|k-1φk{φTk Pk|k-1φk + Qw}-1 Pk|k-1 ) Pk-1|k-1 + Q
(8)
Now let us consider the control problem to be the minimization of the following multistep objective function involving output deviations from a prespecified reference trajectory over a finite horizon into the future:
min E{
In the above equation, wk represents the measurement noise and yˆ k the noise-corrupted output. For simplicity of the presentation, however, we will assume that Φ is an identity matrix and both external signals (wk and vk) are independent, identically distributed Gaussian random vectors, i.e.,
Pk|k ) Pk|k-1 -
Ik ) {θk|k, Pk|k, Xk}
(6)
yˆ k ) φTk θk + wk
Pk|k-1φk{φTk Pk|k-1φk
probability density function of the Gaussian parameter vector and the past input sequence Xk ) {uk-N, ..., uk-1}:
(yk+l - rk+l)2|Ik} ∑ l)1
(9)
where E denotes the mathematical expectation with respect to all of the underlying random variables θk, {wi, vi|i ) k, ..., k + p}. Depending on the assumption about the nature of the future input decisions (uk+1, uk+2, ...), solutions to the above minimization problem differ drastically. (a) Closed-Loop Optimal Solution. If we assume that future inputs will be adjusted in a feedback-optimal sense, the future inputs must be treated as stochastic variables because future measurements are stochastic. Derivation of the optimal control policy in this case requires solving a stochastic dynamic program.19 The resulting solution is said to be “closed-loop optimal”. Even though this approach defines the problem we want to solve for optimal feedback control, the stochastic program is not solvable in most cases, including the specific case we are considering herein. (b) Open-Loop Optimal Solution. Alternatively, if we assume that the future inputs will be adjusted on the basis of Ik only, the input variables can be optimized as deterministic variables. We will show that we can evaluate the expectation analytically and reduce the problem to a deterministic optimization problem. While open-loop control is rarely used in practice, open-loop optimal solutions can still be used in feedback control by implementing them in a receding horizon fashion after each feedback update (as in MPC). 3.4. Closed-Loop Optimal Control. Define
Vk+l-1,k+p ) min E{(yk+l - rk+l)2 + Vk+l,k+p|Ik+l-1} uk+l-1
as the minimum expected loss for the part of the control horizon from k + l - 1 to k + p given the data up to k + l - 1. For l ) p we have
Vk+p-1,k+p ) min E{(yk+p - rk+p)2|Ik+p-1} uk+p-1
ν
+
Qw}-1φTk Pk|k-1 (7)
Here the pair (θk|k, Pk|k) defines the posterior density of θk and represents the best estimate and the associated uncertainty of the parameter vector at time k. 3.3. Control Problem Formulation. First let us define the information state Ik of the system, an important concept in stochastic optimal control theory. Ik is defined as the information about the system that is available to the controller at time k, after the measurement yk is taken but before the control moves {uk, uk+1, ...} are computed. For the above system Ik can be represented by the following finite dimensional vector containing the first two moments for the a posteriori
Note that this is a one-step-ahead minimization problem. Substituting from the model equations and evaluating the expectation and retaining terms involving uk+p-1, we have T θk+p|k+p-1 - rk+p)2 + Vk+p-1,k+p ) min {(φk+p uk+p-1
T Pk+p|k+p-1φk+p} φk+p
The objective function is a polynomial in uk+p-1 and the optimization can be solved by using any root solver. An important observation to be made is that the optimal solution and the optimal cost Vk+p-1,k+p are functions of the parameter uncertainty characterized by Pk+p|k+p-1 as well as the parameter estimate θk+p|k+p-1 and the state χk+p-1.
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2013
Going to the next stage, for l ) p - 1, dynamic programming yields the following single-stage minimization:
In the above equations, T yk+l|k ) φk+l θk|k
Pk+l|k ) Pk|k + lQν
Vk+p-2,k+p ) min E{(yk+p-1 - rk+p-1)2 + uk+p-2
Vk+p-1,k+p|Ik+p-2} Note that the optimal cost for the previous stage (which is referred to as the “cost-to-go”) appears in the objective function. In principle, by continuing the same procedure recursively, we arrive at the following minimization at time instant k (from which the optimal value of uk can be found):
Vk,k+p ) min E{(yk+1 - rk+1) + Vk+1,k+p|Ik} 2
uk
Unfortunately, the above minimization does not yield itself to an analytical solution. The main difficulty is that the dependence of the cost-to-go function on the information vector cannot be expressed as a closed-form function. Without it, the dynamic program for a subsequent stage cannot be formulated. Numerical solution is not a practical option because of the computational complexity and the “curse of dimensionality” involved in evaluating the cost-to-go function. Despite this, the solution procedure provides an important insight into the control problem; i.e., the optimal controller is a dual controller. In other words, the optimal controller not only would try to regulate the system close to the prespecified setpoint trajectory but also would spend some energy to elicit dynamic information from the system and actively reduce future uncertainty. This can be seen from the second stage optimization for l ) p - 1. The cost-to-go function for the previous stage optimization (Vk+p-1,k+p) that appears as part of the objective function depends on Pk+p|k+p-1, which belongs to Ik|k+p-1. The “smaller” Pk+p|k+p-1, the smaller Vk+p-1,k+p. Note that Pk+p|k+p-1 is effected by uk+p-2, which appears as part of the output matrix φk+p-1. This connectivity induces active probing that is optimally balanced with the regulation objective. 3.5. Open-Loop Optimal Feedback Control. Now let us consider the scenario where we make the decisions for all of the future inputs on the basis of Ik only. Under this scenario, the future inputs are treated as deterministic variables and the optimal control problem can be stated as p
min uk,...,uk+p-1
E{
(yk+l - rk+l)2|Ik} ∑ l)1
(10)
This implies that we are considering an open-loop control starting at time k. Feedback control results if we implement only the first control input uk from the computed control sequence and repeat the optimization at subsequent time steps (after shifting the window by one time step at each time). This defines a suboptimal feedback control strategy referred to as OLOFC. We expand each term in the open-loop objective given above as follows:
E{(yk+l - rk+l)2|Ik} ) E{[(yk+l|k - rk+l) + (yk+l - yk+l|k)]2|Ik} T Pk+l|kφk+l ) [yk+l|k - rk+l]2 + φk+l
(11)
(12)
Thus, the open-loop optimal control problem to be solved at each time step is the following constrained nonlinear program (NLP): p
min
T Pk+l|kφk+l] ∑[(yk+l|k - rk+l)2 + φk+l
(13)
uk,...,uk+p-1 l)1
This has the look of classical MPC in that the objective function consists of quadratic terms for the output errors and input terms (appearing in φk+l). To control the size of the optimization, the number of inputs considered can be set as m e p, with the remaining p - m input values constrained as uk+m-1 ) uk+m ) ... ) uk+p-1. Additional constraints on uk, ..., uk+p-1 and yk+1, ..., yk+p can be included in the optimization:
umin e uk+l-1 e umax ymin e yk+l|k e ymax, 1 e l e p
(14)
Comments. (a) In the completely deterministic case, both the closed-loop and open-loop formulations result in the same solution. However, in the case of a stochastic parameter vector, the two approaches produce very different results. The open-loop approach results in a tractable problem but yields only a suboptimal feedback strategy. (b) The OLOFC involves minimization of a multistep quadratic objective, similar to that found in the classical MPC formulation. The difference, however, is that various linear combinations of nonlinear input terms contained in φk+l are now penalized in proportion to the corresponding parametric uncertainty as expressed by matrix Pk+l|k. Changes in the model quality are reflected as changes in the covariance matrix Pk|k, which is recomputed at each time step using the Kalman filter. In the context of conventional MPC, it is well-known that the input weighting term has a significant effect on the closed-loop speed and robustness and a constant diagonal input weighting matrix is not always the best choice. However, there has not been any systematic way to choose a nondiagonal input weighting matrix that varies according to changing uncertainty. The proposed formulation results in an automatic, time-varying, and directional tuning of input weights, based on the conditional covariance of the parameter vector computed under the assumed stochastic parameter model. It is noteworthy that the input penalty term penalizes combinations of both past and future input terms. (c) In the adaptive control context, the OLOFC formulation represents a way to relax the “certainty equivalence” (CE) assumption. In CE control, the estimated parameters are treated as though they were exact; i.e., uncertainties in parameter estimates are not considered. This assumption is reasonable in cases where the number of parameters is small and the parameter variation is slow compared to that of the system state. However, it is not justifiable when parameter variances are large (e.g., during initial phases
2014
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000
of adaptation) or when rapid and sudden process changes occur requiring reidentification. In both cases the uncertainties would be significant because of the large number of parameters in typical NMA models. Overall, the result is a cautious controller with a “multistep” prediction horizon, which can be viewed as an extension of the classical “one-step” cautious control.19 (d) Some readers may question the reason for restricting our formulation to NMA structures. With more general structures involving autoregressive terms, the future outputs are nonlinear in parameters and explicit evaluation of the expectation operator becomes intractable (except for the case when the horizon is 1). Numerical evaluation (e.g., based on a quadrature approximation of the integral) is an option but would be computationally expensive. (e) Even though the stochastic problem has been converted into an equivalent deterministic problem, the optimization one must solve is still a NLP. However, this is a problem common to all MPC techniques involving nonlinear system models, including those that do not consider model uncertainty. The complexity of the NLP would depend on the choice of nonlinear basis functions, and solution algorithms tailored to a specific choice can potentially be developed. However, this is beyond the scope of the present paper. 3.6. Active Probing and Suboptimal Dual Control. One well-known attractive feature of the optimal feedback solution for systems with uncertain parameters is the active probing which acts to reduce future parameter uncertainty in a control-friendly manner. However, the optimal dual control problem is intractable even for the simplest of practical problems. The OLOFC, though cautious against parametric uncertainty, does not result in an active probing because the open-loop control assumption made in the optimization does not reward generation of useful feedback information that can reduce the uncertainty at future time steps. In fact, the mere “cautiousness” of the OLOFC can cause it to turn off when the covariance matrix becomes very large. Practically, this means that off-line identification should be performed first to reduce the uncertainty to a reasonable level before the OLOFC is activated. To remove this limitation, there are several ways to endow the passive-learning OLOFC with the activelearning feature. The simplest would be to superimpose an external perturbation signal onto the OLOFC solution. However, for nonlinear model structures, a priori specification of an appropriate input perturbation sequence that would excite all of the nonlinear modes may be difficult.7 (See refs 6 and 7 for some results on secondorder Volterra systems.) Moreover, this approach does not consider the fact that information may be lacking (or needed) in only a few parameter directions, and any additional probing would be useful only if it is in the uncertain parameter directions. One could alternatively design the external probing signal on-line based on the updated covariance matrix and some optimization criterion based on a measure of the covariance matrix. Let us consider some of the options and their relative merits and demerits. 3.6.1. Including Excitation Constraints in the OLOFC. For a model structure with a regressor vector φ, robust convergence of parameter values requires the sequence of regressor vectors φk to persistently span the parameter space. This spanning condition can be satis-
fied if there exists an integer T such that k+T
F1 I e
φjφTj e F2I ∑ j)k
∀k
(15)
where F1, F2 > 0. For all systems exhibiting an affine dependence on the parameter vector θ, the lower bound of the above persistent excitation (p.e.) condition is a fundamental requirement for robust identifiability of the parameters. For the linear time invariant (LTI) case, the p.e. conditions are well understood and can be specified in terms of the frequency content of the input signal. However, for nonlinear systems, the system dynamics can potentially play tricks on the input frequencies because of the capability of certain nonlinear systems for subharmonic generation. For instance, N/2 distinct frequency components, which suffice for a n parameter LTI case, may not be sufficient for an n parameter nonlinear system. The converse may also be true; i.e., with n/2 distinct input frequencies we may be able to identify a nonlinear system with more than n number of parameters, which will not be possible in a LTI case. Note that the positive semidefiniteness condition in (eq 15) can be reformulated in terms of the eigenvalues k+T φjφTj - F1I and of the information matrix. Because ∑j)k k+T T F2I - ∑j)k φjφj are real, symmetric matrixes, all of the corresponding eigenvalues are real. The constraints k+T φjφTj can be required to satisfy the bound on F2I - ∑j)k eliminated provided one has constraints on the input, which generally is the case in practice. Bounded inputs automatically guarantee the upper bound in eq 15. Now the active learning can be induced into OLOFC by solving the OLOFC problem with a set of additional constraints on the minimum eigenvalue of the information matrix as follows: p
min uk,...,uk+p-1
T [(yk+l|k - rk+l)2 + φk+l Pk+l|kφk+l] ∑ l)1
(16)
subject to
umax g uk+l-1 g umin ymax g yk+l|k g ymin ∆umax g |uj - uj-1|, j ) k, ..., k + m - 1 k
λi(
φjφTj - FI) g 0, ∑ j)k-T
i ) 1, ..., n
· · · · · · · · · k+m-1
λ i(
φjφTj - FI) g 0, i ) 1, ..., n ∑ j)k+m-1-T
All sequences of feasible solutions for the above problem guarantee persistency of excitation in the parameter space. In practice, however, one needs to specify some of the parameters in the above algorithm such as the excitation horizon T, the bound on the minimum singular value of the information matrix F. The following guidelines may be useful for choosing these parameters.
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2015
(a) Choose T g n where n is the dimension of the parameter space. (b) Choice of F reflects the minimum amount of information in the least informative direction. The parameter convergence will depend on the choice of F. Theoretically any F > 0 suffices for guaranteeing p.e. However, excessively small values may result in slow convergence. At the same time large values may result in unnecessary probing leading to deterioration in closed-loop performance. In addition, the constraints can become infeasible. The choice will also depend on the level of measurement noise variance. This idea has been explored by Vuthandam and Nikolaou20 in the context of linear MPC where simultaneous identification and control are achieved by including a weak p.e. condition as a constraint in the on-line optimization problem. 3.6.2. Optimal Design of the Probing Signal. In certain cases, rather than constraining the eigenvalues of the information matrix, it may be more practical to design a magnitude-bounded probing signal directly. For example, as shown in Figure 12, one can consider the probing signal generation in conjunction with the closedloop control problem and solve for it separately. This can be formulated as a two-stage optimization problem. In the first stage, the OLOFC problem (13) is solved. Let us denote the OLOFC solution as follows: {ukc, ..., uk+m-1c}. In the second stage, a probing signal within a certain magnitude limit is generated, which when implemented along with the OLOFC solution, {ukc, ..., uk+m-1c}, satisfies the overall constraints on the inputs. The probing signal is generated as the solution of the following optimization problem involving a measure on the covariance of the parameter vector as follows:
min
ukp,...,uk+m-1p
Ψ(Pk+m-1|k)
(17)
subject to the constraints on inputs and outputs
umax g uk+ic + uk+ip g umin, i ) 0, ..., m - 1 ymax g yk+i|k g ymin, i ) 1, ..., p ∆umax g |(uk+ic + uk+ip) - (uk+i-1c + uk+i-1p)|, i ) 0, ..., m - 1 p
p
∆umax g |uk+i |, i ) 0, ..., m - 1
(18)
Here Pk+m-1|k is the error covariance matrix obtained by propagating Pk|k through the Kalman filter Riccati equation with the input sequence {(ukc + ukp), ..., (uk+m-1c + uk+m-1p)}. Ψ is a measure of the covariance matrix such as determinant, largest eigenvalue, or trace. Note that an upper bound on the magnitude of the probing signal umaxp is also specified. umaxp can be chosen to prevent excessive probing. The solution is implemented in a receding horizon fashion, i.e., only ukc + ukp is implemented, and at the next sampling time, both eqs 13 and 17 are resolved. The constraints included in eq 18 guarantee that the plant operates within the specified constraints and the probing is carried out within the allowable plant perturbations. 3.6.3. Additional Perspective on the Active Probing Problem. An important shortcoming of the above two approaches to simultaneous identification and control is that the interaction and balance between the two are not considered in a rigorous and systematic way.
Figure 1. Actual (s) vs identified (+) parameters for the secondorder Volterra system in example 1.
A viable alternative is to optimize the current input under the assumption that some fixed control law (e.g., the CE control) would be used at future time steps.18 Hence, for the optimization of the current input, future inputs are assumed to be generated in a feedback fashion, just not in an optimal manner. This assumption greatly simplifies the original optimal dual control problem, but one, nevertheless, has to evaluate the expectation of future errors numerically, for example, by Monte Carlo simulation. Another way to balance the control and identification efforts is to incorporate control relevancy into the objective function used for the probing signal generation (Ψ(Pk+m-1|k). Refer to ref 21 for some ideas and inherent difficulties related to control-relevant input signal generation for NMA structures. 4. Numerical Examples Now let us consider a few numerical examples to illustrate various aspects of the proposed algorithms. 4.1. Example 1: A Second-Order Volterra System. Let us start with the following numerical example to illustrate the cautious characteristics of the proposed method. Consider a SISO second-order Volterra system and a second-order Volterra model identified from an input/output experiment. The input employed in the identification was a white noise of very small variance which essentially excites the linear dynamics. As can be seen from Figure 1, the estimates for the first-order terms (the first 10 terms) are close to the actual parameters, but the estimates for the second-order parameters are substantially off for obvious reasons. This model was implemented in a closed loop, using both OLOFC and the CE control strategy. The CE control strategy can be interpreted as a special case of OLOFC with Pk+l|k set to a zero matrix in eq 13. To make the comparison more transparent, we consider a simple scenario where we set Qw ) 0 (no measurement noise) and Qν ) 0 (constant parameter). The Kalman filter was initialized with the off-line leastsquares estimates of the parameters and the corre-
2016
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000
Figure 2. Comparison of OLOFC (solid line), modified OLOFC (dashed line), and the CE controller (dot-dashed line) for a secondorder Volterra system in example 1.
sponding covariance matrix. The following controller settings were used.
p ) 17; m ) 10 The CE controller and OLOFC are simulated for a step change in the setpoint, and the results are plotted in Figure 2. As can be observed, OLOFC is “cautious” compared to the CE controller, especially as the system moves away from the initial region of identification. The detuning of a CE controller is typically done by adding a quadratic penalty term on the manipulated inputs. However, the detuning effects are not transparent, and it is not straightforward to detune for the changing uncertainty. In OLOFC, on the other hand, the built-in automatic detuning mechanism takes care of this and results in a cautious controller. OLOFC can be overly cautious depending on the level of uncertainties. In these situations, the objective (13) can be modified as below to obtain a relatively more aggressive controller depending on the choice for λ. p
min uk,...,uk+m-1
T [(yk+l|k - rk+l)2 + λφk+l Pk+l|kφk+l] ∑ l)1
The above objective interpolates between CE and OLOFC for 0 e λ e 1. For λ ) 0 we have the CE control and for λ ) 1 the OLOFC. By choosing an appropriate value for λ, one can reduce the conservativeness of OLOFC or improve the level of robustness of CE controllers. λ should be viewed as an “on-line” tuning parameter. λ ) 1 would be a good starting value because it represents the most conservative (but safe) setting. Once the controller is on-line, the user can adjust it depending on the closed-loop responses. The above example for λ ) 0.1 is solved, and the closed-loop response is shown in Figure 2. Notice that the level of performance of this modified OLOFC is somewhere between those of the OLOFC and CE controllers. 4.2. Example 2: van de Vusse Reactor. Now let us consider the van de Vusse reaction system1,2 and demonstrate the cautious and directional detuning
Figure 3. van de Vusse reactor: steady-state characteristics. Table 1. Nominal Operating Conditions for the van de Vusse Reactor in Example 2 variable
value
variable
value
F0 V k1
25 L/h 1L 50 h-1
k2 k3 CAf
100 h-1 10 L mol-1 h-1 10 mol L-1
characteristics of the OLOFC. We will be comparing the algorithm with CE control. The following reaction is carried out in a CSTR:
AfBfC 2A f D The process dynamics can be described by the following material balances for components A and B:
dCA F ) -k1CA - k3CA2 + (CAf - CA) dt V dCB F ) k1CA - k2CB - CB dt V The design and kinetic parameters are given in Table 1. The control problem is to regulate CB, the concentration of B, using the inlet flow rate F as the manipulated input. F0 represents the steady-state flow rate. The corresponding output value, CB at steady state, is 1 mol/ L. The steady-state behavior of the reactor is shown in Figure 3 over the entire range of the input variable. As can be seen, the system is significantly nonlinear. There is a change in the sign of steady-state gain. To the left of the peak conversion, the system also exhibits inverse response (unstable zero dynamics). Note that the operating point considered here falls in the region corresponding to the nonminimum phase behavior of the process. A second-order Volterra model of the following form was identified from an input/output experiment around the nominal operating point. N
y′k )
N
N
aiu′k-i + ∑ ∑ bi,ju′k-iu′k-j ∑ i)1 i)1 j)i
(19)
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2017
Figure 4. Model 1 for the van de Vusse reactor: steady-state comparison (solid, system; dashed, model).
where y′k is the output deviation from its steady-state value of 1 mol/L and similarly u′k is the input deviation from its steady-state value of 25 L/h. A sampling time of 0.002 h was chosen. The choice for the sampling time was primarily motivated by the
inverse response characteristics of the system at the nominal operating conditions. We considered two scenarios for input/ouput experimentation. These two cases differ in the level of measurement noise at the process output. In the first case we considered very small measurement noise, Gaussian white noise with variance σ2 ) 10-10. A sequence of uniformly distributed input perturbations between the range 20 and 30 L/h were considered, and a total of 1000 samples were collected. The purpose of this ideal experimentation is to determine an appropriate model order N in eq 19. With N ) 20 (in eq 19), we performed a least-squares fit of the input/output data. A total of 230 model parameters is estimated, 20 first-order parameters ai and 210 second-order parameters bi,j. Note that an upper diagonal form of the second-order kernel is estimated here. It was observed that the first-order parameters (the impulse response) are close to zero at a delay of 20 sampling times. The choice of N ) 20 for the model order seems to be adequate for the range of input variations considered in the experiment. This is confirmed by comparing the model and process steadystate gains over the input range of interest and testing the model responses for low-frequency-type and highfrequency-type input variations. The process and model responses are compared in Figures 4-6. As is demonstrated, the identified second-order Volterra model and process responses are very close. Because the objective of this simulation study is to investigate the cautious and directional detuning capa-
Figure 5. Model 1 for the van de Vusse reactor: low-frequency validation (solid, system; dashed, model).
2018
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 Table 2. Controller Parameters for Example 2 parameter p m λy ymax ymin umax
Figure 6. Model 1 for the van de Vusse reactor: high-frequency validation (solid, system; dashed, model).
Figure 7. Model 2 for the van de Vusse reactor: steady-state comparison (solid, system; dashed, model).
bilities of the OLOFC algorithm, we next develop a model under a more realistic experimental condition. In this experiment we consider a higher level of measurement noise, Gaussian white noise of variance σ2 ) 0.0025. Input perturbations similar to those in the previous experiment are considered. From the collected data, we identified another secondorder Volterra model with the structure given in eq 19 (with N ) 20). By comparing these estimates with the estimates from the previous experiment, we observed that the measurement noise significantly affected the parameter estimates, especially the second-order parameters, bi,j (Figure 7). Figures 8 and 9 give the model responses for low- and high-frequency input perturbations. As can be observed, the model results in comparatively less accurate responses than the ones resulting from the previous model obtained from an idealized
setting 10 5 1 1.25 0 250
parameter umin ∆umax Qw Qv θ0|0 P0|0
setting 0 5 0.000001 0 θ0 P0
experiment. On the basis of the identified model in the second experiment, we now explore the cautious and directional detuning capability of OLOFC. An OLOFC controller was designed with the controller and estimator parameters given in Table 2. Here θ0|0 and P0|0 are the initial settings for the parameter vector and the associated covariance matrix in the Kalman filter, respectively. θ0 and P0 were chosen as the offline estimates of the parameter and the covariance matrix from the least-squares identification. Now the responses of the closed-loop system for perturbations in x1 are evaluated. The closed-loop system is also simulated with a CE controller with parameter settings identical with those for the OLOFC controller. Constr.m routine from MATLAB was used to solve the nonlinear optimization problems in all of the closed-loop simulations considered here. The performance of the OLOFC controller is compared with the CE controller in Figure 10 for a state perturbation in x1 of 0.5. As can be expected, OLOFC results in a much better response than the CE controller. OLOFC is cautious and robust with respect to the underlying uncertainty compared to the CE controller, and this characteristic can be established by observing the relatively large input variations in the CE controller. However, note that the CE controller can be modified with an additional term in the objective function, a term penalizing weighted inputs as in the classical MPC formulation. We explored this in the next set of simulations. Here, the performances of OLOFC and the CE controller with different amounts of input weighting are compared in Figure 11 for the same perturbation in x1. Again, the OLOFC response is superior to the modified CE controller response for the tuning scenarios considered. As can be observed from the response characteristics, the closed loop can be detuned considerably by choosing a large value for the input weight, λu ) 0.01. This gives a relatively low overshoot; however, the response is very sluggish. The response can be made faster by decreasing the input weight; however, for the small value of λu ) 0.000 01, the closed-loop system is highly oscillatory and eventually unstable. For an intermediate value of λu ) 0.001, the system responds faster but with a relatively large overshoot. Compared to the conventional detuning through the input weight, the inherent detuning in OLOFC is time-varying and directional. This can be observed from the response characteristics of OLOFC in Figure 11. The immediate response of the controller is as fast as that of the CE controller with little or no input weighting. This is because the high-frequency response of the model matches well with the system. In other words, the uncertainty in the first few parameters of the model is relatively small. As time progresses, the larger uncertainty in the other parameters assumes importance. However, the time-varying, nonlinear, and directional detuning inherent to the OLOFC results in a more
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2019
Figure 8. Model 2 for the van de Vusse reactor: low-frequency validation (solid, system; dashed, model).
Figure 9. Model 2 for the van de Vusse reactor: high-frequency validation (solid, system; dashed, model).
Figure 10. OLOFC (solid) vs CE control (dashed) of the van de Vusse reactor.
conservative response compared to the CE control after the initial faster response. This example demonstrates the directional detuning and cautious characteristics of the OLOFC framework. The integral square errors (ISE) over a time period of 100 sampling instants for all of the controllers are compared in Table 3. As can be seen, OLOFC results in the lowest ISE.
Even though not shown here, we have carried out several other simulations under various perturbations (e.g., a perturbation of -0.5 in x1), and the results obtained were similar. 4.3. Example 3: OLOFC with a Bounded Probing Signal. Now let us consider an example to illustrate the benefit of the on-line probing signal generation scheme presented earlier. We will consider a second-
2020
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000
Figure 11. OLOFC (solid) vs CE control with input weighting (dashed) of the van de Vusse reactor.
Figure 13. Closed-loop performance in example 3: setpoint (solid), OLOFC (+), and OLOFC with probing (dashed).
Figure 12. On-line generation of the probing signal. Table 3. Controller Performance for a State Disturbance in x1 of 0.5 in Example 2 controller CE controller NLMPC with λu ) 0.00001 NLMPC with λu ) 0.0001 NLMPC with λu ) 0.005 NLMPC with λu ) 0.008 NLMPC with λu ) 0.01 OLOFC
integral square error unstable 1.0947 0.2255 0.0614 0.0516 0.0493 0.0360
Figure 14. Parameter estimation in example 3: system (solid), OLOFC (+), and OLOFC with probing (dashed). Table 4. Controller and Estimator Settings for Example 3 parameter
order Volterra system yk ) φTk θp whose parameters are given as
θp ) [0.6000 0.2140 0.200 0.1500 0.10] φTk ) [uk-1 uk-2 uk-12 uk-1uk-2 uk-22] (20) The following parameter vector is assumed to have been identified:
θm ) [0.6441 0.2176 0.1936 0.1545 0.0950] As can be observed, we have assumed very good parameter estimates. The estimator was initiated with
p m Qν Qw ymax ymin
value
parameter
value
10 5 0 10-3 2.5 0
umax umin ∆umax umaxp λy
2.5 0 1 1 10
a covariance of P0 ) 0.001I. Now we will consider a scenario where the system parameters jump suddenly to the following parameter values:
θp,new ) [0.1000 0.0140 0.400 0.4500 -0.40] (21) In practice, this could be due to the system moving to a different operating point. This change is considered to
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2021
Figure 15. Closed-loop performance in example 4: OLOFC (dashed) and OLOFC with probing (solid).
occur after five sampling units from the beginning of the operation. A detection scheme is assumed to be in place to sense this change. At this time, the covariance in the Kalman filter is reset to P ) 0.9I. The probing signal is computed according to eq 17 with tr(Pk+m-1|k) as the minimization criterion. The choice of the criterion tr(Pk+m-1|k) was primarily motivated by numerical complexity. The controller and estimator parameters are given in Table 4. The probing is carried out for five sampling units starting 10 min after the jump is detected and is switched-off afterward. In Figures 13 and 14, we show the performance of the OLOFC controller with and without the probing signal generator. As can be seen in Figure 14, the parameter estimation is significantly better when a probing signal is explicitly computed and the proposed scheme implemented. 4.4. Example 4: OLOFC with Excitation Constraints. Now let us consider an example to illustrate the effect of adding the excitation constraints on the closed-loop estimation and control. We will consider the following second-order Volterra system:
yk ) φTk θp where
θp ) [0.26000 0.2140 0.1200 0.1500 -0.10 -0.05 -0.01 -0.001 0.02] φTk ) [uk-1 uk-2uk-3 uk-12 uk-1uk-2 ... uk-32] The following parameter is assumed to have been
identified from the above system:
yk ) φTk θm where
θm ) [0.1864 0.2603 0.1460 0.1490 -0.04 -0.07 -0.005 -0.1 0.085] As can be observed, there is significant mismatch between the model and the plant parameters. The estimator was initiated with a covariance of P0 ) 0.05I. The controller and estimator parameters are given in Table 5. Now let us consider the scenario where the probing characteristics of the OLOFC algorithm with excitation constraints can be demonstrated. For the first 27 sampling instants, OLOFC is used to control the system output at the steady-state value of 0.603. As can be observed from the estimator performance in Figure 16, at the beginning, the parameter estimates start drifting because of the lack of excitation. Then the excitation conditions are activated after 27 sampling instants from the beginning of the simulation. The performances of OLOFC and OLOFC with the excitation constraints are compared in Figures 15 and 16. As can be observed in Figure 16, the estimator performance is significantly better when the excitation constraints are added. The effect of improved parameter estimates on closed-loop performance is illustrated in Figure 15 for a step change in the setpoint from 0.603 to 1. The performance of OLOFC with excitation constraints is comparatively superior than when OLOFC
2022
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000
Figure 16. Parameter estimation in example 4: OLOFC (dashed) and OLOFC with probing (solid). Table 5. Controller and Estimator Settings for Example 4 parameter p m T Qν Qw
value 5 5 27 0 10-3
parameter ymax ymin umax umin F
value 10 -10 10 0 0.00000021
Acknowledgment The authors gratefully acknowledge the support of DuPont, Aspen Technology, and the National Science Foundation (Grant CTS-9357827). Literature Cited
is implemented without the excitation constraints. The overshoot is smaller, and the system settles faster.
(1) Doyle, F. J., III; Ogunnaike, B. A.; Pearson, R. K. Nonlinear model based control using second-order Volterra models. Automatica 1995, 31, 697-714.
5. Conclusions
(2) Hernandez, E. Control of nonlinear systems using inputoutput information. Ph.D. Dissertation, Georgia Institute of Technology, Atlanta, GA, 1992.
In this paper, we presented a novel nonlinear MPC algorithm that minimizes the expectation of a multistep quadratic loss function for NMA models. The algorithm was engineered to be cautious against parametric uncertainties by treating parameters as random variables with the mean and covariance given by the recursive least-squares identification. Using a series of examples, we demonstrated the claimed “cautious” characteristics of the proposed algorithm toward model uncertainty for the prototypical case of a second-order Volterra structure. In essence, the proposed scheme can be interpreted as a systematic way of relaxing the model inversion through the input weighting by explicitly including information about parameter uncertainty. We also presented additions to the OLOFC formalism that can potentially incorporate an active-learning feature into the otherwise passive-learning algorithm.
(3) Leontaritis, I. J.; Billings, S. A. Input-Output parametric models for nonlinear systems: parts I & II. Intl. J. Control 1985, 41, 303-344. (4) Diaz, H.; Desrochers, A. A. Modeling of nonlinear discretetime systems from input-output data. Automatica 1988, 24, 629641. (5) Pearson, R. K. Nonlinear input-output modeling. Preprints of ADCHEM’94, Kyoto, Japan, 1994. (6) Pearson, R. K.; Ogunnaike, B. A. Nonlinear Process Identification. In Nonlinear Process Control; Henson, M. A., Seborg, D. E., Eds.; Prentice-Hall: Upper Saddle River, NJ, 1993. (7) McCullough, G.; Maner, B. R.; Doyle, F. J., III; Pearson, R. K. Identification of second-order Volterra models from plant data for use in model predictive control. AIChE Annual Meeting, San Francisco, CA, 1994. (8) Eskinat, E.; Luyben, W. L. Use of Hammerstein models in identification of nonlinear systems. AIChE J. 1991, 37, 255-268.
Ind. Eng. Chem. Res., Vol. 39, No. 6, 2000 2023 (9) Barmish, B. R.; Hollot, C. V.; Kraus, F. J.; Tempo, R. Extreme point results for robust stabilization of interval plants with first-order compensators. IEEE Trans. Autom. Control 1992, AC-37, 707-714. (10) Bernstein, D. S.; Haddad, W. M. Robust controller synthesis using Kharitonov’s theorem. IEEE Trans. Autom. Control 1992, AC-37, 129-132. (11) Campo, P. J.; Morari, M. Robust model predictive control. Proceedings of the Automatic Control Conference, Atlanta, GA, 1987; pp 1021-1026. (12) Veres, S. M.; Norton, J. P. Predictive self-tuning control by parameter bounding and worst-case design. Automatica 1993, 29, 911-928. (13) Lee, J. H.; Yu, Z. Worst-case formulation of model predictive control for systems with bounded parameters. Automatica 1997, 33, 763-781. (14) Feldbaum, A. A. Optimal control theory; Academic Press: New York, 1965. (15) Tse, E.; Athans, M. Adaptive stochastic control for a class of linear systems. IEEE Trans. Autom. Control 1972, AC-17, 3851. (16) Genceli, H.; Nikolaou, M. Design of robust constrained nonlinear model predictive controllers with Volterra series. AIChE J. 1995, 41 (9), 2098-2107.
(17) Michalska, H.; Mayne, D. Q. Robust receding horizon control of constrained nonlinear systems. IEEE Trans. Autom. Control 1993, AC-35, 1623-1633. (18) Tse, E.; Bar-Shalom, Y.; Meier, L., III. Wide-sense adaptive dual control of nonlinear stochastic systems. IEEE Trans. Autom. Control 1973, AC-18, 98-108. (19) Åstro¨m, K. J.; Wittenmark, B. Adaptive control; AddisonWesley: Reading, MA, 1989. (20) Vuthandam, P.; Nikolaou, M. Constrained model predictive control and identification: A weak persistent excitation approach. AIChE J. 1997, 43 (9), 2279-2288. (21) Chikkula, Y.; Lee, J. H. Input sequence design for parametric identification of nonlinear systems. Proceedings of the American Control Conference, Albuquerque, NM, 1997; pp 30373041.
Received for review June 4, 1999 Revised manuscript received February 21, 2000 Accepted February 23, 2000 IE990393E