Ind. Eng. Chem. Res. 1990,29, 1647-1657
1647
PROCESS ENGINEERING AND DESIGN Newton-Type Controllers for Constrained Nonlinear Processes with Uncertainty Wei-Chong Li and Lorenz T. Biegler* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
Recently developed Newton-type controllers have been demonstrated for nonlinear processes in both single-step and multistep applications. Derived from concepts related to nonlinear operator theory, these model-based controllers can be extended to deal with both input and output constraints. Also, concepts derived from operator theory, such as descent properties and contraction mappings, are invoked to prove stability properties for these methods. However, Newton-type controller properties were derived under the assumption of a perfect process model. In this paper, we relax this assumption and demonstrate under which conditions stability properties still hold for these systems. In particular, we deal with the cases of additive output uncertainty (such as measurement noise) and general model mismatch. T o deal with the latter case, a parameter estimation algorithm is applied to reduce the plant/model mismatch. The combined Newton-type controller with parameter estimation is applied to a pH control process and is compared to literature results for Nonlinear Inferential Control. In addition, cases with significant measurement noise are also examined. Results on this example problem indicate the effectiveness of this control strategy for nonlinear control. 1. Introduction
Newton-type control strategies for nonlinear processes were first proposed by Economou (1985) and were demonstrated on a small nonlinear reactor problem for which linear control laws became unstable. Li and Biegler (1988) and Li et al. (1990) extended these concepts by adding a relaxation parameter in order to enforce global convergence. In addition, this approach was also adapted through quadratic programming to include constraints on process inputs and (linearized) outputs. More recently, Li and Biegler (1989) extended these Newton-type methods to a multistep formulation that could be implemented within a moving time horizon framework. However, the algorithms developed in these studies are based on a very restrictive assumption, i.e., that the model is perfect. Unfortunately, it is rare that the perfect model can be derived or obtained through experiments without any uncertainties. Morari and Zafhiou (1989) provide an excellent approach to the design of robust linear controllers where the uncertain plant can be represented as a family of linear plants in the frequency domain. In addition, they provide necessary and sufficient conditions for robust performance and robust stability for linear processes with this uncertainty description. However, when nonlinear systems are involved, analysis for robustness becomes more difficult. A major reason for this is that, unless a transformation to a linear structure is apparent, analysis of the nonlinear controller must be done in the state space domain. Consequently, effects of uncertainty have largely been limited to sensitivity analysis (Morari and Zafiriou, 1989), and little work dealing with the analysis of model uncertainties for nonlinear systems has been reported in the process control literature. In this study, a new approach is developed to deal with modeling errors. The approach integrates the constrained
* Author to whom correspondence should be addressed.
Newton-type control algorithm with a nonlinear parameter estimator where the discrepancy between predicted system outputs and measurements is used to estimate unknown or uncertain parameters. Section two provides a brief review of approaches that couple estimation techniques to process control algorithms. In our approach, the model is continuously updated by the information from the estimator; the newly updated model is employed to predict the system outputs in the future time horizon and determine the appropriate control variables to satisfy the control objective. These control and parameter estimation algorithms are described in section three. In addition, contraction mapping (local property) and descent directions (global property) for the stability of nonlinear systems are analyzed in section four when the reference model is imperfect. Here we consider both additive output uncertainty as well as model mismatch and derive conditions under which local and global convergence properties can be satisfied. Section five describes a pH control process with a number of unknown model parameters. Applied by Parrish and Brosilow (1988) to illustrate their Nonlinear Inferential Control (NLIC) approach, this example is used to demonstrate the Newton-type control law and offer a comparison with NLIC approach. In addition, significant measurement noise is added to this example to illustrate the reliability of the controller. Finally, section 6 concludes the paper by summarizing the analysis and results. 2. Review of Estimation and Control Approaches While little has been done to analyze the effect of uncertainty for nonlinear process control algorithms, there are many approaches that incorporate estimation algorithms to minimize adverse effects due to uncertainty. Covering all of these is certainly beyond the scope of this study; a brief review is given here to summarize other research directions. These have been classified as adaptive
0888-5885/ 9012629- 1647$02.50/0 0 1990 American Chemical Society
1648 Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990
control, optimizing control with model identification, and nonlinear inferential control. 2.1. Adaptive Control. In recent years, there has been extensive interest in adaptive control systems that automatically adjust the controller to compensate for unanticipated changes in the process or the environment. An extensive review has been given by Seborg et al. (1986). Adaptive control is usually based on simultaneous model identification and control. Here the basic assumption is that the plant can be represented by a member of a family of linear, finite-dimensional parametric models. The process model is identified on-line, based on input and output measurements of the plant. The model parameters are usually calculated by the recursive least-squares (RLS) method. Depending on whether or not the estimated model appears explicitly in the control law, the design method can be classified as a direct or an indirect method. In the indirect approach, a process model is employed and the control calculation is based on the estimated model parameters. In the direct approach, the original process model is converted to a predictive form that allows the future process output to be predicted from current and past values of the input and output variables. Astrom and Wittenmark (1973) suggested a direct approach, the self-tuning regulator (STR). They assumed that a perfect model is obtained for each execution time. Then a feedback controller is designed to minimize the variance of the output variables. Several studies (Clarke and Gawthrop, 1975, 1979; Clarke, 1981a,b, 1984) later developed the self-tuning controller (STC). Instead of minimizing the variance of output variables, STC minimizes the variance of an auxiliary output, i.e., a linear combination of output variables, setpoint, and control variables. Even though both STR and STC are based on minimizing a quadratic cost function, STC is easier to tune and can be applied to nonminimum-phase plants. Wellstead et al. (1979) and Astrom and Wittenmark (1980) suggested an indirect approach to design self-tuning controllers based on pole placement. Here the controller can be tuned by appropriate shifting of the closed-loop poles, which allows the controller actions to be limited as desired. Unfortunately, the actual closed-loop performance is uncertain before the controller parameters are determined. However, the basic assumption of a linear model in adaptive control design can cause problems. Few plants are truly finite order and linear; modeling error due to unmodeled plant dynamics and nonlinearities can be present even though the model is updated continuously. Therefore, the robustness of the controller is an unresolved question in adaptive control design, even though the philosophy of identifying a model based on measurements of input and output variables can shed light on future research directions. 2.2. Optimizing Control and Model Identification. When nonlinear systems are involved, one can employ a parameter estimator in the control loop in order to minimize the difference between the predicted and measured system outputs in the immediate past time horizon. This is analogous to the indirect approach for adaptive control. Bamberger and Isermann (1978) developed an approach for on-line optimization of the steady-state behavior of slow dynamic processes. A nonlinear dynamic process model is identified on-line, and an updated model is used to determine the steady-state behavior of dynamic processes and search for an optimal operating condition. Also, Prett and Gillette (1979) employed a steady-state model with model parametes that are updated by using available
measurements. This approach was successfully applied to a fluid catalytic cracking operation. Jang et al. (1987) extended this approach to dynamic nonlinear process control by proposing a two-phase algorithm to control nonlinear processes with unknown parameters. In the first phase, unknown parameters and states in the proposed model are estimated on-line by using standard nonlinear programming (NLP) techniques. The updated model is then used to predict the system behavior and determine the appropriate system inputs. The control objective is formulated as an optimization problem, whose sensitivities are calculated through a set of adjoint variables, and a two-point boundary value problem must be solved at each iteration. Moreover, process constraints are not considered in their approach. Rawlings and Eaton (19%) also used an estimation approach to solve the Shell Control Problem (Prett and Morari, 1987). Here the approach is partitioned into two parts. In the control phase, the problem is formulated by using Quadratic Dynamic Matrix Control (QDMC). The difference between measured system outputs and predicted system outputs in the immediate past is used to estimate uncertain system gains and loads. The standard QDMC algorithm estimates a current disturbance in each output by subtracting the model prediction from the available measurement. Future outputs are usually predicted by assuming that the disturbance remains constant. In this way, the QDMC-type controller can control the process satisfactorily, even when some system outputs cannot be measured. Eaton et al. (1988) also discuss a similar strategy for the solution of control systems modeled by nonlinear ODE’S. 2.3. Nonlinear Inferential Control. Finally, Parrish and Brosilow (1988) propose a NLIC approach in order to accommodate inevitable modeling errors. NLIC decomposes the control problem into two parts. The nonlinear model is used to predict system outputs and to determine control variables that satisfy the control objective. Then a linear filter is employed to translate the original setpoint and the cumulative projection errors to desired system responses. Here the original setpoint is transformed to an on-line modified setpoint that can compensate for modeling errors; the controller tracks the modified rather than the original setpoint. One advantage of NLIC is that it is stable even when the control effort is saturated. Another advantage is that its structure focuses attention on the use of available measurements to infer how disturbances influence the process. However, NLIC is a modification of a linear inferential approach for nonlinear systems and, therefore, has several limitations. In particular, since nonlinear models need to be transformed to linear models by the control module in NLIC, the nonlinear model must satisfy the following properties: 1. At all operating points, the model must be linearized so that it can be written as a set of linear difference equations. 2. It must be locally asymptotically stable everywhere in the operating range. 3. The sign of the steady-state gain must not change in the operating region. This approach will be considered in section five and compared to the algorithm developed in the next section. 3. Newton-Type Controller with Estimation: A Two-Phase Approach To begin, we assume that the structure of the system is exactly known but some of the parameters are unknown or uncertain. This assumption will later be relaxed for the stability analysis. Here the parameters can include model
Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990 1649
f
l
-
I
Figure 1. Constrained nonlinear Newton-type controller with estimation.
constants such as equilibrium constants, feed composition, heat-transfer coefficients, etc. They can also represent unmeasured disturbances. The autonomous lumped parameter (MIMO) nonlinear system can be modeled by the following equations: dx/dt = f ( w , p )
(la)
Y = g(x,p)
(1b)
where x E R". is the state of the system u ( t ) E R". is the input, p E R"Pis the set of unknown parameters, and y E R"Yis the system output, with ny I nu. The nonlinear model derived from the fundamental physical and chemical laws is used directly in our control algorithm. Hence, the model will have a wider range of validity and parameters that are physically meaningful. The model must be capable of exhibiting the measured process response despite process modeling errors. This property is referred to as "model consistency" (Parrish and Brosilow, 1988) and is required so that the model response can follow the process. Intuitively, one would expect that if the model were incapable of replicating the process response, then use of the model to predict and specify process performance would not be possible. In this two-phase approach, we consider the control algorithm (phase I) and parameter estimation (phase 11) separately. A block diagram of this two-phase approach is given in Figure 1. Of course, since all the blocks involve nonlinear transformations, standard block diagram manipulations are not valid here. 3.1. Derivation of the Newton-Type Controller. Newton-type controllers based on the nonlinear model in eq 1are implemented by using a moving time horizon as presented in Figure 2. In this discrete-time formulation, the sth sampling interval extends from ts to ts+l where T - t8+l - t8 is the constant sampling interval length, xs is the state at tS,and usis the system input held constant over ( t s , t r + l ) . Here ~ ( t ~ + ~ ; t ~is ,the x ~solution , u ~ ) of eq l a at time ts+l for u ( t ) = us ( t s < t < ts+') and the initial condition x(tS;ts,xB,uS) = xS; xs denotes the state of the system at t = t s + l , i.e., x S + l : x8
= xs+l
=
X(tS+T;t*,XS,US)
(2)
Since (la) is autonomous (x(ts+T;t8,x,u)= x(tS+'+ T;tS+l,x,u), i.e., f does not depend on t explicitly), time will be dropped from the parameter list and the following convention will be used: xs
= x(T;xS,u') = X ( t s + T ; t s , x s , u s )
The derivatives of defined as follows:
xs with
(3)
respect to xs and us will be
9 s = axs/axs
(4)
rs = aXe/aus
(5)
t
/Y +
I
I
I
t s
I s+l
t s+2
s.1
PPrttimehaimn
3
t
Rsent r i
t s+3
F'uturetimehaizon
Figure 2. Moving time horizon notation.
Here y" = g(xs)is the system output at ts and its derivative with respect to xs is defined by cs p a y a + l / a x s = a g ( x s ) / a x s (6) Equations for the sensitivity functions asand recan be found in Li and Biegler (1988) and Li et al. (1990). Moreover, they can be calculated efficiently by the algorithm of Caracotsios and Steward (1985). In order to simplify the notation of a Newton-type control law, we assume for the moment that nu = ny,i.e., system (1)is square. Also, we estimate to first order that the control objective is met at time s 2, i.e., p+2- y* = g(xS+Z) - y* = g(x(xs+l,uS+1))- y* = 0 (7)
+
and expand this expression in a Taylor series around the state x = x s , u = us: 0 = g(x(T;xs+',us+')) - y* = g(x(T;x"uS)) - y* + (xs+l
- x s) +
s+lo [ l l ( x s + l - x s ,u
s T 2
u )
II 1
(8)
To compute uS+l,we solve (8)to first order and set the higher order terms in (8) to zero. Furthermore, x S + l is substituted from (2), yielding o = y+1-y* + c s w ( x 8 - XS) + (csrs)(us++1-US) (9) By solving (9) for us+1,the following Newton-type control law is obtained: us+l = us + ( c s r s ) - 1 [ ( c w ( xs x s ) + (y* - p++')](io) Moreover, by including a relaxation factor X E (0,l) in the control law, (10) becomes us+l = us x ( c s r S ) - l [ ( c w ( x s- xs) + (y* - ys++')l (11) The relaxation factor not only changes the shape and speed of the system response but also expands the region of stability of the closed-loop system. In fact, we see in section 4 that this relaxation factor can also be used to control the effect of model mismatch on the controller. The Newton-type control law can also be extended to deal with input and output constraints. Here the control objective described above can be written as follows:
+
min Ily(T;xs+l,u8+1) - y*l12 ua+l
(12)
If the first-order approximation based on x s and us is used to represent y ( T ; ~ ~ + ~ , uy~becomes +l), y(T;xs+1,u8+1) = ys+2 ys+l + C S ~ ~ ~ ( T ; P , U ~ )+ - X( c ~ e) r s ) ( u s + l -US) (13)
1650 Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990
Substituting (13) into (12) and rearranging yields the objective function 1 min CTAus + - ( A U ~ ) ~ H A U ~ JUS 2 Aus =
us+l
-
us
(14)
where
C T = - [ C T ( x S - xs) + (y*- y s + l ) ] T ( ~ s r s )(15) H = (csre)T(csrs) (16) Now the above Quadratic Programming (QP) problem without any constraints is identical with the Newton-type control law. Note that the Hessian of the objective function is positive semidefinite, which means that a global minimum value can be found. Moreover, within a QP formulation, nonsquare systems can also be controlled by this algorithm; i.e., ny does not need to be equal to nu. With quadratic programming and the quadratic form of the objective function, linear equality and inequality constraints on control variables can be handled easily. In this case, problem (14) becomes 1 min CTAu8+ - ( A U ’ ) ~ H A U ~ Au3 2
subject to
aLIA(us + Au) Ia’ (17) As with the unconstrained control law of eq 11, a relaxation parameter is also used so that the controller move is given by uS+l= us + X(Aus),and X is chosen so that the deviation from the setpoint is reduced. Li and Biegler (1989) also extended this control law to the multistep case, where optimal controller moves are predicted over a moving horizon. Now when modeling errors exist, nominal parameters in the model differ from those in the process. Since the structure is assumed to be exact, these are the only source of modeling error, and differences in predicted and measured system outputs are used to update the nominal values of these parameters. Further we assume that the model parameters in (1)vary slowly and remain constant over the time horizon considered. If, on the other hand, variation in the parameters is sufficiently large over the time horizon, then one could attempt to postulate some dynamic model for the variation. As seen above, similar approaches are proposed by other researchers (Bamberger and Isermann, 1978; Prett and Gillette, 1979; Jang et al., 1987). Here, if the model structure of the system is exactly known, this approach can exploit the measurement data to get an optimal set of parameters that is as close as possible to the actual plant values. Details of the parameter estimation algorithm will be discussed in the next subsection. Applying the controller in a moving horizon framework, parameter estimation can then be implemented over past measurements, which consist of m equally spaced sampling intervals as shown in Figure 2. The system outputs a t the end of every subinterval in that horizon, Y+l, y,y-l, ..., ysf2-’”, are measured, recorded, and used by the estimator to update the nominal parameters. As time progresses, replaces the measured a new measurement at time tS= F1 system outputs of the furthest subinterval in the past time horizon for parameter estimation. After implementing parameter estimation for a few time horizons, the model parameters may be sufficiently accurate to predict the system outputs if these parameters are varying very slowly. For example, the following criterion can be employed to decide whether the estimation phase can be eliminated:
m
(ys+2-[- p+2-1)TA~(ys+2-1 - p+2-1 < ) -
t
118)
L=l
Here, c is a small positive number, A, is a weighting matrix to reflect the importance of this discrepancy in the past time horizon, p+2-r is the predicted system output, and p+z-r is the measured system output. If the above condition is satisfied, the model parameters are accurate enough to be used in the control phase and the parameter estimation phase can be suppressed. Finally, the sampling interval in the past time horizon need not be equal to the sampling interval in the predictive time horizon and depends on the available measured data. A more detailed description of the parameter estimation formulation is considered in the next subsection. 3.2. Parameter Estimation. The objective of parameter estimation is to find a set of optimal parameters as close as possible to the actual values in the process. However, since the actual unknown parameters are not available and only the system outputs can be measured, an alternative is to find a set of parameters that minimizes the difference between the predicted and the measured system outputs in the past time horizon. If the model structure is exact, this discrepancy is an indication of how well the model with the nominal values of these parameters represents the process. The following residuals, for m steps in the past time horizon, are defined to represent this difference: eS-i = ys+2-i - y + 2 - i i = 1, 2, ..., m (19) ys+2-i , the predicted system outputs, can be calculated from model (1)based on a set of nominal parameters. The on-line measured system outputs Ys+2-i, are recorded. To further simplify the problem in this study, the discrepancy is assumed to be caused by inaccurate values of these model parameters and measurement noise. Most of the measurement noise consists of random variables that can be represented by different types of distribution functions. However, as indicated by Bard (1974), the most popular distribution function used to represent the measurement noise is a normal distribution function, due to the fact that it often describes measurement noise simply and accurately. Here the mean vector, e, and covariance matrix, P, are sufficient to define a multivariate normal distribution function:
The most general objective for parameter estimation is the maximum likelihood estimate. Here the random variable e is substituted in the joint probability distribution function by the actual residual values, es. As the residual values are functions of the model parameters, the maximum likelihood estimate of the parameters can be obtained directly by maximizing the resulting objective function. To simplify this problem, it can be shown (Bard, 1974) that if the residuals can be represented by a normal distribution and the covariance matrix (P = W-l) is known, the values of parameters can be determined by a weighted least-squares method m
min @(p)= min ~ ( e s - i ) T W i e s - i
P
P
i-1
(21’1
where Wi is a positive definite matrix of dimension ny x nY. Finally, if the covariance matrix is diagonal and has equal elements, then the unweighted least-squares method
Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990 1651 possesses optimal statistical properties (Bard, 1974). The mathematical representation of the unweighted leastsquares method is
Moreover, meaningful parameter estimates in engineering problems are usually constrained to upper and lower bounds, and consequently, the parameter estimation problem can be written in the following form: min + ( p )
P
subject to (23) p* Ip I p* Caracotsios (1986) developed a Successive Quadratic Programming type method to solve problem (23). In his approach, the objective function is expanded up to second order, and the quadratic approximation function serves as a temporary objective function subject to rectangular Constraints on parameters. Solving these quadratic programs successively then leads to a global minimum value. This approach was implemented in the GREG (General REGression) computer code, which was used in this study. Since the code is designed specially to solve nonlinear parameter estimation problems such as problem (7),it has some specialized features: In order to detect linear dependencies among observed system responses, an eigenvalue-eigenvector analysis on the moment matrix of residuals is performed. To avoid calculating the second derivatives, the Hessian is constructed in a similar manner as in the Gauss-Newton method; i.e., only the first-order sensitivity functions are used. The first-order sensitivity function is efficiently calculated by decoupling the model equation and sensitivity function equations (see Caracotsios and Stewart (1985) for details). To handle a singular Hessian at any stage of optimization, the numerical rank of the Hessian matrix is carefully determined by introducing lower bounds on the normalized diagonal pivots and an appropriate positive scalar quantity is selected for the pseudo-Hessian. To prevent the trust region from being too conservative, the region is automatically adjusted (enlargement or contraction) by checking descent properties. Since the quadratic expansion is not valid for all of the values of the parameter vector p in the feasible region, one can impose constraints on the parameter search direction vector to assure the validity of the quadratic expansion. Thus, additional rectangular constraints are imposed that form the trust region, which varies over the course of the optimization. To increase the computational efficiency, a generalized pseudoinverse is constructed to find a search direction for the next iteration. Having described both the control strategy and the parameter estimation procedure, we are now ready to state the combined algorithm for this two-phase approach. 3.3. Two-PhaseAlgorithm. 1. (a) Choose the number of steps in the predictive time horizon, I, and sampling interval, T. (b) Choose the number of steps in the past time horizon, m. (c) Guess the values of the unknown model parameters. 2. If the number of executed steps is larger than m (2np), go to step 3. (It is necessary to accumulate m measurements before the parameter estimator can initiate.) Otherwise, go to step 4.
3. (a) Determine the controller move from either eq 11 or eq 17 by using the updated model. Determine the step size that leads to a reduction in the predicted output deviation from the setpoint. Execute the (first step) controller move. (b) Measure the system output p+2and replace the furthest past measured system outputs p+l-"' in the record. 4. Check the model accuracy. If C$l(e"i)TAi(es-i) It, where 6 is small positive number, and Ai is a weighting factor matrix, go to step 3. Otherwise, go to step 5. 5. (a) Use recorded system output measurements in the past time horizon to estimate unknown model parameters in the estimation phase. (b) Update the model by using newly estimated model parameters. (c) Go to step 3, with s=s+l. 4. Stability Analysis Stability analysis for nonlinear systems in the presence of model uncertainty is generally quite difficult and has largely been confined to sensitivity analyses. For nonlinear controllers, Economou et al. (1986) used a nonlinear IMC structure in order to demonstrate a zero offset property if the steady-state gains of the plant and the model are equal. This controller, however, requires a nonlinear model inverse, which may not generally be available. Newtontype controllers, for example, only provide an inverse of the model linearization. They also showed more general stability results by applying Zames' small gain theorem to the IMC structure. Li and Biegler (1988) also applied the small gain theorem to constrained Newton-type controllers and showed that input constraints can be substituted for the filter block in the IMC structure in order to maintain stability. However, while these results lead to useful qualitative arguments, the computation of nonlinear gains is generally too difficult to serve as a quantitative analysis tool. Using concepts for convergence of methods for solving nonlinear equations, Li et al. (1989) and Li and Biegler (1989) established stability results for single-step and multistep Newton-typecontrollers, respectively. First they showed that if a sufficiently large sampling interval is chosen for open-loop-stable systems, a region of attraction always exists, ]ly(t)- y*ll 5 r, for which the contraction mapping theorem applies. Consequently, monotonic convergence to the setpoint is guaranteed from any point within this region. We term this result a local property for Newton-type controllers. For linear systems, this region of attraction, when it exists, is R n y , while for nonlinear systems it may be confined to only a small neighborhood around the setpoint. Consequently, Li et al. (1990) included a relaxation parameter in their Newton-type controller that ensures that setpoint deviations, IIys+' - y*ll, decrease with each controller move. Existence of a nonzero relaxation parameter that ensures this decrease can be guaranteed if the controller move is a descent direction for Ily(t) - y*ll. For Newton-type controllers, Li et al. (1990) proved that descent directions can always be found if the sampling interval is sufficiently large, the process is open loop stable, and Ip+'- y*ll 2 6. Here 6 is an arbitrarily small, positive constant. We term this property the global property for Newton-type controllers. Thus, with the global property, the setpoint deviation is reduced for all y ( t )where ily" - y*ll 2 6, while the local property guarantees convergence to the setpoint within the region of attraction, Ib(t)- y*ll I r. Consequently, stability can be guaranteed by combining these two properties and choosing 6 well within the region of attraction. Li and Biegler (1989) showed that these stability properties apply
1652 Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990
-e- r' and regions defining the local and global convergence properties may not overlap for sufficiently large M , thus leading to instability. This concept is illustrated in Figure 3. From the proof, we also note that by reducing the step size, the terms describing the uncertainty also can be reduced. Consequently, the region of attraction can be increased by taking smaller controller moves (effectively detuning the controller) and sacrificing control performance. 4.2. General Model Uncertainty. We now extend these results to general model uncertainty. In this case, less precise results are shown because of the uncertain structure of the model. Here we lump the model mismatch into two types. Additive output error, unmeasured disturbances, and other terms that lead to steady-state mismatch are lumped into the 11 term, which is assumed to be bounded. Dynamic mismatch terms, which vanish a t steady state, are observed in the calculation of the sensitivity functions. Here we require that
We now consider the effect of this uncertainty description on the stability properties. The following theorem gives a relaxation of the global property. Theorem 1: Global Property under Output Uncertainty. If an open-loop nonlinear system is asymptotically stable in the large and has bounded additive output uncertainty, i.e., 111111 IM , where 7 is a random variable and M is some positive number, then the descent direction of the closed-loop system is guaranteed for IlY" - y*11,1 6 + M (where 6 > 0) with a sampling interval sufficiently large. Here we have llys+2
- y*ll -
IIyS+l -
6h -7
y * J (I
(24)
for some nonzero value of X E (0, 1). A proof of this theorem can be found in Li (1989; theorem 5.1). Now, in order for the system output to converge to the setpoint, the region of attraction must be large enough to include 6 M . In the next theorem, we consider the effect of additive uncertainty on the local
+
(25a)
rs = Ts
(25b)
lim
T--
+
Y = d X , P ) + 17
lim (Ps= @
T-..
lim Cs = C s
(25~)
T--
where y , P,P,and Csare output and sensitivjty functions calculated from the model and 3, @,P,and Csare output and sensitivity functions from the actual process. Note that the terms in eq 25a tend to zero, while those in (25b) and (25c) tend to their steady-state values. In addition, we assume that higher order sensitivity functions for the process model also converge to the actual steadystate values for the process. This classification now allows us to make some observations about the stability of Newton-type controllers for the general mismatch case. We consider first the relaxation of the global property. Theorem 3: Global Property for General Model Mismatch. If an open-loop nonlinear system has the properties (1) it is asymptotically stable in the large, (2) the model sensitivity functions approach the actual values when T approaches infinity, and (3) 11q11 I M, where 9 is random noise as well as other steady-state mismatch terms, and M is a positive number, then the descent direction of the closed-loop system is guaranteed for IIy"" - y*ll 2 6 -b M , where 6 > 0, with a sampling interval sufficiently large. That is, 11ys+2
- ?/*I1- I l j i S + '
- y*ll 5
6X -2
(26)
for some nonzero value of h E (0, 11.
Here the bar is used to indicate the variables representing the actual quantities in the process. The proof of this theorem can be found in Li (1989; theorem 5.2). This
Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990 1653 theorem shows that there exists a sufficiently large but finite sampling interval for which the descent property is guaranteed. Note here that when model mismatch exists, we require for lP+l-y*ll < 6 + M that the region of attraction must again be increased to include the additional region affected by the measurement noise and model mismatch. The effect of model mismatch on the region of attraction, however, is difficult to determine unless the actual process model is known. Here the region of attraction for the local property is defined by the following problem: Find r'such that for all in llvll IM
I
-
Flowrate 0 basesolu.
E Figure 4. Neutralization process.
- y*ll 5 r'. Here the measurements and for all ~ ( tin) control law are given by
Q=Y+v us+l
= us
+ X ( C ~ ~ ~ ) - - ~ [ (-C ~ @ +~ y* ( X -~ p+1] XS+~)
Differentiating this control law leads to the following controller terms in eq 27
x(csrs)-l[cs@s - Cs@] A(
c srs)-l[csrs- cT I
It is interesting to note here that both T and X have an effect on the local property. First, increasing T will cause the differences in the dynamic mismatch terms (given in square brackets) to vanish. Also, by reducing the step size, the contribution of model mismatch terms is lessened. Finally, general model mismatch can be reduced through parameter estimation. In this way both steady-state and dynamic mismatch terms are minimized for a given T. Although, the effect of this task is hard to assess in a quantitative manner, parameter estimation can be very beneficial in order to maintain a large enough region of attraction for stability. In summary, we see that the stability analysis for the general model mismatch case confirms qualitatively some intuitive conjectures. First, robust stability properties can be improved by increasing the sampling interval and decreasing the step size for the controller moves. Both actions serve to detune the controller and lead to trade-offs between performance and stability. Second, parameter estimation reduces the model mismatch terms and thus helps restore the stability properties determined for the nominal plant model. In the next section, we consider a nonlinear process with parametric mismatches and output uncertainty in order to demonstrate the performance of the constrained Newton-type controller. 5. Test Problem
In this section, a pH control problem is simulated to demonstrate the effectiveness of the nonlinear constrained Neyton-type controller with parameter estimation. This single-input, single-output (SISO) problem was taken from Parrish and Brosilow (1988) and is two dimensional, having one ODE and one algebraic equation with three unknown parameters and one unknown state in the model. Only the
outlet pH value measurement is available; the control objective is to maintain the pH value of the outlet flow at a desired level. Here, two cases are examined. In the first case, the process is simulated by our proposed control algorithm under the conditions described by Parrish and Brosilow. The simulation results are plotted with those simulated by the NLIC strategy. In the second case, white noise is introduced in the measurement. It is demonstrated that the control algorithm still works even though it takes more time to reach the setpoint. On the other hand, NLIC fails to satisfy the control objective. The pH process considered in this example is illustrated in Figure 4. The process consists of a CSTR into which a feed stream of unknown composition flows. The solution in the CSTR is maintained at the desired pH by addition of a basic solution. The feed to the CSTR is assumed to be a monoprotic acid (HA) and its soluble salts (A-). The concentrations of the acids and the salts in the CSTR are unknown. The neutralization is carried out by using a solution of a strong base of concentration C,at a flow rate of m (the control variable). A model of the system is derived based on the following assumptions: 1. The CSTR is perfectly mixed. 2. The acid-base equilibrium reaction is fast. 3. The equilibrium of the mixture of acids and salts is effectively represented by a single expression. The third assumption is based on the fact that mixtures of acids and bases often produce titration curves similar to those produced by a single acid and its salts (Gusstafson, 1982; Gusstafson and Waller, 1983). The process has been modeled by the following equations: H+ + N+ = A- + OH- (charge balance) (28a) HA
+ A- = U
(an unknown)
(28b)
(H+A-)/ (HA) = K, (acid-base equilibrium constant) (2W
Equations 28a-28d may be reduced to
The parameters of the process are N+,the concentration of the cation (base) in the CSTR; Kw,the equilibrium constant for water; K,, the equilibrium constant of the acid mixture; F, the feed rate to the tank; m, the flow rate of the base (control variable); Cr, the concentration of the
1654 Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990 Table I. Values parameters F K,
u c, v
of Parameters real values 13 L/min 0.001 0.13 mol/L 1.0 equiv/L 90 L
for the Test Problem NLIC NLCNC 13 L/min 10 L/min 0.0005 0.0005 0.10 mol/L 0.10 mol/L
15
rate)
+ ( p z + x 1 ) x q 2 + ( X G ~ - 10-l~- p 2 p 3 ) x 210-'4p, =
o
(29b)
11 10
9 t 10
0
1
le3
..-..--.-p2
-
(Equil. cons.
*'
(30)
A new set of parameter values is obtained by the estimator, which is utilized to update the model. Also, in (29) the initial value of cation concentration, xl, is unknown and needs to be identified. Therefore, four model param-
(estimated) (real)
10
Min.
Figure 6. Estimated equilibrium constant of acid and base versus time.
1 .....-...
-
0 I u I3.0
0.5pg Ip 2 I2 . 0 ~ 8
RP2
4
4.94 0
(29~)
where xl is N+, x q is H+,p1 is F , pq is K,, p3 is U , and u is m. This model is highly nonlinear, with a gain variation of 2 or more orders of magnitude over the operating region. The gain variation is due to several factors, among them changes in operating point pH, buffering level U , and equilibrium constant K,. The NLIC strategy directly utilizes the process model described above. The details of implementing NLIC on this model are described in Parrish and Brosilow (1988) and are not considered here. It should be noted that this algorithm is subject to tuning parameters that can affect performance; we used the tuning parameters in their study. It is also worth noting that in their work the acid flow rate of the process is initially set to 10 L/min. However, if the acid flow rate in the model is significantly different from the process, we found that NLIC fails to converge. Therefore, the true value of acid flow rate (13 L/min) is used by NLIC analysis in this study. In our control algorithm, the GREG parameter estimation package (Caracotsios, 1986) is used to estimate the unknown model parameters. The sampling interval is set to 0.1 min. In the past time horizon, six measurements ( m = 6) are recorded before estimating these parameters. Here the superscript "0" is used to indicate that these are the nominal values of parameters in the current model. The following set of constraints is used by the estimator: 0.5pp 5 p 1 5 1.5pp
Min
Figure 5. Estimated acid flow rate versus time.
0.15
y = -(log xz)
0 . 5 ~ :Ip3 I1 . 5 ~ :
Pi (ertimated) RP1 (real)
P1
base solution; U , the unknown concentration of the acid and its salts; and V , the volume of the tank. All parameters must be positive, and upper and lower bounds are imposed on the control variable, m. Process model parameters are given in Table I. To represent model mismatch, the model parameters in this table differ from the true values. In column NLIC, the values of the parameters are used to simulate the model and determine the control trajectory with the NLIC strategy. In column NLCNC, the model parameters are farther away from the true values than with NLIC. Because of the estimation feature within our control algorithm, these values become the initial guesses of the model parameters for the Newton-type controller. As our algorithm proceeds, newly estimated parameters are employed to update the model on-line. Using the numerical data given in Table I, model (28) can be written as follows: dxl/dt = 0.0111(~- (p1 + u)x1) (294 x23
- - - _-- -.
-
0.12 0.11
-
1:
p3 RP3
(ertimated) (real)
j
i
0
10
Min.
m
Figure 7. Estimated equilibrium constant of the mixture of acids and salts versus time.
eters need to be estimated by using a single set of measurements, Le., the pH value in the tank. Equation 29 is a system of differential-algebraic equations (DAE's), which can be converted to an ODE system. Here, the desired value of x2 (H+ concentration) can be calculated from (29c). After the model is updated by a set of newly estimated parameters, the value of xl, the concentration of cation, can be computed from (29b). Here the single-step constrained Newton-type control strategy is used. Figures 5-7 show estimated values of three estimated parameters versus time. The actual parameter values used to simulate the process are also plotted as references. We observe that the estimated values of these parameters converge to the actual values after a few sampling times and remain there. Figure 8 plots system outputs simulated by using our control algorithm and NLIC versus time. From this figure, it is clear that our control algorithm performs significantly better than NLIC because the former totally eliminates oscillations and reaches the setpoint much faster. After a few time intervals, an almost perfect model is obtained in our control algorithm (see estimated values of the model parameters in Figures 5-7) and a nonlinear controller is
Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990 1655
- ---
"1 Y
' P1 (estlmrtrd)
P1
-.-*-
(pH vduo In tho tank)
(Acid flow rate) NLIC NLCNC
r
i
I i
-Setpoint lo 6
4
0
.
.
,
.
.
2
,
.
.
4
10
, 6
. Min.
Min.
Figure 10. Estimated acid flow rate versus time (with white noise).
Figure 8. pH value in the tank versus time.
'1
-.-.-
NLIC NLCNC 0.002
P2
U
@qui. cons.
(base solution flow rate)
of icida P bw:)
0.001
1 0
10
Min. 0
Figure 9. Basic solution flow rate versus time.
used, which fully represents the system dynamic characteristics. Therefore, a control variable can be found to satisfy the control objective, Le., minimize the norm of the differences between the linearized system outputs and setpoint. This argument is supported by the fact that the flow rate of base (the control variable) used by our algorithm always remains a t the maximum feasible value, 3.0 L/min, until the system output reaches the setpoint showed in Figure 8. On the other hand, NLIC does not utilize the pH measurement to estimate the unknown model parameters. In fact, a model with significant modeling error is always used to determine the appropriate control variable to satisfy the control objective. The filter block in NLIC only partially compensates for the modeling error. This is why the true value of acid flow rate, F (see Table I), needs to be used in order to make NLIC work. Since the filter can only locally compensate for modeling error, it fails to compensate for the modeling error if the initial difference between the true value and model value of acid flow rate is too large. However, when parameter changes are slow, the filter with local compensation still allows convergence to the setpoint. In addition, the discrete form of the model used by NLIC cannot fully represent the dynamic characteristics of the process. Therefore, NLIC is unable to use the maximum control power (in this case, the maximum value of the base flow rate), as seen in Figure 9, and the system output converges to the setpoint a t a slower rate. Measurement Noise Case. The pH process is also examined with measurement noise. Here the process model remains the same except that white noise is added to the H+concentration measurement: x Z m = ~ ~ ( 1+. 0N ( p , P ) )
(31)
2
1
B
Min.
Figure 11. Estimated equilibrium constant of acid and base versus time (with white noise). 0.16
1
-
P3 (estimated) P3 (real)
-.I-.
P3
0.14
;
' d ' ,
Mln.
Figure 12. Estimated equilibrium constant of the mixture of acids and salts versus time (with white noise).
where N is a normal distribution function, p = 0, and h2 = 0.01. Here, x Z m is a measured concentration, and x 2 is the actual value of concentration in the process. For this case, NLIC is unstable and is not plotted. Figures 10-12 show estimated values of three model parameters. Compared to Figures 5-7, it can be clearly observed that it takes longer for the estimated parameters to converge to the actual values. Note that the estimated values of parameters remain constant and only have a slight offset from the actual values after about 3-4 min. Figure 13 shows actual and measured pH values versus time. In the presence of measurement noise, it takes longer
1656 Ind. Eng. Chem. Res., Vol. 29, No. 8, 1990
Y (pH vslue)
a
2
4
6
Min.
Figure 13. pH value in the tank versus time (with white noise).
3i 41
,(base :.,no. solution
+ u
2
tainty, but because little can be said in general about the uncertain structure of the actual process, these results are less precise. In both cases, however, stability properties can be improved by increasing the sampling interval or reducing the relaxation parameter. However, both also serve to detune the controller. Also, we note that adverse effects of model mismatch on stability can be minimized by incorporating a parameter estimation algorithm. Finally, a small nonlinear control example is considered to demonstrate the effectiveness of this two-phase strategy. Here, a two-dimensional (SISO)pH control problem with three unknown model parameters and an unknown state is simulated by the proposed control algorithm. It should be noted, however, that MIMO applications have also been considered with these controllers (Li and Biegler, 1988, 1989; Li et al., 1990). In this example, only one measurement is available for estimating unknown model parameters and identifying the unknown state variable. The simulation results from using our control algorithm are compared with those from an NLIC strategy. It was observed that our algorithm performs much better than does NLIC. In addition, when the measurement has white noise, NLIC fails to converge, while our control algorithm drives the system to the setpoint and gets offset free control. Acknowledgment
O t
a
. .
, 2
.
.
, 4
.
. ,
.
6
Mln.
Financial support from the National Science Foundation through a Presidential Young Investigator award is gratefully acknowledged.
Figure 14. Basic solution flow rate versus time (with white noise).
to reach the setpoint, and oscillations may result because of modeling error. However, the estimator updates the unknown model parameters on-line, and the model is increasingly accurate. Thus, the controller is able to drive the system to the setpoint and get offset free control, which can be seen in Figure 13. Figure 14 shows the basic solution flow rate, control variable m, versus time. 6. Conclusions In this study, the assumption of a perfect model for Newton-type controllers is relaxed. Here a two-phase approach is developed to deal with modeling error. In the first phase, the updated model is used to predict system outputs and to determine appropriate control variables to optimize the control objective. In the second phase, the discrepancy between predicted and measured system outputs in the past time horizon is used by a parameter estimator to determine optimal values for the unknown process parameters. The Newton-type algorithm is also analyzed for stability properties in the presence of model uncertainty. In previous studies, nonlinear operator theory was used to show local and global convergence properties as long as a perfect process model is assumed. Local convergence assumes a contraction mapping in some region of attraction about the setpoint, while the global property guarantees improvement toward the setpoint by providing descent directions for the controller moves. Here these properties are relaxed by treating two cases of model uncertainty. First, only additive output uncertainty is considered. Here descent directions can be found only in a region farther away from the setpoint. Moreover, the region of attraction shrinks for nonlinear systems. Consequently, regions where both global and local properties hold may not exist for large additive uncertainties and instability may result (see Figure 3). Similar results are observed in the second case for general model uncer-
Appendix Proof of Theorem 2. The region of attraction for which a contraction mapping exists is defined by the following problem: Find r such that
for all y ( t ) in I l y ( t ) - y*ll 5 r. Here I I ( X ~ , ~ ~ += ~ ,uSfl U ~ )is the control algorithm. For the Newton-type controller with a perfect model, the terms in (A.l) are given by ax( T;X ,u) ax
ax ( T;W)
=
@Is
=
rs
For additive output errors where y+'= ys+' + 7: the following terms in ( A . l ) become
Ind. Eng. Chem.Res., Vol. 29, No. 8, 1990 1657
a[ (csrs)-l] ax
ax8
a[ (csrs)-l] aus
v v
and the presence of the random variable leads to the following definition for the region of attraction: Find r'such that for all 7 in 11711 IM.
for all y ( t ) in Ily(t)- y*ll Ir'. Since the region of attraction, r, given by (A.l) applies only for 9 = 0, we see that r'cannot be greater than r , as (A.2) must hold for all bounded values of 9. Hence, this region of attraction must be more conservative than the one defined by (A.1).
Literature Cited Astrom, K. J.; Wittenmark, B. On Self-Tuning Regulators. Automatica 1973, 9, 185. Astrom, K. J.; Wittenmark, B. Self-Tuning Controllers Based On Pole-Zero Placement. IEE Proc. 1980, 127, 120. Bamberger, W.; Isermann, R. Adaptive On-Line Optimization of Slow Dynamic Processes. Automatica 1978, 14, 223. Bard, Y.Nonlinear Parameter Estimation; Academic Press: New York, 1974. Caracotsios, M. Model Parametric Sensitivity Analysis and Nonlinear Parameter Estimation: Theory and Application, Ph.D. Dissertation, University of Wisconsin-Madison, March 1986. Caracotsios, M.; Stewart, W. E. Sensitivity Analysis of Initial Value Problems with Mixed ODES and Algebraic Equations. Comput. Chem. Eng. 1985,9, 359-365. Clarke, D. W. Introduction to Self-Tuning Controllers. Self-Tuning and Adaptive Control: Theory and Applications; Peregrinus: London, 1981a. Clarke, D. W. Implementation of Self-Tuning Controlers. SelfTuning and Adaptive Control: Theory and Applications; Peregrinus, London, 1981b. Clarke, D. W. Self-Tuning Control of Nonminimum-Phase systems. Automatica 1984,20, 501. Clarke, D. W.;Gawthrop, P. J. Self-Tuning Control. R o c . IEE 1975, 122, 929.
Clarke, D. W.; Gawthrop, P. J. Self-Tuning Control. Proc. IEE 1979, 126. 633. Eaton, J. W.; Rawlings, J. B.; Edgar, T. F. Model-Predictive Control and Sensitivity Analysis for Constrained Nonlinear Processes. Proceedings of the IFAC Workshop on Model Based Process Control, Atlanta, June 1988. Economou, C. G. An Operator Theory Approach to Nonlinear Controller Design. Ph.D. Dissertation, California Institute of Technology, June 1985. Economou, C. G.; Morari, M.; Palsson, B. 0. Internal Model Control. 5. Extension to Nonlinear Systems. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 403. Gusstafson, T. K. Calculation of the pH Value of a Mixture of Solutions-An illustration of the Use of Reaction Invariants. Chem. Eng. Sci. 1982, 39 (9). Gusstafson, T. K.; Waller, K. V. Fundamental Properties of Continuous pH Control. ISA Trans. 1983,22 (1). Jang, S. S.; Joseph, B.; Mukai, H. On-Line Optimization of Constrained Multivariable Chemical Processes. AIChE J. 1987,33 (l), 26. Li, W . 4 . A Mathematical Programming Approach for Control of Constrained Nonlinear Systems. Ph.D. Dissertation, Carnegie Mellon University, Pittsburgh, PA, 1989. Li, W. C.; Biegler, L. T. Process Control Strategies for Constrained Nonlinear Systems. Ind. Eng. Chem. Res. 1988, 27, 1421. Li, W. C.; Biegler, L. T. Multistep, Newton-type Control Strategies for Constrained, Nonlinear Processes. Chem. Eng. Res. Des. 1989, 67, 562. Li, W. C.; Biegler, L. T.; Economou, C. G.; Morari, M. A Newton Type Control Strategy for Nonlinear Systems. Comp. Chem. Eng. 1990, in press. Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989. Parrish, J. R.; Brosilow, C. B. Nonlinear Inferential Control. AIChE J. 1988, 34 (4), 633. Prett, D. M.; Gillette, R. D. Optimization and Constrained Multivariable Control of a Catalytic Cracking Unit. Presented at the AIChE National Meeting, Houston, TX, 1979. Prett, D. M.; Morari, M. Shell Control Problem. In Shell Process Control Workshop;Prett, D. M., Morari, M., Eds.; Butterworths: Stoneham, MA, 1987, pp 225-250. Rawlings, J. B.; Eaton, J. W. Optimal Control and Model Identification Applied to the Shell Standard Control Problem. Presented at the Second Shell Process Control Workshop, 1990. Seborg, D. E.; Edgar, T. F.; Shah, S. L. Adaptive Control Strategies for Process control: A Survey. AIChE J. 1986, 32 (6), 881. Wellstead, P. E.; Prager, D.; Zanker, P. Pole-Assignment Self-Tuning Regulator. Proc. IEE 1979, 126, 781.
Received f o r review July 5 , 1989 Revised manuscript received March 19, 1990 Accepted March 28, 1990
