Adaptive Model Predictive Control of Multivariable Time-varying Systems

2708

Ind. Eng. Chem. Res. 2008, 47, 2708-2720

Adaptive Model Predictive Control of Multivariable Time-varying Systems Srinivas Karra, Rajesh Shaw, Sachin C. Patwardhan,* and Santosh Noronha Dept. of Chemical Engineering, Indian Institute of Technology, Bombay, Powai, Mumbai, 400076, India

In this work, we develop a novel adaptive model predictive control (AMPC) formulation for multivariable time-varying systems. A two-tier modeling scheme is proposed in which the deterministic and stochastic components of the model are updated on-line by two separate recursive pseudolinear regression schemes. To incorporate good long-range prediction capability with respect to manipulated inputs, we propose to identify a parametrized form of an output error (OE) model using the input-output data. To account for unmeasured disturbances, the residuals generated by the OE model are modeled as ARMA processes. To address the admissibility issue and the parameter-drift problem in on-line estimation, we propose to use a combination of a novel constrained recursive estimation scheme and the conventional recursive methods for model parameter estimation. An interesting feature of the proposed modeling scheme is that it allows separation of stationary and nonstationary components of the unmeasured disturbances. The deterministic and stochastic components of the model are then combined to form a linear time-varying state-space model, which is then used to formulate the predictive control problem at each sampling instant. A novel feature of the proposed controller is that we develop an infinite-horizon MPC formulation with a time-varying objective function in the adaptive-control framework for dealing with the grade-transition problems in continuously operated plants. The applicability of the proposed approach to the control of semibatch and continuously operated large-scale industrial processes is demonstrated by carrying out simulation studies on the production of penicillin and on the benchmark Tennessee Eastman control problem. The proposed modeling and control scheme is also validated by conducting experiments on a benchmark quadruple-tank setup. The analysis of the simulation and experimental results reveals that the proposed two-tier modeling scheme successfully captures time-varying system dynamics with respect to manipulated inputs as well as unmeasured disturbances over a wide operating range. The proposed AMPC scheme is able to achieve tight control of time-varying semibatch processes and is capable of managing large transitions in the operating conditions of a continuously operated complex multivariable process. 1. Introduction Model predictive control (MPC) has established a new era in the control of complex unit operations in the process industry over the past two decades.1,2 The key component of any MPC scheme is the model used for carrying out on-line predictions of future plant behavior. In many continuously operated plants, an MPC is implemented as a part of an on-line optimization scheme, where the operating point keeps shifting depending on the market demands, raw-material availability, energy consumption, and so forth. The linear prediction models used in most of the industrial MPC schemes, on the other hand, are typically developed only once in the neighborhood of a nominal operating point, at the beginning of the controller implementation. Such linear-model-based MPC schemes may prove ineffective if large transitions have to be made in the operating conditions. Also, as the plant dynamics changes with time, a large mismatch develops between the model and the process. Though MPC formulations are equipped to deal with a certain degree of plantmodel mismatch, the presence of large plant-model mismatch can significantly deteriorate the controller performance. In recent years, there has been a tremendous growth in the use of batch and semibatch reactors in the chemical industry due to the demand for a large number of specialty chemicals. Control of a batch/semibatch reactor along the optimal set-point trajectory is one of the key problems in this area. These systems typically exhibit highly nonlinear and time-varying behavior, which makes the associated control problems difficult to handle. The majority of control schemes proposed in the literature to * To whom correspondence should be addressed. E-mail: sachinp@ iitb.ac.in. Fax: +91-22-25726895.

deal with these systems are based on models derived from first principles.3 However, the development and validation of a mechanistic model is, in general, a difficult and time-consuming process and requires a domain expert. On the other hand, employing a control scheme based on a black-box model, which can capture time-varying characteristics of a fed-batch system, can prove to be a relatively easy and economically attractive alternative option in many situations. The problem of maintaining MPC performance in the face of changing operating conditions and/or plant characteristics has been dealt mainly by (a) incorporating robustness at the controller design stage,1 (b) employing nonlinear models for prediction,4 or (c) updating the parameters of the linear prediction model on-line. Whereas the last alternative is attractive from the viewpoint of implementation,14 it has received relatively less attention in the MPC literature when compared to the first two alternatives. Qin and Badgwell,2 in their recent review of industrial MPC algorithms, point out that only two adaptive MPC algorithms have reached the market place despite a strong market incentive for self-tuning MPC. Over last two decades, a number of adaptive MPC schemes have been proposed in the literature. The most prominent among these is the generalized predictive control (GPC) proposed by Clarke et al.5 Adaptive GPC formulations have been employed for controlling continuous as well as semibatch processes.7-9 While controlling a continuously operated system over a wide operating range, one of the major concerns is handling of the nonstationary component of the unmeasured disturbances. Similarly, while modeling a semibatch process using linear perturbation models, there is need to account for that fact that there is no single steady state as the process drifts continuously.

10.1021/ie070823y CCC: $40.75 © 2008 American Chemical Society Published on Web 03/26/2008

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2709

The use of models with integrated white-noise disturbance (ARIMAX models) is quite typical in most of the adaptive GPC formulations developed for such systems. Whereas this assumption helps in introducing integral action in the controller formulation,15 it poses some practical difficulties at the model identification stage. Under the integrated white-noise assumption, the parameter estimation problem has to be formulated using the differenced variables, that is, using ∆y(k) ) y(k) y(k - 1) and ∆u(k) ) u(k) - u(k - 1). This high-pass filtering of input-output data can lead to large bias errors at low frequencies while identifying the model. This translates to poor model predictions with respect to manipulated inputs and can result in poor controller performance. Another important contribution to the area of adaptive MPC has been by Nikolaou and co-workers.10-12 The simultaneous MPC and identification (MPCI) approach proposed by Genceli and Nikoloau10 solves an on-line optimization problem, which includes additional constraints that ensures persistent excitation (PE) necessary for on-line parameter estimation while minimally disturbing the process. The resulting NLP is computationally demanding, and to alleviate these difficulties Vutanadam and Nikolaou12 reformulated the PE constraints in the frequency domain. Genceli and Nikoloau10 have proposed to identify a model with finite impulse response (FIR) structure. Maiti and Saraf13 identified a subset of impulse response coefficients online and subsequently used them in MPC formulation. From the viewpoint of time-series modeling, the FIR representation has output error (OE) structure, and, as a consequence, these models generate excellent long-range predictions with respect to manipulated inputs. However, one drawback of the FIR representation is that the number of model parameters to be estimated is large. Thus, particularly for multivariable systems, the problem of estimating FIR model coefficients on-line can become highly ill conditioned. Also, the FIR model cannot capture the dynamics of the unmeasured disturbances affecting the plant. Thus, when an FIR model is used to formulate an MPC scheme, the effects of fast-changing unmeasured disturbances in the past on the future plant behavior cannot be predicted systematically. Vutanadam and Nikolaou12 and Shouche et al.11 formulate MPCI based on the DARX model, which is better suited for capturing the effects of stationary as well as nonstationary disturbances. Oshima et al.6 have proposed a novel method for on-line modeling of plant-model mismatch and unmeasured disturbances. By this approach, the residual signal (or prediction errors) generated by comparing the model predictions (developed once in the beginning) with the measured output are modeled as an ARX process driven by the past prediction errors and the manipulated inputs. The parameters of the ARX model are identified on-line using recursive leastsquare estimation. Whereas this approach facilitates on-line updation of deterministic as well as stochastic components, the use of FIR-type representation for the deterministic component and AR representation for the noise model again requires a large number of model parameters to be estimated on-line. In the certainty equivalence adaptive-control approach, the controller is designed under the assumption that the identified linear model gives an exact representation of the process dynamics in the neighborhood of the current operating point.14,15 In the review of the adaptive control presented at CPC-V, Ydstie15 identifies the admissibility problem as one of the key issues that must be addressed in certainty equivalence control to achieve bounded input bounded output (BIBO) stability and robustness with respect to unmodelled dynamics and disturbances. This implies that the estimated model must be well

behaved (controllable or stabilizable). The use of a Laguerre filter-based model is advocated to deal with this problem because these models are stable and robustly stabilizable, provided that the unmodelled dynamics are stable. One drawback of Laguerre or any other orthonormal basis filter-based model is that these models have fixed poles.16 As a consequence, the parameter adaptation is restricted to adaptation of system zeros and steadystate gains. When a process is operated over a wide operating range, time constants can change significantly and adaptive control based on fixed-pole models may not provide uniformly satisfactory performance over the entire operating range. Ydstie15 identifies the instability of the parameter estimator or the parameter drift as another important issue that must be addressed while developing an adaptive-control scheme. This drift can be avoided either by adding a deliberate perturbation to the set point or manipulated inputs or by using parameter projection to ensure that parameters do not go outside a bounded region.15 In this work, we develop an adaptive MPC formulation for multivariable time-varying systems based on a two-step modeling scheme. The key issue while developing an adaptive MPC scheme is choosing a model structure that can generate good long-range predictions and is capable of capturing the dynamics of stationary as well as nonstationary unmeasured disturbances. To incorporate excellent long-range prediction capability with respect to manipulated inputs, we propose to identify a parametrized form of an OE model, which are parsimonious in parameters, from the input-output data. The parameters of the OE component of the model are updated on-line by the recursive output error (ROE) method.17,18 An additional advantage of using ROE for parameter estimation is that it can facilitate control of multirate systems without requiring any significant modifications.20 However, the prediction error approach for parameter estimation requires that the one-step predictor be stable. As a consequence, the use of OE structure restricts the applicability of our approach to open-loop stable processes unless modifications are made in the predictor formulation for accommodating unstable poles, as suggested by Forselle and Ljung.19 To extract the model for the unmeasured disturbances, the residuals generated by the ROE estimators are further modeled as ARMA processes. The deterministic and stochastic components of the model are then combined to form a linear time-varying state-space model, which is further used to formulate the adaptive MPC (AMPC) problem at each sampling instant. Two novel features of the proposed AMPC formulation are highlighted below. (1) To deal with the problem of parameter drift that can cause excursion of system poles outside unit circle, we propose to use a recursive constrained estimation algorithm developed by Singh et al.,25 which ensures that the poles of the identified model are restricted to the unit circle. This also implicitly addresses that the admissibility issue as the resulting model is always open-loop stable and stabilizable. (2) We propose a two-tier modeling strategy that facilitates separate estimation of nonstationary and stationary components of the unmeasured disturbances. (3) In the objective function of AMPC, we include a terminalstate weighting term with a time-dependent weighting matrix, which is similar to the infinite-horizon formulation proposed by Muske and Rawlings.21 The use of the proposed recursive constrained estimation scheme ensures that the formulation proposed by Muske and Rawlings21 for open-loop stable systems can be used at every sampling instant. Thus, the terminal-state weighting matrix is recomputed at each time instant by solving

2710

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

a discrete Liapunov equation formulated using the identified state-space model. The applicability of this method to the control of semibatch processes and for transition control of large-scale industrial processes is illustrated by conducting simulation studies on (a) a fed-batch bioreactor process model involving multivariable control of dissolved oxygen and pH (b) control of large transitions in the benchmark Tennessee Eastman (TE) problem. The efficacy of the proposed adaptive MPC scheme is also established by conducting experimental studies on a benchmark quadruple-tank system.22 The rest of this article is organized as follows. The proposed model and the methodology used for parameter identification is described in section 2 and section 3, respectively. The formulation of a state-space model and MPC scheme is presented in section 4. The results of the simulation and experimental investigations are discussed in sections 5 and 6, respectively. The main conclusions reached based on the analysis of these results are presented in section 7. 2. MPC Relevant Model Development 2.1. Choice of Model Structure and Model Predictions. The primary concern in any MPC formulation is the quality of model prediction generated by the internal model. The quality of model predictions depends on (a) the fidelity of the deterministic component of the model and (b) the strategy used to model unmeasured disturbances and to predict their effects on future plant behavior. The quality of model predictions is linked with the model structure employed while identifying the MPC relevant model. Consider a general SISO process of the form

δy(t) ) G(q-1)δu(t) + V(t)

(1)

where δy ∈ R and δu ∈ R represent perturbations in measured output and manipulated input, respectively, G(q-1) represents the true process transfer function, and V(t) represents the (zeromean) unmeasured disturbance affecting the output. As the disturbances affecting chemical plants are often of drifting nature, a popular approach to cope with such disturbances is to employ an ARIMAX-type model5,7,8,15,9 of the form

δy(t) ) G ˆ (q-1)δu(t) +

H ˆ (q-1) e(t) ∆

(2)

where ∆ ) 1 - q-1. Here, G ˆ (q-1) represents the proposed model (typically of lower order) for capturing dynamics with respect to known input, and ∆V(t) ) H ˆ (q-1)e(t) represents the model of the unmeasured disturbances. For on-line identification, the above model has to be rearranged as follows

ˆ (q-1)e(t) ∆y(t) ) G ˆ (q-1)∆u(t) + H

∆y(t) ) y(t) - y(t - 1) and ∆u(t) ) u(t) - u(t - 1) (4) As the model is identified by minimizing the one-step prediction error, using Parseval’s theorem, the expression for asymptotic bias error in the estimation of model parameters θ can be expressed as follows18,23

θ

δyˆ (t) )

B(q-1) A(q-1)

q-tdδu(t)

(6)

δy(t) ) δyˆ (t) + V(t)

(7)

where A(q-1) and B(q-1) are polynomials in the shift operator q-1 and td represents time delay. In this work, we propose to use the above parametrized form of an OE model for capturing dynamics of the deterministic component of the system under the following assumptions: • Process under consideration is open-loop stable and can be modeled using an OE model with all of the poles inside the unit. In particular, if a system has integrating modes, then it is possible to approximate the behavior of these modes using stable poles that are arbitrarily close to the unit circle. • Time-delay matrix is time invariant. In the context of a continuously operated system, δy(t), δyˆ (t), and δu(t) represent perturbation variables defined with respect to some steady state, say (ys,us), as follows

δy(t) ) y(t) - ys; δyˆ (t) ) yˆ (t) - ys ; δu(t) ) u(t) - us We can rearrange the above OE model as follows

A(q-1)yˆ (t) ) B(q-1)u(t - td) + η

(8)

y(t) ) yˆ (t) + V(t)

(9)

where η represents a bias term defined as

(3)

where

θ* ) min

where Φu(ω) and Φ∆V(ω) represent the spectrum of known inputs signal, u(t), and disturbance signal, ∆V(t), respectively. As is evident from the above expression, the distribution of bias error |G(ejω) - G ˆ (ejω)| is influenced by the spectrum of the high-pass prefilter |∆(ejω)|2. Filtering of inputs and outputs through the high-pass operator ∆ pushes the model fit to the high-frequency range.18 As a consequence, the identified G(q-1) can have a significant bias at low frequencies, and this can translate to poor predictions when used in the MPC formulation. Thus, it is necessary to develop a modeling strategy that can capture drifting nature of unmeasured disturbances without requiring a compromise of the quality of model predictions at low frequencies. 2.2. Models for Deterministic and Stochastic Components. The most widely used FIR model in MPC formulations has OE structure. Because OE models are functions of past manipulated inputs alone, these models are ideally suited for long-term predictions. The FIR model, however, is nonparsimonious in parameters and, therefore, not suitable for on-line parameter estimation. A better choice for on-line parameter estimation is the parametrized form of an OE model, which can be represented as follows

∫-ππ |Hˆ (ejω1 ,θ)|2 [|G(ejω) G ˆ (ejω,θ)|2|∆(ejω)|2Φu(ω) + Φ∆V(ω)]dω (5)

η ) [A(1)ys - B(1)us]

(10)

The advantage of the above form is that we can work with absolute values of variables {y(t),u(t)} without requiring explicit knowledge of the steady-state operating point (ys,us).18 If it is desired to use this model for capturing dynamics of a timevarying system, such as a (a) batch or semibatch process or a (b) continuously processing plant operated over a wide operating range, then it becomes difficult to define perturbations with respect to a single steady state or to use a constant bias term η. In the context of such time-varying systems, we can modify the OE model (6-7) by introducing perturbations in the neighborhood of a time-varying trajectories {ys(t),us(t)} as


follows

yˆ (t) - ys(t) )

B(q-1,t) A(q-1,t)

A(i)(q-1,t)yˆ i(t) ) -d

q [u(t) - us(t)]

y(t) ) yˆ (t) + V(t)

(11) (12)

This model can be rearranged as

yˆ (t) )

B(q-1,t) A(q-1,t)

q-du(t) +

1 η(t) A(q-1,t)

y(t) ) yˆ (t) + V(t)

(13) (14)

m

∑ B(ij)(q-1,t)uj(t) + ηi(t) j)1

yi(t) ) yˆ i(t) + vi(t)

D(i)(q-1,t)vi(t) ) C(i)(q-1,t)ei(t) (15)

represents a time-varying unknown bias term. Note that argument (t) has been added to each polynomial in (q-1) to emphasize the fact that their coefficients are changing as a function of time. For the purpose of parameter estimation, the above model can be rearranged as follows:

V(t) ) y(t) - [1 - A(q-1,t)]yˆ (t) - B(q-1,t)u(t - d) - η(t) This assumes that • parametrization of the deterministic component is such that there is no correlation between the residual sequence {υ(t)} and the manipulated inputs sequence {u(t)} • the nonstationary component of the unmeasured disturbances is captured through {η(t)}, and sequence {υ(t)} represents a stationary stochastic process Under these assumptions, we can model its variation as an ARMA process as follows,

υ(t) )

C(q-1) D(q-1)

D(i)(q-1,t)vi(t) )

B(q-1,t) A(q-1,t)

u(t) +

V(t) )

1 η(t) + V(t) A(q-1,t)

C(q-1,t) D(q-1,t)

r

∑ C(ij)(q-1,t)ej(t) j)1

(23)

where e ∈ Rr represents the vector of innovations. The resulting MIMO model can be expressed in the following generic form

A(q-1,t)yˆ (t) ) B(q-1,t)u(t) + η(t)

(24)

y(t) ) yˆ (t) + v(t)

(25)

-1

-1

D(q ,t)v(t) ) C(q ,t)e(t)

(26)

where

A(q-1,t) ) I + A1(t)q-1 + A2(t)q-2 + ... + Ana(t)q-na (27) B(q-1,t) ) B1(t)q-1 + B2(t)q-2 + ... + Bnb(t)q-nb (28) C(q-1,t) ) I + C1(t)q-1 + C2(t)q-2 ... + Cnc(t)q-nc

(29)

D(q-1,t) ) I + D1(t)q-1 + D2(t)q-2 + ... + Dnd(t)q-nd (30) Here, Ai, Bi, Ci, and Di represent matrices of dimensions (r × r), (r × m), (r × r), and (r × r) , respectively. 3. Constrained Recursive Parameter Estimation

e(t)

(17)

where {e(t)} is a white-noise sequence and C(q-1) and D(q-1) are monic polynomials. Thus, the final form of the proposed model for time-varying processes can be expressed as follows,

y(t) )

(22)

or through MISO models

y(t) ) [1 - A(q-1,t)]yˆ (t) + B(q-1,t)u(t - d) + η(t) + V(t) (16) The time-varying bias term {η(t)} in the OE model given by eq 13 can also be viewed as a slowly drifting or nonstationary component of the unmeasured disturbance affecting the output. We develop a model for the residual signal, {V(t)}, defined as follows

(21)

where y ∈ Rr represents the controlled-variable vector, u ∈ Rm represents the manipulated variable vector, v ∈ Rr represents the residual vector for the OE model, and η ∈ Rr represents the bias vector. Here, A(i)(q-1,t) and B(ij)(q-1,t) represents polynomial functions of (q-1) operating on ith output and jth input, respectively. The dynamics of the residual signals can be captured either through SISO models of the form

where

η(t) ) [A(q-1,t)ys(t) - B(q-1,t)q-dus(t)]

(20)

e(t)

(18)

(19)

where y(t) and u(t) represent absolute values of the measurement and manipulated input. Also, when d > 0, the first d coefficients of B(q-1,t) will be zero. Because the deterministic component of the model has output error structure, the above formulation can be extended to an r × m multivariable system by developing r MISO models of the form (for i ) 1,...,r)

In this section, we propose a two-step identification strategy for on-line modeling that can deal with nonstationary disturbances while retaining the ability to generate good predictions at low frequencies. We develop the identification scheme for a SISO model. Extension to the MISO case is obvious and, hence, not discussed separately. In the development that follows, it is assumed that the time-delay matrix for the multivariable process under consideration is known a priori from off-line analysis of the data. 3.1. Constrained Parameter Estimation. The coefficients of polynomial A(q-1,t), B(q-1,t), and time-varying bias term η(t) can be estimated using the prediction error method, which minimizes variance of residual signal υˆ (t) defined as

V(t) ) y(t) - φ(t)Tθd(t)

(31)

φ(t) ) [-yˆ (t - 1) ... -yˆ (t - na) u(t - td - 1) ... u(t - td - nb) 1]T (32) θd(t) ) [a1(t) ... ana(t) b1(t) ... bnb(t) η(t)]

(33)

2712


In the present work, we assume that the dynamics of the process under consideration can be adequately represented by a set of local linear models with poles strictly inside the unit circle. Thus, while performing on-line parameter estimation, it is necessary to ensure that local models do not become unstable due to variance errors caused by low signal-to-noise ratios. Following Ydstie,15 the problem of on-line parameter estimation can be posed as a constrained optimization problem over a moving time window [t - N,t] as follows, t

[R(i)(υˆ (i))2 + ∑ θ (t) i ) t - N

θˆ d(t) ) min d

(∆θd(t))TPd(t - N)-1(∆θd(t))] ∆θd(t) ) θd(t) - θˆ d(t - N) subject to

θd(t) ∈ Θ where Θ represents constraint set for parameter estimates, R(i) represents the discount factor, and Pd(t - N) g I > 0 represents a positive definite covariance matrix, which defines a tradeoff between emphasis placed on new measurements and the old estimate. The constraints on the estimated parameters can be derived using the Jury’s stability criteria or the modified Routh stability criteria (Seborg et al.27). It may be noted that the resulting formulation is qualitatively similar to the nonlinear dynamic data reconciliation (NDDR),26 which is employed online for improving the quality of the operating data. The main difficulty with this formulation is that it results in a highly nonlinear constrained optimization problem and solving it online in real time may not be feasible. 3.2. Unconstrained Pseudolinear Regression. In adaptive control, the recursive solution to the unconstrained version of the above optimization problem is popularly employed for online computations mainly due to its computational simplicity. Even the unconstrained problem is a nonlinear optimization problem, which can be solved on-line using a recursiVe pseudolinear regression (RPLS) method called recursiVe output error (ROE) as follows17,18

Pd(t) )

λd(t - 1) + φ(t + 1)TPd(t - 1)φ(t + 1)

(34)

1 [I - Kd(t)φ(t + 1)T]Pd(t - 1) (35) λd(t - 1)

Here, Pd(t - 1) represents the covariance matrix of model parameters and λd(t - 1) is an exponential forgetting factor. To avoid loss of sensitivity and parameter blow-up in a region where the covariance of the parameters becomes very small, the forgetting factor can be updated as suggested by Fortescue et al.,24 λd(t) ) 1 -

[1 - φ(t + 1)TKd(t)][y(t + 1) -φ(t + 1)Tθd(t)]2

∑

θd(t + 1)

(δθd(t + 1))TPd(t)-1(δθd(t + 1))] δθd(t + 1) ) θd(t + 1) - θˆ d(t)

θd(t + 1) ∈ Θ

where gain Kd(t) is computed as follows

Pd(t - 1)φ(t + 1)

θˆ d(t + 1) ) min [[υˆ (t + 1)]2 +

subject to

θˆ d(t + 1) ) θˆ d(t) + Kd(t)[y(t + 1) - φ(t + 1)Tθˆ d(t)]

Kd(t) )

number of significant recent samples to be considered for model building (N0d). The parameters of the ARMA model can be estimated using an extended recursive least-square (ELS) algorithm, which also belongs to the family of recursive pseudolinear regression (RPLR) methods.17 Convergence results for RPLS scheme for models with OE and ARMA structure have been discussed by Ljung.18 3.3. Constrained Recursive Estimation. The unconstrained recursive formulation requires much less computational efforts when compared to the moving window constrained parameter estimation. However, this strategy is susceptible to parameter drift and the predictor can become unstable due to variance errors. Thus, if it is desired to use on-line recursive parameter estimation, admissibility and drift issues need to be addressed to achieve BIBO stability.15 The moving window formulation suggested by Ydstie15 as a remedy to these problems has a very high computational cost. Recently, in analogy with the recursive nonlinear dynamic data reconciliation approach proposed by Vachhani et al.,26 Singh et al.25 have proposed a novel constrained recursive formulation for on-line parameter estimation. This approach combines the efficiency of on-line recursive estimation and the ability of dynamic data reconcilation techniques to handle algebraic inequality and equality constraints. In this work, we propose to adopt this approach to address the admissibility and drift issues. To understand how this scheme works, let us assume that parameter estimates θˆ d(t) ∈ Θ are available at instant (t + 1). We first compute θ˜ d(t + 1) using the unconstrained recursive estimation scheme discussed above. If estimated θ˜ d(t + 1) ∈ Θ (i.e., roots of A(q-1,t + 1) polynomials are inside the unit circle) then we set θˆ d(t + 1) ) θ˜ d(t + 1) and proceed with the control-law calculation. However, if θ˜ d(t + 1) ∉ Θ then we formulate a constrained optimization problem over one sampling period as follows

(36)

0d

where, ∑0d is the product of expected measurement noise variance based on real knowledge of the process (σ0d2) and the

The well-known Jury’s stability criterion is used to impose constraints in the parameter space. In particular, when polynomial A(q-1) is of first or second order, Singh et al.25 have shown that the above constrained optimization problem can be reduced to a QP problem. After solving the constrained optimization problem, updation of the covariance matrix is carried out using eqs 35 and 36 as in recursive estimation. This formulation ensures that the poles of A(q-1) remain within the unit circle at significantly less computational cost when compared to the moving window formulation. Though the proposed constrained optimization approach ensures the stability of the estimated model, the procedure is computationally cumbersome when compared to the recursive estimation techniques. Thus, we propose to combine the two approaches to make the parameter estimation scheme computationally efficient. A check for stability is performed every time the parameters are estimated using the unconstrained recursive estimation procedures outlined above. We define the spectral radius of A(q-1) as


F ) max |pi|

a state realization for the stochastic model can be expressed as follows9

i

where pi are zeros of A(q-1). The constrained optimization procedure is invoked only when F exceeds unity. While developing an ARMA model using the ELS algorithm, it is required that the poles of polynomial D(q-1) remain within the unit circle to meet the requirement that the predictor be stable. If poles of D(q-1) are found to move outside the unit circle, then a constrained optimization problem can be formulated and solved in a similar manner. 4. Adaptive Predictive Control Formulation 4.1. MPC Relevant State Space Model. While developing an MPC scheme based on the above model, it is convenient to work with a state realization of the identified model (24-26). The deterministic component of the proposed model can be expressed as

yˆ (t) ) -A1(t)yˆ (t - 1) - ... -Ana(t)yˆ (t - na) + B1(t)u(t - 1) + ... + Bnb(t)u(t - nb) + η(t) (37) y(t) ) yˆ (t) + v(t) Thus, defining state vector

xd(t)

at time instant t as

T

T

T

u(t - nb + 1)T]T (39)

(48)

v(t) ) Hsxs(t)

(49)

where w(t - 1) ) e(t) and

Φs(t) ) -D1(t) I : [0] [0] [0] : [0]

-D2(t) ... -Dnd(t) C1(t) [0] [0] ... [0] : [0] ... I [0] [0] [0] ... [0] [0] ... [0] I : [0] [0] ... [0]

where

yˆ (t) ) Hdxd(t)

[

-A1(t) I : Φd(t) ) [0] [0] [0] :

-A2(t) [0] : ... [0] [0] :

... ... : I ... ... :

-And(t) B2(t) [0] [0] : : [0] [0] [0] [0] [0] I : :

.... [0] : [0] [0] [0] :

Bnc(t) ... : ... ... ... :

]

(41)

(42)

Ψη ) [I [0] ... [0] ... [0]]T

(44)

Hd ) [I [0] ... [0]]

(45)

v(t) ) -D1(t)v(t - 1) - ... -Dnd(t)v(t - nd) + e(t) + C1(t)e(t - 1) + ... + Cnc(t)e(t - nc) (46) Similarly, defining the state vector xs(t) as

xs(t) ) [v(t) ... v(t - nd + 1) e(t) ... e(t - nc + 1) ]

T T

I [0]

]

(50)

(51) (52)

Note that the system matrices (Φ,Γ) appearing in the state realizations constructed above differ due to structural differences between the OE and the ARMA models. We define an augmented state vector as follows:

x(t) )

[ ] xd(t) xs(t)

A combined state-space model can be obtained at tth sampling instant by augmenting the two state realizations,

(54)

where

(43)

T

...

y(t) ) Hx(t)

Γd ) [BT1 (t) [0] ... [0] I [0] ... [0]]T

T

... [0] ... [0] ... [0]

x(t) ) Φ(t)x(t - 1) + Γ(t)u(t - 1) + Ψηˆ (t) + Kw(t - 1) (53)

Here, [0] and I represent null matrices and identity matrices of appropriate dimensions. The stochastic component of the model can be represented as

T

[0] [0] [0]

Hs ) [I [0] ... [0]]

the deterministic process model in discrete state-space form can be given by

xd(t) ) Φd(t)xd(t - 1) + Γd(t)u(t - 1) + Ψηηˆ (t) (40)

C2(t) ... Cnc(t) [0] ... [0]

Γs ) [I [0] ... [0] I [0] ... [0]]T

(38)

x (t) ) [yˆ (t) ... yˆ (t - na + 1) u(t - 1) .... d

[

xs(t) ) Φs(t)xs(t - 1) + Γs(t)w(t - 1)

(47)

Φ(t) )

[

]

[ ] [ ] [ ]

Ψη Φd(t) [0] Γ (t) ; Γ(t) ) d Ψ ) ; Φs(t) [0] [0] [0] [0] K ) Γ ; H ) [Hd Hs ] (55) s

4.2. AMPC Formulation with Time-Varying Control Objective. We briefly describe the predictive control formulation based on the state-space model developed above. At each sampling instant, the linear time-varying state-space model given by eqs 53-54 is used for predicting future behavior of the plant over a finite future time horizon of length p (prediction horizon) starting from current instant t. To carry out predictions based on the above model, it becomes necessary to introduce a model for the behavior of the bias term. Because η(t) is expected to capture slowly varying disturbances, we assume that the bias term remains constant over the prediction horizon, that is,

ηˆ (t + j + 1|t) ) ηˆ (t + j|t)

(56)

ηˆ (t|t) ) ηˆ (t)

(57)

for j ) 0,1,....p - 1. It is further assumed that only q (control horizon) future manipulated input moves can be chosen freely with following input blocking constraints

2714


u(t + j|t) ) u(t|t) for j ) 0,1,...c1 - 1

(58)

u(t + j|t) ) u(t + c1|t) for j ) c1,...,c2 - 1

(59)

........ u(t + j|t) ) u(t + cq-1|t) for j ) cq-1,...p - 1

(60)

where cj are selected such that

0 < c1 < c2 < .... < cq-1 Given a desired final set point r(t), a filtered future set-point trajectory can be generated at instant t as follows

xr(t + j + 1|t) ) Φrxr(t + j|t) + Γr[r(t) - y(t)]

QL(t) ) (Hd)TWE(Hd) + [Φd(t)]TQL(t)[Φd(t)]

The inclusion of the terminal weighting as described above is similar to the infinite-horizon formulation of linear MPC proposed by Muske and Rawlings.21 When number of manipulated inputs exceeds the number of controlled outputs, timevarying target state xj(t) appearing in terminal-state error (t + N|t) can be computed by solving the following optimization problem

min us(t)TWUus(t)

(72)

[I - Φ(t)]xj(t) ) Γ(t)us(t) + Γηηˆ (t)

(73)

r(t) ) H(t)xj(t)

(74)

us(t)

subject to

(3.16)

yr(t + j|t) ) y(t) + Crxr(t + j|t)

(61)

L

for

u e us(t) e u j ) 0,1,...p - 1

with initial condition xr(t|t) ) 0. Here, the set-point filter matrices are tuning parameters and are selected such that the steady-state gain of the filter Cr(I - Φr)-1Γr ) I. The adaptive MPC problem at the sampling instant t is formulated as a constrained optimization problem as follows

{

(t + p|t) QL(t)(t + p|t) +

min Uf(t)

T

p-1

ef(t + j|t) WEef(t + j|t) ∑ j)1 T

cq-1

+

∑ j)1

∆u(t + j|t)TW∆U ∆u(t + j|t)

}

(62)

(t + p|t) ) x(t + p|t) - xj(t)

(64)

xˆ (t + j|t) ) Φ(t)xˆ (t + j - 1|t) + Γ(t)u(t + j - 1|t) + Ψ(t) (65) yˆ (t + j|t) ) Hxˆ (t + j|t)

(66)

∆u(t + l|t) ) u(t + l|t) - u(t + j - 1|t)

(67)

y e yˆ (t + j|t) e y

(68)

∆u e ∆uf(t + j|t) e ∆u L

(69) H

min [r(t) - yj(t)]TWE[r(t) - yj(t)]

(75)

[I - Φ(t)]xj(t) ) Γ(t)us(t) + Γηηˆ (t)

(76)

yj(t) ) H(t)xj(t)

(77)

us(t)

u e us(t) e u

(63)

uL e u(t + j|t) e uH

If the number of controlled outputs exceeds or is equal to the number of manipulated inputs, the target state xjd(t) can be computed by solving the following optimization problem

L

ef(t + j|t) ) yr(t + j|t) - yˆ (t + j|t)

H

H

subject to

subject to following constraints

L

(71)

H

The above optimization formulations can be transformed into QP problems for improving on-line computational efficiency. The controller is implemented in a moving horizon frame work. Thus, after solving the QP problem, only the first move uopt(t|t) is implemented on the plant, that is,

u(t) ) uopt(t|t)

(78)

and the optimization problem is reformulated at the next sampling instant based on the updated information from the plant and updated model. It may be noted that the AMPC formulation with the terminal weighting proposed above can be employed only for controlling a continuously operated process. The infinite-horizon formulation (i.e., with terminal-state weighting) does not hold much significance in the context of batch or semibatch processes, which are operated over a finite time interval.

(70)

j ) 1,2,....p together with the input blocking constraints given by eqs 5860. Here, WE and W∆U represent the error weighting matrix and the input move weighting matrix, respectively. It may be noted that, unlike conventional MPC formulations, the objective function of the proposed adaptive MPC formulation has a timevarying weighting matrix QL(t) appearing in the terminal weighting term. Because the identified model has all of its poles inside the unit circle, assuming that the certainty equivalence principle holds, the time-varying weighting matrix QL(t) in the terminal-state weighting is computed at each sampling instant by solving the following discrete Liapunov equation:

5. Simulation Case Studies The applicability of the proposed adaptive-predictive control scheme to the control of the semibatch process, and the transition control of large-scale industrial processes is illustrated by simulation studies on following systems: • Multi-variable control of dissolved oxygen and pH in fedbatch bioreactor. • Control of large transitions in Tennessee Eastman (TE) problem. The first problem involves maintaining pH and dissolved oxygen (DO) in a bioreactor at certain desired levels throughout the entire operation, whereas the second problem involves large transitions from one operating conditions to another. Note that


both of the control problems require continuous and significant changes to be made to the manipulated inputs to meet their respective goals; and this, in turn, ensures the persistency of excitation. In both of the cases, the parameters of the MISO OE models are identified using a combination of the unconstrained and the proposed constrained recursive parameter estimation strategy. In addition, we implement a fixed deadzone15 and turn off parameter estimation when the norm of the modeling error is lower than a threshold. The standard deviation of measurement noise in each output is used as a basis for selecting the threshold for the corresponding MISO parameter estimator. 5.1. Penicillin Production in a Fed-Batch Bioreactor. An unstructured model for penicillin production in a fed-batch fermentor given by Birol et al.28 has been used for simulation studies. The model consists of nine differential equations relating fermenter states (biomass, substrate, product, dissolved oxygen, CO2 and hydrogen ion concentrations, fermenter temperature, and fermenter volume) with six inputs (substrate feed and coolant flows, acid and base flows, agitation rate, and aeration). The model equations, initial conditions, and nominal values of the model parameters are given in Birol et al.28 While carrying out the simulations, we assume that the fermenter operates isothermally and that the fermentor temperature is maintained at the optimum value (25 °C) by a perfect temperature controller. We also assume that the substrate feed rate is fixed at a constant optimum value (0.042 liter/h), which reproduces a predecided growth pattern. The control objective in this case study is to maintain the dissolved oxygen concentration and pH of the fermentation medium at the set-point by manipulating aeration, agitation, and acid/base flows. As the culture makes transitions from one growth phase to the other (lag phase f exponential phase f stationary phase), the oxygen uptake rate and acid secretion rates change substantially. With an increase in the batch time, the biomass accumulation rate increases and this results in higher oxygen consumption. Similarly, as time progresses, the acid secretion rates increases. Thus, the manipulated inputs have to be changed continuously to keep pH and DO at desired set points. As the input-output dynamics changes significantly with time, control of pH and DO at their respective set points is an interesting control problem. Simulations are carried out for a batch time of 100 h, with a sampling time of 0.5 min. Measurement noise is introduced in pH (σpH ) 0.01) and DO (σDO ) 0.02) by adding a zero-mean noise sequence with Gaussian distribution. To initialize the parameter estimator, perturbations are introduced in the manipulated variables under open-loop conditions for the initial 310 sampling instants. Two MISO OE models and two SISO ARMA models of the form -1 -2 [1 + a(i) + a(i) ˆ i(t) ) 1 q 2 q ]y 4

-1 -2 [1 + b(ij) + b(ij) ∑ 1 q 2 q ]uj(t) + ηi(t) j)1

(79)

yi(t) ) yˆ i(t) + vi(t)

(80)

-1 (i) -1 [1 + d(i) 1 q ]vi(t) ) [1 + c1 q ]ei(t)

(81)

for i ) 1,2 are identified on-line and used to implement AMPC scheme (without terminal weighting) starting from the 311th sampling instant. The constraints on the inputs and input rates are given in Table 1, and the tuning parameters used in AMPC formulation are listed in Table 2.

Figure 1. Bioreactor control: Variation of controlled outputs.

Figure 2. Bioreactor control: Varition of manipulated inputs. Table 1. Bioreactor Control: Input and Input-Rate Constraints manipulated variable

input constraints

input rate constraints

aeration rate (1/h) agitation power (W) acid flow rate (1/h) base flow rate (1/h)

0 to 30 2 to 30 0 to 8 × 10-5 0 to 3 × 10-3

-0.1 to 0.1 -1 to 1 -3 × 10-7 to 6 × 10-7 -4 × 10-6 to 3 × 10-6

Table 2. Bioreactor Control: AMPC Tuning Parameters tuning parameter

value

prediction horizon control horizon input blocking penalty on error in DO penalty on error in pH penalty on aeration penalty on agitation penalty on acid flow penalty on base flow

50 4 [1 1 1 47] 1 20 1 1 1 1

The regulatory control response with the AMPC strategy for the given set points in %DO (DOset ) 90%) and pH (pHset ) 5) are plotted in Figure 1. The corresponding profiles of biomass, substrate, and penicillin concentrations generated using AMPC are shown in Figure 3. It is evident from these plots that even at the time of large state transitions, the variations in the controlled variables (%DO and pH) from the set point are negligible. This is achieved by continuously changing the manipulated inputs in accordance to the changing dynamics in oxygen uptake rate and acid secretion rate (Figure 2). Figures 4 and 5 show the variations of spectral radii (magnitude of maximum magnitude root) of A(i)(q-1) with respect to batch time for DO and pH models, respectively. It may be noted that the constrained recursive estimator has to be

2716

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 Table 3. TE Process: Local PI Loops Tuning Parameters control variable reactor level (Lr) stripper level (Lstr) separator level (Lsep) reactor pressure (Pr)

Figure 3. Bioreactor control: Variation of biomass, substrate, and product concentration.

Figure 4. Bioreactor control: Variation of spectral radius of a polynomial in the MISO model for DO.

Figure 5. Bioreactor control: Variation of spectral radius of a polynomial in the MISO model for pH.

invoked at several sampling instants in both the cases to ensure that the roots of A(i)(q-1) are maintained within the unit circle. Typical scenarios where the unconstrained recursive estimator causes excursions of the spectral radii outside the unit circle and subsequent projection of the parameters using constrained recursive formulation to maintain the spectral radius within the unit circle are also highlighted in these figures. As is evident from these figures, the proposed constrained recursive formulation successfully maintains the poles of the identified model inside the unit circle. The significant variations of spectral radii with time also reflect the fact that the DO and pH dynamics

manipulated variable E feed (F3) stripper bottoms outlet, (F11) separator bottoms outlet, (F10) purge, (F9)

Kc

τI (hr)

2.0162 -3.156

3.63 0.1

-2.975

0.13

-10.768

0.001

are changing continuously in the exponential growth phase as well as in the stationary phase. 5.2. Tennessee Eastman (TE) Process. To demonstrate the ability of the proposed AMPC algorithm to manage large transitions in the operating point in a large-scale system, we present a case study of the benchmark TE problem. The TE plant-wide control problem, developed by Downs and Vogel,29 has been a useful platform for testing strategies for plant-wide control, multivariable control, on-line optimization, predictive control and system identification. The process is highly nonlinear and is open-loop unstable with a large number of measured and manipulated variables. The plant produces two products (G and H) from four reactants (A, C, D, and E). Two other byproducts are present in the form of an inert (B) and a byproduct (F). The process has five major unit operations: reactor, product condenser, vapor/liquid separator, recycle compressor, and product stripper. The feed to the process is in the gaseous phase. The gaseous products from the reactor pass through a condenser and subsequently the vapor/liquid separator, where the condensed products are separated from the unconverted reactants. This stream of unconverted reactants then passes through the recycle compressor before re-entering the reactor. The desired products (G and H) are obtained from the base of the stripper column, whereas the inert (B) and the byproduct (F) are purged as vapors from the top of the vapor-liquid separator. Ricker and Lee30 have developed a reduced-order nonlinear, mechanistic model consisting of 26 states and 10 manipulated inputs for the TE plant. We use this mechanistic model together with the base case bias parameters reported in Ricker and Lee30 to simulate the TE plant dynamics in our simulation studies. Note that the measurements are corrupted with zero mean and normally distributed measurement noise. Details of noise model parameters used in the simulation studies can be found in Ricker and Lee.30 To control the TE process, it is required to stabilize the reactor and satisfy the product specifications. The process parameters like reactor level, stripper level, separator level, and reactor pressure exhibit faster dynamics and are therefore controlled at a faster rate (at 3 s) in our simulation studies. To make the process open-loop stable, four local PI control loops are employed on these controlled variables. The specifications of these loops are summarized in Table 3. The important or primary-controlled variables for AMPC are product concentration and production rate. Some of the critical intermediate stream concentrations (A in feed, E in feed, and B in purge stream) are also required to be within an acceptable range. During transitions, it is important to have a strict control on reactor pressure as it is quite sensitive to changes in other process variables. While implementing the AMPC scheme, we assume that PI loops on the reactor level, stripper level, separator level, and reactor pressure are closed and that their set points are available for manipulation to the master AMPC in a cascade control configuration. In addition to these set points, flow rates of feed 1, feed 2, feed 4, feed 8, reactor temperature and separator temperature are also used as manipulated inputs by AMPC. These are tabulated in Table 4 along with the penalties


Figure 6. Transition control in the TE plant: Primary controlled outputs.

Figure 7. Transition control in the TE plant: Secondary controlled outputs (continuous line) and setpoints (dashed line).

Table 4. TE Process: Weights on Error and Input-Rate in AMPC Objective Function controlled variable

penalty on error

% G in product

5

production rate, F11

5

% B in purge % A in feed % E in feed reactor pressure

10 2 5 5

manipulated input

penalty on input moves

reactor level set point stripper level set point separator level set point reactor pressure set point F1 F2 F4 F8 reactor temperature separator temperature

1 1 1 0.2 1 1 1 1 1 1

Table 5. TE Process: Input and Input-Rate Constraints in AMPC manipulated variable

input constraints

input rate constraints

reactor level set point (%) stripper level set point (%) separator level set point (%) reactor pressure set point (kPa) F1 (kmol/h) F2 (kmol/h) F4 (kmol/h) F8 (kmol/h) reactor temperature (°C) separator temperature (°C)

10 to 90 10 to 90 10 to 90 2650 to 2850 0.1 to 45 1 to 181 1 to 681 850 to 1460 115 to 128 75 to 88

-5 to 2 -0.2 to 0.2 -0.2 to 0.2 -0.2 to 0.2 -0.1 to 0.1 -1 to 1 -2 to 2 -2 to 2 -0.01 to 0.01 -0.01 to 0.01

used in the AMPC objective function (i.e., WE and W∆U matrices). Table 5 gives the nominal constraints on manipulated variables. The supervisory control action of AMPC is performed at every 0.1 hour interval. It is assumed that all of the concentration measurements are available at 0.1 hour intervals. The prediction horizon is chosen as p ) 80 (7.5 hrs) and control horizon is chosen as q ) 5 with following input blocking constraints

Input Blocks ) [5 10 15 20 25 ]

(82)

For deterministic modeling, six MISO OE and SISO ARMA models of the form 79-81 are recursively identified and then used to formulate AMPC formulation (with terminal-state weighting). The control objective of this case study is to achieve a grade transition starting from the base case (i.e., 50% G in product) to 10/90 case (i.e., 10% G in product) and then bring the process back to base case. Figures 6 and 7 present the profile of the controlled variables for a grade transition from the base case to 10/90 and back again to the base case using AMPC with terminal weighting. Because

Figure 8. Transition control in the TE plant: Manipuated inputs.

Figure 9. Transition control in the TE plant: Manipulated inputs.

the terminal weighting matrix QL(t) is recomputed each time the parameters are updated, time varying QL(t) is reported in terms of the ratio |QL(t)|2/|WE|2 as a function of time (Figure 10). It may be noted that the set points have to be changed gradually when the terminal-state weighting is used in the AMPC formulation. If the current set point specified (i.e., r(k)) is far off from the current operating point, then the terminalstate xj(t) predicted using the local linear model can be significantly different from the terminal state at the final steady state. This difficulty can be alleviated by introducing a set-point filter. It can be seen that AMPC is able to achieve smooth transition of the important controlled variables like product concentration, production rate, and reactor pressure to the new set points. Also, these variables are tightly controlled at their

2718


Figure 10. Transition control in the TE plant: Variation of the condition number of the gain matrix and |WL|/|WE| as functions of time.

Figure 12. Quadruple tank experimental setup: Schematic diagram.

Figure 11. Transition control in the TE plant: Spectral radii of A(q-1) for MISO OE models.

respective set points once the set points are achieved. The transition of the remaining controlled variables (%B in purge, %A in feed, and % E in feed) takes a relatively longer duration and is achieved by maintaining them well within the shut-down limits. The corresponding manipulated variable profiles are plotted in Figures 8 and 9. Figure 10 reports variation of the condition number of the steady-state gain matrix K(t) defined as K(t) ) H[I - Φ(t)]-1Γ(t). Figure 11 depicts the spectral radii of MISO OE models as a function of time. These remain inside the unit circle at all the sampling instants. It is evident from Figures 10 and 11 that the deterministic model parameters undergo significant variations during grade transitions. 6. Experimental Studies The efficacy of the proposed AMPC scheme is also demonstrated by conducting experimental studies on the benchmark quadruple-tank system. The quadruple-tank system used for the experimental verification is a scaled and modified version of the benchmark setup presented by Johanson.22 The schematic of this process is shown in Figure 12, and details of the experimental setup are presented in Table 6. The controlled variables are the levels of tanks 1 and 2, and the manipulated variables are control-valve positions. The system exhibits interacting multivariable dynamics with each of the valves affecting both of the outputs. It has an adjustable multivariable zero that can be set to a right-half or left-half plane value by changing the valve settings. In addition to this, in our setup spherical bulbs are placed in the upper-two tanks to introduce additional nonlinearity in the dynamics (Figures 12 and 13).

Figure 13. Quadruple tank experimental setup: Operating range in tank 3 and tank 4.

Table 6. Quadruple Tank Setup: Specifications S. No.

variable

value

1 2 3 4 5 6 7 8 9 10 11 12 13

inner diameter of tanks height of bottom tank over flow height of bottom tank over flow height of top tank orifice diameter of top tanks orifice diameter of bottom tanks sphere diameter cylindrical connection between spears clearance between sphere and top tank bottom pipe line inner diameter pipe line outer diameter pump rating the control-valve coefficient

15 cm 60 cm 55 cm 43 cm 8.5 mm 10.5 mm 13 cm 4 cm 4.5 cm 15 mm 21.6 mm 0.25 Hp 4.7

The water in the upper-two tanks is allowed to flow in the annular space between the tank and the spherical bulbs. As the height of water in these two tanks changes, the cross section for the flow also changes, resulting in changing time constants and time-varying behavior. The flow ratios (0.43, 0.34) are adjusted in such a way that maximum flow is diverted to the upper tanks. In such a situation, this process is characterized by a right-half plane zero and exhibits non-minimum phase behavior. The control objective in this case study is to compare the tracking ability of the proposed AMPC scheme with that of the conventional MPC (CMPC), which is based on a fixed OE


Figure 14. Quadruple tank setup: Comparison of controlled outputs.

Figure 15. Quadruple tank setup: Manipulated input profiles. Table 7. Quadruple Tank System: DMC/AMPC Tuning Parameters prediction horizon control horizon penalty on level 1 penalty on level 2 penalty on input 1 penalty on input 2 set-point filter

60 1 12 10 1 1 0.3

Table 8. Quadruple Tank System: Model Parameters parameter a(11) 1 a(11) 2 b(11) 1 b(11) 2 b(12) 1 b(12) 2 η1 d(1) 1 c(1) 1

(min, max) value -1.7, -1.37 4.12 × 10-1, 7.42 × 10-1 1.63 × 10-2, 6.19 × 10-2 -3.45 × 10-2, -1.95 × 10-2 -2.83 × 10-2, -1.7 × 10-2 7.71 × 10-2, 9.94 × 10-2 1.57 × 10-2, 3.43 × 10-2 -8.63 × 10-1, -8.35 × 10-1 2.17 × 10-1, 2.35 × 10-1

parameter a(21) 1 a(21) 2 b(21) 1 b(21) 2 b(22) 1 b(22) 2 η2 d(2) 1 c(2) 2

(min, max) value -1.7, -1.5 5.06 × 10-1, 7.31 × 10-1 -6.7 × 10-2, -3.53 × 10-2 1.44 × 10-2, 1.94 × 10-2 9.12 × 10-2, 1.15 × 10-1 -5.28 × 10-3, 9.31 × 10-3 2.15 × 10-2, 3.36 × 10-2 -8.53 × 10-1, -8.33 × 10-1 3.94 × 10-1, 4.54 × 10-1

model developed at the operating point, for large changes in the set points. To make a fair comparison, identical tuning parameters have been used for CMPC and AMPC controllers (Table 7). The operating point at which the model used in CMPC formulation is developed is also reported in Table 7. In AMPC

Figure 16. Quadruple tank setup: Variation of model parameters in AMPC.

formulation, two MISO OE models of the form (79-80) and two SISO ARMA models of the form (81) are recursively identified and then stacked to generate a MIMO model. Figures 14 and 15 compare the servo responses of CMPC and the AMPC strategy for large set-point changes. The variations of normalized OE and ARMA model parameters in AMPC formulation are reported in Figure 16, whereas the minimum and maximum values of the model parameters are listed in Table 8. Figure 13 shows the operating range in terms of level variations in tanks 3 and 4 when the set-point changes are introduced in the levels of tanks 1 and 2. As can be seen from Figure 14, the CMPC leads to large overshoots in level responses mainly due to time-varying dynamics in tanks 3 and 4. The proposed AMPC approach, on the other hand, is able to adapt to time-varying behavior of the system and achieve the desired transition without any overshooting. 7. Conclusions In this work, we develop a novel adaptive model predictive control (MPC) formulation for multivariable time-varying systems. A two-tier modeling scheme is proposed in which the deterministic and stochastic components of the model are updated on-line by two separate recursive pseudolinear regression schemes. To address the admissibility issue and the parameter-drift problem in on-line estimation, we propose to use a combination of a novel constrained recursive estimation scheme and the conventional recursive method for model parameter estimation. The residuals generated by the OE models are modeled as ARMA processes to account for unmeasured disturbances. An interesting feature of the proposed modeling scheme is that it allows separation of stationary and nonstationary components of the unmeasured disturbances. The deterministic and stochastic components of the model are combined to form a linear time-varying state-space model, which is further used to formulate the predictive control problem at each sampling instant. A novel feature of the proposed controller formulation is that we develop an infinite-horizon formulation with a time-varying objective function in the adaptive-control framework for dealing with grade-transition problems in continuously operated plants. The applicability of this method to the control of the semibatch and continuously operated largescale industrial processes is illustrated by carrying out simulation studies on a penicillin production process and on the benchmark Tennessee Eastman (TE) control problem. The proposed modeling and control scheme is also validated by conducting experiments on the benchmark quadruple-tank setup. The

2720


analysis of the simulation results reveals that the proposed twotier modeling scheme successfully captures system dynamics with respect to manipulated inputs as well as unmeasured disturbances over a wide operating range. The proposed adaptive MPC (AMPC) scheme is able to achieve tight control of timevarying semibatch processes and is capable of managing large transitions in the operating point of a continuously operated complex multivariable process. The experimental verification demonstrates that the AMPC scheme performs considerably better than the conventional nonadaptive MPC scheme for the servo control problems. Supporting Information Available: Complete ref 25 in PDF format. This material is available free of charge via the Internet at http://pubs.acs.org. Literature Cited (1) Morari, M.; Lee, J. H. Model Predictive Control: Past, Present and Future. Comput. Chem. Eng. 1999, 23, 667-682. (2) Qin, S. J.; Badgwell, T. A. A Survey of Industrial Model Predictive Control Technology. Control Engineering Practice 2003, 11, 733-764. (3) Berber, R. Control of Batch Reactors - A Review. Trans. Inst. Chem. Eng. 1996, 74 (Part A), 3-20. (4) Henson, M. A. Nonlinear Model Predictive Control : Current Status and Future Directions. Comput. Chem. Eng. 1998, 23, 187-202. (5) Clarke, D. W.; Mohtadi, C.; Tuffs, P. S. Generalized Predictive Control - I. The Basic Algorithm. Automatica 1987, 23, 137-148. (6) Oshima, M.; Ohno, H.; Hashimoto, I.; Sasajima, M.; Maejima, M.; Tsuto, K.; Ogawa, T. Model Predictive Control with Adaptive Disturbance Prediction and Its Application to Fatty Acid Distillation Column Control. J. Process Control 1995, 5 (1), 41-48. (7) Diaz, C.; Dieu, P.; Feuillerat, C.; Lelong, P.; Salamie, M. Adaptive Predictive Control of Dissolved Oxygen Concentration in a LaboratoryScale Bioreactor. J. Biotechnol. 1995, 43, 21-32. (8) Rodrigues, J. A. D.; Toledo, E. C. V.; Filho, R. M. A Tuned Approach of the Predictive Adaptive GPC Controller Applied to a FedBatch Bioreactor Using Complete Factorial Design. Comput. Chem. Eng. 2002, 26, 1493-1500. (9) Camacho, E. C.; Bourdons, C. 1999, Model PredictiVe Control; Springer-Verlag: London. (10) Genceli, H.; Nikolaou, M. New Approach to Constrained Predictive Control with Simultaneous Model Identification. AIChE J. 1996, 42, 10, 2857-2867. (11) Shouche, M.; Genceli, H.; Vuthandam P.; Nikolaou, M. Simultaneous Constrained Model Predictive Control and Identification of DARX Processes. Automatica 1997, 34, 12, 1521-1530. (12) Vuthandam, P.; Nikolaou, M. Constrained MPCI: A Weak Persistent Excitation Approach. AIChE J. 1997, 43, 9, 2279-2288.

(13) Maiti, S. N.; Saraf, D. N. Adaptive Dynamic Matrix Control of a Distillation Column with Closed-Loop Online Identification. J. Process Control 1995, 5 (5), 315-327. (14) Seborg, D. E.; Edgar, T. F.; Shah, S. L. Adaptive Control Strategies for Process Control: A Survey. AIChE J. 1986, 32 (6), 881-913. (15) Ydsti, B. E. Certainty Equivalence Adaptive Control: What’s New in the Gap. In Chemical Process Control - V; Kantor, J. C., Garcia, C. E., Carnaham, B., Eds.; CACHE-AIChE: Tahoe City, California, 1997, 9-23. (16) Patwardhan, S. C.; Shah, S. L. From Data to Diagnosis and Control Using Generalized Orthonormal Basis Filters. Part I: Development of State Observers. J. Process Control. 2005, 15, 7, 819-835. (17) Astrom, K. J.; Wittenmark, B. AdaptiVe Control, 2nd ed.; Pearson Education Asia, 2001. (18) Ljung, L. 1999, System Identification: Theory for the User, 2nd ed.; Prentice Hall, Inc.: Eaglewood Cliffs, New Jersey. (19) Forsselle, U.; Ljung, L. Identification of unstable systems using output error and Box-Jenkins model structures, Proceedings of the 37th IEEE CDC, Tampa, Florida, 1998; 3932-3937. (20) Lu, W.; Fisher, D. G. Least Square Output Estimation with MultiRate Sampling. IEEE Trans. Autom. Control 1989, 34, 6, 669-672. (21) Muske, K. R.; Rawlings, J. B. Model Predictive Control with Linear Models. AIChE J. 1993, 39, 262-287. (22) Johansson, K. H. The Quadruple-Tank Process: A Multivariable Laboratory Process with an Adjustable Zero. IEEE Trans. Control Systems Technology 2000, 8 (3), 456-465. (23) Jorgenson, S. B.; Lee, J. H. Recent Advances and Challenges in Process Identification, Sixth International Conference on Chemical Process Control, Tucson, Arizona, 2002; Rawlings, J. B., Ogunnaike, B. A., Eaton, J. W., Eds.; CACHE-AIChE: New York; 55-74. (24) Fortescue, T. R.; Kershenbaum, L. S.; Ydstie, B. E. Implementation of Self-tuning Regulators with Variable Forgetting Factors. Automatica 1981, 17 (6), 831-835. (25) Singh, A.; Badwe, A.; Patwardhan, S. C. Constrained Recursive Parameter Estimation for Adaptive Control, Proceedings of the 8th International Symposium On Dynamics And Control Of Process Systems (DYCOPS-2007), Cancun, Mexico, June 6-8, 2007. (26) Vachhani, P.; Rengaswamy, R.; Gangwal, V.; Narasimhan, S. Recursive Estimation in Constrained Nonlinear Dynamical Systems. AIChE J. 2005, 51, 3, 946-959. (27) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control; John Wiley: New York, 1989. (28) Birol, G.; Undey, C.; Cinar, A. A Modular Simulation Package for Fed-Batch Fermentation: Penicillin Production. Comput. Chem. Eng. 2002, 26, 1553-1565. (29) Downs, J. J.; Vogel, E. F. A Plant-Wide Industrial Process Control Problem. Comput. Chem. Eng. 1993, 17 (3), 245-255. (30) Ricker, N. L.; Lee, J. H. Nonlinear Modeling and State Estimation for the Tennessee Eastman Challenge Process. Comput. Chem. Eng. 1995, 19 (9), 983-1005.

ReceiVed for reView June 15, 2007 ReVised manuscript receiVed August 8, 2007 Accepted January 25, 2008 IE070823Y

Adaptive Model Predictive Control of Multivariable Time-varying Systems

Recommend Documents