Generalized Predictive Control Algorithm with Real-Time

Article pubs.acs.org/IECR

Generalized Predictive Control Algorithm with Real-Time Simultaneous Modeling and Tuning Yong Kuen Ho,*,† Hak Koon Yeoh,† and Farouq S. Mjalli‡ †

Chemical Engineering Department, Faculty of Engineering, University of Malaya, 50603 Kuala Lumpur, Malaysia Petroleum & Chemical Engineering Department, Sultan Qaboos University, Muscat 123, Oman

‡

ABSTRACT: A dual adaptation generalized predictive controller (DA-GPC) capable of synchronous modeling and tuning in real time is proposed. The distinctive feature of the proposed controller is the utilization of a second order discrete time controlled auto-regressive integrated moving average (CARIMA) model to capture all possible variations of numerator dynamics via the variable forgetting factor recursive least squares (VFF-RLS) algorithm. The evolving dynamic model is then used to tune, in real time, the value of the move suppression weight. Simulation studies on two benchmark process control problems, i.e., the Van de Vusse nonlinear isothermal stirred tank reactor (with dead time) and the strong acid−base neutralization process, revealed the improved performance of the DA-GPC scheme over the traditional model adaptive GPC and the conventional GPC schemes in terms of set point tracking. Further simulation studies revealed that the DA-GPC scheme was adept in rejecting load disturbances in the process. The good performance attained was attributed to the judicious adaptation of the move suppression weight. Moreover, a user-friendly aspect of the DA-GPC scheme is that the only parameter requiring manual tuning is the VFFRLS parameter (σ), which may span a large range without materially affecting the efficacy of this adaptive controller.

1. INTRODUCTION The performance of the linear generalized predictive controller (GPC)1,2 is closely dependent on two major components embedded within its structure, i.e., the model parameters and the tuning parameters. For strictly linear and time-invariant processes, the tasks of modeling and tuning are straightforward when designing the GPC controller. However, such trivial processes most often do not exist in reality. Instead, real-world processes are usually nonstationary in time and nonlinear in nature. Worse, unanticipated changes in the process dynamics can occur due to reasons which may not be obvious at the instance where urgent troubleshooting is required to recover the process. For such situations, one of the good solutions is to employ an adaptive GPC controller such that the real-time dynamics of the process can be accounted for in the computation of the optimized controller moves. The construction of the adaptive GPC controller to date is primarily achieved through the use of the recursive least squares (RLS) algorithm.3 In this approach, the RLS algorithm is used as a model estimator to identify the parameters of a discrete time linear input/output model at every time step based on the given data of the inputs and outputs of the process. The identified model parameters are then used to update and hence keep the internal model of the GPC controller in pace with the everevolving dynamics of the process. Hereafter this form of controller is referred to as the adaptive generalized predictive controller (A-GPC). Although A-GPC is more adept in controlling realistic nonstationary processes, it is limited by the degree to which the controller is properly tuned. In fact, the necessity to appropriately tune the predictive controller according to the different unique process models employed is evident in the multiple model predictive controller (MMPC) framework.4−6 Such crucial linkage between the model and the tuning © 2014 American Chemical Society

parameters of predictive controllers in the adaptive framework, however, is generally not found in RLS-based adaptive GPC controllers. One early exception is found in the work of Diaz et al.,7 in which the RLS algorithm was employed to update the model parameters of the GPC and the move suppression weight was concurrently tuned based on a simple decreasing law. Despite this attempt to “model and tune” in parallel, the basis of change in the tuning parameter was ad hoc in nature, requiring some understanding of the process dynamics and was not directly correlated with the estimated model parameters. A non ad hoc strategy was not available in the literature until recently, where Ho et al.8 used the RLS algorithm in conjunction with the first order plus dead time (FOPDT) tuning strategy9 to tune a GPC with an internal FOPDT model structure. The control scheme was shown to perform well on processes which can be approximated by the FOPDT model. Other than these exceptions, in general the implementations of RLS-based adaptive GPC controllers to date are restricted to model adaptation alone; see, e.g., Ho et al.,10 Moon et al.,11,12 and Corrêa et al.13 At the other extreme, studies on autotuning of predictive controllers (which are not based on the RLS algorithm)14−21 did not attempt to also update the model parameters. Motivated by the above, in this work we attempt to address the gap between the modeling and tuning components in the RLS-based adaptive GPC framework by extending our previous work8 beyond the application boundaries of FOPDT processes. As a proof of concept, only stable single input−single output (SISO) processes will be considered in this work. We name the Received: Revised: Accepted: Published: 9411

June 17, 2013 November 29, 2013 April 18, 2014 April 18, 2014 dx.doi.org/10.1021/ie401905w | Ind. Eng. Chem. Res. 2014, 53, 9411−9426

Industrial & Engineering Chemistry Research

Article

adaptive GPC controller which uses the RLS algorithm for synchronous updating of the model as well as the tuning parameters at every time step as the dual adaptation generalized predictive controller (DA-GPC). Considering the need for any proposed control scheme to be practical and beneficial to a wide range of processes, a few important characteristics are embedded in the design of the DA-GPC, namely the following: (i) a discrete time second order plus dead time (SOPDT) equivalent of the controller’s internal model is used so that a broader range of process dynamics beyond those of FOPDT systems can be covered, including those which exhibit inverse response and overshoot phenomena; (ii) model adaptation by the RLS algorithm enables the DA-GPC to handle commonly encountered process nonlinearities, i.e., slowly time-varying process dynamics and nonlinearities induced by differences in operating points; (iii) an autotuning strategy is formulated based on the tuning correlations proposed by Shridhar and Cooper9 to enable the DA-GPC to retune itself to obtain good process response and actuator moves; (iv) all the above objectives are achieved with relative implementation simplicity, with minimal user knowledge and intervention. As with all other advanced controller design methods, the shortcoming of the proposed algorithm is that some form of offline work is required to initiate the controller implementation. Nonetheless, such offline procedures involved minimal efforts. The organization of this paper is as follows. Section 2 gives a brief theoretical background of A-GPC, followed by the detailed structure of the DA-GPC in section 3. Section 4 compares the performance of the DA-GPC, A-GPC, and conventional GPC schemes, based on the results of two different benchmark process control case studies. Lastly, significant conclusions of this work are given in section 5.

(4)

ωk = Pk − 1 − γkψkTPk − 1

(5)

⎧ if tr(ωk /λk ) ≤ C ⎪ ωk / λk Pk = ⎨ ⎪ otherwise ⎩ ωk

(6)

θk̂ = θk̂ − 1 + γkεk

(7)

In eqs 2−7, ε is the prediction error, γ is the Kalman gain, λ ∈ (0, 1] is the forgetting factor, λmin is the minimum forgetting factor, and P is the covariance matrix, whereas C and σ (which adjusts the adaptivity of the VFF-RLS algorithm toward process parameter changes) are design parameters. For deterministic systems, usually a small amount of noise (σw, e.g., ∼σ/103) is added to the process output data (yk) entering the VFF-RLS algorithm to ensure that the system is sufficiently excited for successful closed loop parameter estimation.24,25 The regressor vector (ψ) and the vector of estimated model parameters (θ̂) with d as the bias parameter are given as ψ = [−yk − 1 −yk − 2 uk − 1 − D uk − 2 − D 1]

T

θ ̂ = [ a1 a 2 b1 b2 d ]T

(9)

y = H Δ u k − 1 + K Δ u k − 1 + Qy ⃗

⃗k

⃖

⃖k

H ∈ ℜ(N2 − N1) × M K ∈ ℜ(N2 − N1) × (D + 1)

and

(1)

Q ∈ ℜ(N2 − N1) × 3

where u is the process input, y is the process output, v is a stochastic noise variable with normal distribution and zero mean, D ≥ 0 is the discrete dead time expressed as integer multiples of the sampling time (ts), T(z−1) is the GPC design polynomial selected as unity for simplicity, and Δ = 1 − z−1. In eq 1, the subscript k is a nonnegative integer which denotes the sampling instance (k = 0, 1, 2, ...), while a1, a2, b1, and b2 are the model parameters which are adapted in real time by the variable forgetting factor recursive least squares (VFF-RLS) algorithm. The latter is an improved version of the RLS algorithm with a time-varying forgetting factor22 and a strategy to bound the trace of the covariance matrix at all times:23

while the positive integers N1, N2, and M ≤ (N2 − N1) are the minimum prediction horizon, maximum prediction horizon, and the control horizon, respectively. The vectors of past and future process variables at sampling instance k are defined according to the convention of Rossiter:28 y = [ yk + 1 yk + 2 ··· yk + N ]T 2

⃗k

(12)

Δ u k − 1 = [Δuk Δuk + 1 ··· Δuk + M − 1]T

(13)

Δ u k − 1 = [Δuk − 1 Δuk − 2 ··· Δuk − (D + 1)]T

(14)

⃖k

γk = Pk − 1ψk /[1 + ψkTPk − 1ψk]

(2)

(11)

y = [ yk yk − 1 yk − 2 ]T ⃗

εk = yk −

(10)

In eq 10,

[1 + a1z −1 + a 2z −2]yk

T θk̂ − 1ψk

(8)

The bias parameter d was estimated so that unbiased estimates of a1, a2, b1, and b2 can be generated.26 In addition to the above, the VFF-RLS algorithm used in this study was implemented through Bierman’s P = UDUT factorization algorithm27 to ensure the numerical stability of the recursive computations. Having obtained the instantaneous model parameters via the VFF-RLS algorithm, eq 1 is augmented to predict the future output trajectory by including variables for future instances:

2. ADAPTIVE GENERALIZED PREDICTIVE CONTROL (A-GPC) The A-GPC scheme used in this work assumes a discrete time SOPDT equivalent controlled auto-regressive integrated moving average (CARIMA) model,1 which is represented in the z-domain as

= [b1z −1 + b2z −2]uk − D + T (z −1)vk /Δ

T T ⎧ ⎪1 − {εk εk /[σ (1 + ψ Pk − 1ψ )]} if λk > λ min k k λk = ⎨ ⎪ otherwise ⎩ λmin

⃖

Imposing the constraints of umin ≤ Δ u k − 1 ≤ umax and Δumin ⃗ ≤ Δ u k − 1 ≤ Δumax, as well as assuming unity weights on the

(3)

⃗

9412

dx.doi.org/10.1021/ie401905w | Ind. Eng. Chem. Res. 2014, 53, 9411−9426


(

)

Article

parameters estimated by the VFF-RLS algorithm (i.e., a1, a2, b1, and b2) are used for the online computation of the move suppression weight (R) in addition to producing updates for the controller’s internal model (details below). The move suppression weight is chosen due to its efficacy in affecting closed loop performance.7,9,20,31 This is in line with the perspective of optimal control, where the tuning parameters are the weights of the objective function.25 Ideally, if tuning correlations similar to those proposed by Shridhar and Cooper9 based on the continuous time FOPDT model are available for the continuous time second order plus dead time model with numerator dynamics (SOPDTN), the model parameters estimated by the VFF-RLS algorithm can be used directly upon discrete−continuous time conversion to calculate the move suppression weight. However, the derivation of SOPDTN-based tuning correlations for direct calculation of the move suppression weight is nontrivial; hence an alternative route is taken here which enables retuning of the move suppression weight using the FOPDT-based tuning correlations proposed by Shridhar and Cooper9 without compromising the model order employed within the structure of the controller. In this approach, the estimated discrete time model parameters by the VFF-RLS algorithm are first converted to one of the possible forms of continuous time SOPDTN equivalent. Upon conversion, the corresponding FOPDT model parameters are obtained by minimizing the integral of squared error between the SOPDTN and FOPDT unit step responses. The solutions to the minimization problem are used to calculate the move suppression weight of the controller by reverting to the FOPDT-based tuning correlations proposed by Shridhar and Cooper.9 Note, however, that this conversion to the FOPDT equivalent is used solely for tuning, as depicted in Figure 1. For the controller’s internal model, the structure of a discrete SOPDT equivalent is preserved. Sections 3.1−3.4 elaborate the procedures involved in the DA-GPC algorithm. 3.1. Primary Screening of Estimated Model Parameters. The model parameters identified by the VFF-RLS algorithm can be arranged in the form of the following discrete time SOPDT transfer function:

output residuals r k − y , where r k = [rk+1 rk+2 ··· rk+N2]T is ⃗

⃗k

⃗

the vector of future set points (ri), the prediction equation in eq 10 is used to set up the following quadratic programming (QP) constrained optimization problem: min

{ 12 Δu

⃗

T

[HTH + R̅ ]Δ u k − 1

⃗ k−1

Δ uk −1

⃗

}

+ [HT(K Δ u k − 1 + Qy − r k)]T Δ u k − 1 ⃖

⃖k

⃗

⃗

subject to

ΩΔu k − 1 − μ ≤ 0 ⃗

(15)

In the above, R̅ ∈ ℜM × M is a positive definite weighting matrix containing the move suppression weight (R) in the diagonals, i.e., R̅ ij = Rδij. The construction of matrix Ω and vector μ can be found in the standard literature.28,29 The solution to the QP problem above yields the optimal Δ u k − 1, of which only the first ⃗ element is implemented in the process. In the A-GPC scheme, the procedures of estimating the model parameters, constructing the prediction equation, and computing the optimized control moves are repeated at every time step. However, the tuning parameters, i.e., N1, N2, M, and R, are not adapted throughout the entire implementation; rather, the one-off selection of their values is done offline by selecting N1 = D, while the rest are selected according to the tuning strategy proposed by Hinde and Cooper30 as well as Shridhar and Cooper.9

3. DUAL ADAPTATION GENERALIZED PREDICTIVE CONTROL (DA-GPC) A distinguishing feature of the DA-GPC scheme from the AGPC scheme is the ability to “model and tune” in parallel at every time step. Figure 1 further highlights this point by showing that the A-GPC is a subset of the DA-GPC scheme. To enable the controller to retune itself in real time, the model

G(z −1) = z −D{[b1z −1 + b2z −2]/[1 + a1z −1 + a 2z −2]} (16)

Prior to using these discrete time model parameters in the DAGPC or the A-GPC schemes, the stability of this instantaneous model must be screened by evaluating the poles of eq 16, given as z1 = −0.5a1 + 0.5 a12 − 4a 2

(17)

z 2 = −0.5a1 − 0.5 a12 − 4a 2

(18)

Only model parameters with stable poles (i.e., |z1| < 1 and |z2| < 1) are used for model adaptation in the controller. Additionally, the gain of eq 16 as well as the absolute feedback error (e) is checked by the following equations:

Figure 1. Simplified diagram illustrating the salient differences between A-GPC and DA-GPC schemes. In the figure, u is the process input, y is the process output, a1, a2, b1, and b2 are the model parameters of a SOPDT equivalent CARIMA model, and R is the move suppression weight.

K p = [b1 + b2]/[1 + a1 + a 2]

(19)

e = |r − y|

(20)

where the difference between the set point (r) and the output (y) at the current instance is the feedback error. For processes in which the sign of the process gain matters, the sign of Kp in eq 19 is also screened. Further, in the event where the absolute 9413



Article

feedback error in eq 20 is smaller than a certain threshold (i.e., e < ξ), denoting a near steady state behavior, the controller is decoupled from the parameter estimator. This decoupling is prudent as it is unnecessary (nor desirable with RLS) to continue modifying the model parameters near steady state. In this work, the threshold ξ was selected as 1% of the significant variation (e.g., the magnitude of typical changes in set point or in the nominal process output) plus the average magnitude of noise. For example, if the process operates between output values of 1.2, 0.8, 0.6, etc. and there is no noise in the process, then the significant variation is of the order 10−1, and ξ = 0.01 × 10−1 + 0 = 1 × 10−3. If the model parameters failed any of the primary screening procedures above, the previous model and tuning parameters are retained in the controller setup until new and permissible estimates are available, while procedures discussed in sections 3.2−3.4 are bypassed. 3.2. Conversion of Discrete Time SOPDT Model to the Continuous Time SOPDTN Model. As the tuning correlations used in this work (to be discussed in section 3.4) are presented in the continuous time domain, it is necessary that eq 16 be converted to its continuous time equivalent: G(s) = K p(τns + 1)e−χs /[(τ1s + 1)(τ2s + 1)]

ζSOPDTN ∈ ℜ+ ⎧>1 (overdamped) ⎪ + ⎪ if τ1 , τ2 ∈ ℜ and τ1 ≠ τ2 ⎪ ⎪ ⎪1 (critically damped) = ⎨ if τ1 , τ2 ∈ ℜ+ and τ1 = τ2 ⎪ ⎪ ξ, values of τ1, τ2, τn, and χ are subsequently calculated: τ1 = −ts/ln( −0.5a1 + 0.5 a12 − 4a 2 )

(22)

τ2 = −ts/ln( −0.5a1 − 0.5 a12 − 4a 2 )

(23)

⎧− t /ln(− b /b ) if ln( − b /b ) ∈ ℜ 2 1 2 1 ⎪ s and ln(− b2 /b1) ≠ 0 ⎪ τn = ⎨ ⎪0 if ln(− b2 /b1) ∉ ℜ ⎪ or ln(− b2 /b1) = 0 ⎩

(24)

χ = Dts

(25)

(28)

In practice, a much narrower range of the values of ζSOPDTN than those given in eq 27 necessitates controller tuning. In our implementation, considering that 0.2 ≤ ζSOPDTN ≤ 5 is the typical range for chemical processes, 33 for values of ζSOPDTN outside this range (including complex values of ζSOPDTN), only model adaptation is implemented while the controller is not retuned. 3.3. Approximating the SOPDTN Model with a FOPDT Model for Tuning. After discrete−continuous model conversion, the order of the SOPDTN model in eq 28 is reduced to obtain its FOPDT equivalent. Specifically, only the FOPDT dead time (θp) and time constant (τp) are sought, as the FOPDT gain remains identical to eq 19. For this purpose, although there are simple methods such as the “half rule”34,35 and approximations by Taylor series expansions,32 their applications are limited to systems with real and distinct poles only (i.e., ζSOPDTN > 1), and the approximations deviate from optimality in some cases.34 Since we were unsuccessful in uncovering any published general methods to reduce the model order for all cases of ζSOPDTN, here we propose to minimize the integral of squared error between the unit step responses for the SOPDTN and FOPDT models to reduce the order of eq 28: min{Ij(τp , θp)} τp , θp

The dimensionless damping factor (ζSOPDTN) is calculated using the following expression: ζSOPDTN = [τ1 + τ2]/[2 τ1τ2 ]

Ij(τp , θp) = (26)

t = θp

∫t =χ +

Only positive real values of ζSOPDTN which fulfill one of the following conditions are useful for autotuning:

(29)

[3−1{Gj(s)/s}]2 dt

t = t∞

∫t =θ

[3−1{Gj(s)/s} − 3−1{H(s)/s}]2 dt

p

(30) 9414



Article

where 3−1{·} is the inverse Laplace transform of {·}, t is the physical quantity denoting time, t∞ is the time for both the process responses to reach steady state and beyond, and H(s) is the FOPDT model: H(s) = K pe−θps /[τps + 1]

technique which is well-known for its quadratic convergence property36 is used. Defining h = 1, 2, 3, ... as the number of iterations and given the initial guess values τp,0 and θp,0, eq 34 is iterated until ∥xh − xh−1∥ < 10−6 and ∥Fh−1∥ < 10−6: xh = xh − 1 − Jh − 1−1Fh − 1

(31)

(34)

where:

The choice of the objective function presented in eq 30 is best illustrated with a graphical example. Figure 2 shows an

xh − 1 = [ τp, h − 1 θp, h − 1]T Jh − 1

(35)

⎡ ∂F1(τp, h − 1 , θp, h − 1)/∂τp ∂F1(τp, h − 1 , θp, h − 1)/∂θp ⎤ ⎥ =⎢ ⎢ ∂F (τ ⎥ , θ )/ ∂ τ ∂ F ( τ , θ )/ ∂ θ p 2 p, h − 1 p, h − 1 p⎦ ⎣ 2 p, h − 1 p, h − 1 (36) T

Fh − 1 = [ F1(τp, h − 1 , θp, h − 1) F2(τp, h − 1 , θp, h − 1)]

The analytical expressions of J and F can be obtained using symbolic manipulation engines such as Mathematica (http:// www.wolfram.com/mathematica/). The algorithms for finding good initial guess values τp,0 and θp,0 as well as a fault detection scheme are presented in the Appendix. By virtue of the form of the denominator of eq 21, the analytical expressions for eqs 32, 33, and 36 are less bulky and are hence less prone to programming errors. Similarly, as shown in the Appendix, more compact and elegant expressions for the value of the minimum point of the SOPDTN model can be obtained, i.e., eqs A.1−A.3. 3.4. Automated Tuning Strategy. As alluded to in the beginning of section 3, only the move suppression weight (R) is autotuned while the rest of the tuning parameters are fixed; hence the optimal FOPDT model parameters obtained are used only in the computation of R. The fixed values of minimum prediction horizon (N1) and maximum prediction horizon (N2)9 are calculated as follows:

Figure 2. Unit step responses of SOPDTN (τn > 0) and FOPDT models, where χ is the dead time of the SOPDTN model, θp is the dead time of the FOPDT model, t∞ is the time required for the transients to reach steady state, 3−1{·} is the inverse Laplace transform of {·}, Gj(s) is the SOPDTN transfer function, and H(s) is the FOPDT transfer function.

attempt to fit a SOPDTN unit step response with τn > 0 by a FOPDT unit step response. For the transients presented in Figure 2, a good approximation of θp should begin with θp = χ as the value of θp should never be smaller than χ to avoid underestimation of the process dead time. In general, if θp ≥ χ, the resultant integral of squared error between the two curves up to t = χ (region I) is zero and can be removed from the cost function. For the remaining transients, it is desirable that the integral of the squared error between the two curves from t = χ to t = t∞ be minimized by adjusting values of τp and θp. Since a zero output for the FOPDT model is obtained from t = χ to t = θp (region II) due to the criterion of θp ≥ χ, the effect of the H(s) term can be safely excluded from the cost function in this region. As such, full integration of the squared error is only evaluated for t = θp to t = t∞ (region III) and the final objective function to be minimized is the sum of all the objective functions for regions I−III. Generally, the cost function developed here applies also to cases where τn ≤ 0 since θp ≥ χ applies to these cases (i.e., considering the inverse response behavior exhibited for τn < 0 and the sluggish initial transients for τn = 0). The optimal values of τp and θp which minimize the cost function in eq 30 can be obtained by solving for the roots of the following first order partial derivatives: F1(τp , θp) = ∂{Ij(τp , θp)}/∂τp = 0

(32)

F2(τp , θp) = ∂{Ij(τp , θp)}/∂θp = 0

(33)

(37)

N1 = D

(38)

N2 = CEIL[(5τp*/ts) + D*]

(39)

where D* is defined as D* = CEIL[θp*/ts] and CEIL[·] is the nearest higher integer to [·]. In addition, the fixed value of the control horizon (M)30 is calculated as M = CEIL[τp*/ts]

(40)

In the above, τ*p and θ*p (not to be confused with τp and θp) are the base case FOPDT model parameters identified from offline system identification32 (e.g., by using the open loop process reaction curve method). The automated tuning of R is done by using the tuning correlations proposed by Shridhar and Cooper9 for M > 1: f = [M /500][(3.5τp/ts) + 2 − (M − 1)/2]

(41)

R = fK p2

(42)

In eqs 41 and 42, we do not present the case for M = 1 because the sampling time (ts) in eq 40 is usually selected to be smaller than the base case FOPDT time constant (τ*p ), resulting in M > 1. In summary, Figure 3 shows the sequence of evaluations of the model and tuning parameters in the DA-GPC algorithm presented throughout sections 3.1−3.4. Note that the A-GPC scheme without automatic tuning is implemented in this work using a similar structure, but with the exclusion of the

In principle, any root-finding algorithm can be used to solve the above equations. In this work, the Newton−Raphson (NR) 9415



Article

move suppression weight. As such, the internal model assumed by the GPC, i.e., a discrete SOPDT equivalent of the CARIMA model, remains identical regardless of whether A-GPC or DAGPC is deployed.

4. SIMULATION RESULTS Two benchmark case studies for process control were used to validate the performance of the control algorithms presented in this work. The first is a Van de Vusse isothermal continuous stirred tank reactor (CSTR),37 which was chosen mainly for its inherent inverse response characteristics within a specified range of process input, although the process is also nonlinear across different operating regions. The second case study concerns a highly nonlinear neutralization process involving strong acid and strong base,38 which is the standard benchmark test for any control algorithm focusing on nonlinearity. 4.1. Case Study I: Concentration Control of Nonlinear Van De Vusse Isothermal Continuous Stirred Tank Reactor (CSTR). Consider the following series-parallel Van de Vusse reaction scheme:37 k1

k2

(43)

T→U→V k3

(44)

2T → W

where ki (i = 1, 2, 3) are the reaction rate constants for the three reactions. These reactions take place in an isothermal CSTR, where the model equations are39 dC T(t )/dt = [FR (t − Φ)/VR ][C TF(t ) − C T(t )] − k1C T(t ) − k 3C T 2(t ) Figure 3. Flowchart of the DA-GPC algorithm presented throughout sections 3.1−3.4. The procedures enclosed within the double dotted lines (corresponding to the same in Figure 1), excluded from the AGPC implementation, are the core of the DA-GPC scheme.

(45)

dC U(t )/dt = [−FR (t − Φ)/VR ]C U(t ) + k1C T(t ) − k 2C U(t )

(46)

In the above, FR is the feed flow rate of T, CTF is the concentration of T in the feed stream, CT and CU are the concentrations of T and U respectively in the reactor, VR is the volume of the reactor content which is kept constant during

procedures enclosed within the dotted line. It should be emphasized that the enclosed procedures are distinguishing features of the DA-GPC deployed solely for the tuning of the

Table 1. Initial Values and Nominal Operating Values (Denoted by Additional Subscript “n”) for the Van de Vusse CSTR Reactor Model As Well As Parameter Settings for Initialization of the DA-GPC, A-GPC, and GPC Controllersa,b process model43

VFF-RLS

parameter

value

k1 (s−1) k2 (s−1) k3 (L mol−1 s−1) VR (L) CT(0) (mol L−1) CU(0) (mol L−1) CTFn (mol L−1) FRn (L s−1) CTn (mol L−1) CUn (mol L−1) Φ (s)

5/360 5/180 1/360 700 2.9175 1.10 10 10 3.712 1.226 5.06

parameter θ̂T0 P0 σ (mol2 L−2) σw (mol2 L−2) C D ts (s) λmin

GPC and Variants value

parameter

value

01×5 1 × 104I5×5 1.2 × 10−4 2 × 10−7 4 × 104 1 5 0.9

θGPC R N1 N2 M FRmin (L s−1) FRmax (L s−1) ΔFRmin (L s−1) ΔFRmax (L s−1) ts (s) t∞c (s) ξ (mol L−1)

θ* 1.58 × 10−4 2 43 7 0 15 0.75 −0.75 5 5 × 104 1 × 10−3

θGPC = [ a1 a 2 b1 b2 D ]

θ* = [− 1.61 6.482 × 10−1 − 6.486 × 10−3 7.338 × 10−3 1]

θ* is identified from a step change in FR from 10 L s−1 to 9 L s−1. θ* changes with time for the A-GPC, whereas both θ* and R change with time for the DA-GPC. bAll VFF-RLS parameters, t∞ and ξ are not applicable to the GPC. cDA-GPC only. a

9416



Article

Figure 4. Performance of the DA-GPC, A-GPC, and GPC schemes in tracking a series of set point changes in the concentration of U (CU). The DA-GPC and A-GPC schemes were activated at 1400 s.

Figure 5. Manipulated variable profile, i.e., the feed flow rate of T (FR), produced by the DA-GPC, A-GPC, and GPC schemes in tracking a series of set point changes in the concentration of U (CU). The DA-GPC and A-GPC schemes were activated at 1400 s.

operation, and Φ is the dead time of the process which was not present in the original model equations but is introduced in this work to challenge the controllers. From a control perspective, the dead time could simulate the event that changes in the manipulated variable FR in the reactor feed may not be translated immediately into changes in the controlled variable CU which is measured downstream of the reactor. Here, CTF is a possible source of disturbance. The initial values and nominal operating conditions used to simulate the reactor dynamics as well as the controller settings for the DA-GPC, A-GPC, and GPC schemes are shown in Table 1. This reactor system is not only nonlinear, but also exhibits inverse response behavior (i.e., at FR ≤ 15 L s−1 with parameter values in Table 1) and a notorious change in the process gain39,40 (FR > 15 L s−1). As such, if the operating region of FR is not bounded, the Van de Vusse reactor will exhibit input multiplicity.6,39 This is a very challenging control problem which had been addressed through multiple model based strategies,6,39−41 among others.40,42 The RLS algorithm is in principle able to track such unanticipated and time-varying changes in the dynamics, with the advantage of circumventing the need for offline development of multiple models. However, the target here is more modest, focusing on regions with

inverse response and nonlinear dynamics more commonly seen in chemical process plants. Thus, we restrict the study to FR ≤ 15 L s−1. For subsequent controller design purposes, two sets of base case FOPDT and SOPDTN model parameters were preidentified offline from a step change in FR from 10 to 9 L s−1: base case FOPDT model: C U(s) [mol/L] 2.2527 × 10−2e−48.166s = FR (s) [L/s] 33.156s + 1

(47)

base case SOPDTN model: C U(s) [mol/L] 2.2466 × 10−2(− 40.889s + 1)e−5.2617s = FR (s) [L/s] 532.5941s 2 + 46.1883s + 1 (48)

For the VFF-RLS algorithm, all parameters have default values except for σ, ts, and D. The value of σ was adjusted based on a general guide that it should be decreased for more nonlinear processes such that the forgetting factor λ ≪ 1 is obtained to improve adaptivity, whereas ts was selected as 9417



Article

Figure 7. Manipulated variable profile (i.e., feed flow rate of T, FR) exhibited by the DA-GPC and A-GPC schemes for different values of VFF-RLS design parameter σ.

Figure 6. Temporal evolution of the estimated process gain (Kp), estimated FOPDT time constant (τp), and the corresponding evolution of move suppression weight (R) for the DA-GPC scheme for the Van de Vusse reactor.

different operating ranges. Other than that, chattering behavior was also observed for the manipulated variable FR as the fixed tuning parameters were unable to cater to the wide range of reactor operation. This situation improved with the deployment of model adaptation in the A-GPC scheme, where offsets were no longer observed in the CU profile. Despite the improved response, significant chattering was still observed in the FR profile due to the absence of continuous retuning as the model parameters changed, which consequently produced slight oscillatory behavior in the CU profile. With the DA-GPC scheme, all such shortcomings were overcome with appropriate continuous retuning of the controller, giving a smooth CU profile and near absence of chattering in FR (except for a little between 3400 and 3600 s). Figure 6 shows the temporal evolution of the values of Kp, τp, and R for the DA-GPC scheme. From Figure 6, the values of Kp and τp, which are representative of the dynamical changes in the process, changed significantly over time. The DA-GPC scheme proved to be capable of updating the value of R accordingly. If such a significant variation in the value of R was needed to produce good controller performance, clearly the A-GPC scheme with model adaptation alone and a constant value of Rtuned for the nominal operating value of CUwas inadequate. This contrast in performance can further be explained by considering the evolution of the value of parameter Kp, which has a relatively higher effect on the value of R compared to τp according to eqs 41 and 42. As the value of Kp increased over time (denoting a more sensitive process), the value of R should be increased to avoid overactive controller output. Since the A-GPC was tuned for the low gain

approximately 1/10 of the predetermined FOPDT time constant (τp*), and D (not D*) was calculated as D = CEIL[χ*/ts] (i.e., χ* is the base case SOPDTN dead time identified offline). The internal model of the controller for all schemes was initialized as the discrete counterpart of the base case SOPDTN model, and all controller tuning parameters were initialized using eqs 38−42 with the base case FOPDT model. In addition, t∞ was selected as a large value and the constraints on the controller slew rates, i.e., ΔFRmin and ΔFRmax, were chosen as ±5% of the difference between FRmax and FRmin. Figures 4 and 5 show the performance of the various control schemes in tracking consecutive set point changes for CU. To demonstrate the effect of the different adaptive mechanisms in the DA-GPC and A-GPC schemes over time, all simulation runs were initialized with the GPC controller and the VFF-RLS algorithm was activated at the onset of simulation while the deployment of the DA-GPC and A-GPC schemes were initiated at 1400 s purely for clarity in the display so that a distinct variation in the performance can be seen in the different adaptive strategies employed. From Figures 4 and 5, the response exhibited by the GPC scheme is poor and offsets were observed as the set point was moved away from the nominal CU at 1.226 mol L−1 (e.g., with an absolute maximum of approximately 8% from the set point value between 4000 and 6000 s). The offsets occurred because the internal model of the controller, which was identified within a narrow operating range, was unable to represent the dynamics of the reactor at 9418



Article

Figure 8. Performance of the DA-GPC scheme in (a) rejecting a 5% step increment in the nominal feed concentration of T (CTF), introduced at 7700 s as a disturbance to the process and (b) responding to white noise addition of zero mean and variance = 1 × 10−5 mol2/L2. The concentration of U (CU) is the controlled variable, while the feed flow rate of T (FR) is the manipulated variable.

Table 2. Initial Values and Nominal Operating Values (Denoted by Additional Subscript “n”) of the Neutralization Model As Well As Parameter Settings for Initialization of the DA-GPC, A-GPC, and GPC Controllersa,b process model parameter −1

CAFn (mol L ) CBFn (mol L−1) η (mol L−1) Kw VN (L) FAn (L s−1) FBn (L s−1) pHn

VFF-RLS value −4

6.3096 × 10 3.1623 × 10−4 0 10−14 2 1.667 × 10−3 3.325 × 10−3 7

parameter θT0̂ P0 σ σw C D ts (s) λmin

GPC and Variants value

parameter

value

01×5 1 × 104I5×5 2.5 × 10−6 2 × 10−7 4 × 104 0 30 0.9

θGPC R N1 N2 M FBmin (L s−1) FBmax (L s−1) ΔFBmin (L s−1) ΔFBmax (L s−1) ts (s) t∞c (s) ξ

θ* 5.133 × 109 1 49 10 0 0.005 2.5 × 10−4 −2.5 × 10−4 30 5 × 104 1 × 10−2

θGPC = [ a1 a 2 b1 b2 D ] θ* = [− 1.753 7.659 × 10−1 1.035 × 104 − 9.199 × 103 0 ] η = CA(0) − C B(0)

θ* is identified from a step change in FB from 3.325 × 10−3 L s−1 (199.5 mL min−1) to 3.333 × 10−3 L s−1 (200 mL min−1). θ* changes with time for the A-GPC whereas both θ* and R change with time for the DA-GPC. bAll VFF-RLS parameters, t∞ and ξ are not applicable to the GPC. cDAGPC only. a

region resulting in a tight controller, it was unable to perform well in the high gain regions of the process. While one could contend that tuning the A-GPC for the high gain region of the process instead (resulting in a sluggish “worst case” controller) will possibly mitigate the chattering behavior as well as the oscillatory response observed, this sluggishness risks not meeting process requirements. Having a controller which can retune itself in real time such as the DA-GPC is definitely a better alternative. The closed loop response produced by the DA-GPC and AGPC schemes was further studied for different values of the VFF-RLS design parameter σ. Figure 7 shows the performance of the controllers in terms of the manipulated variable profiles. At values of σ spanning 3 orders of magnitude, the performance of the DA-GPC scheme remained better than that of the A-

GPC scheme. Figure 7 also suggests that adjusting the value of σ alone, while influencing the performance of the A-GPC, is inadequate. Proper determination of the move suppression weight R, for which the value is not obvious in many cases, still must be done by the control engineer. However, with the DAGPC scheme, the move suppression weight is selected automatically: one only needs to adjust the value of σ. Moreover, since the performance of the DA-GPC is relatively insensitive to the values of σ in comparison to the A-GPC scheme, the selection of σ can be left to the nonexpert engineers. The regulatory performance and the robustness of the DAGPC scheme toward noisy measurements on the CSTR were also tested in this study. For the former, a 5% step increment in CTF was introduced at 7700 s as a disturbance to the process 9419



Article

Figure 10. Manipulated variable profile, i.e., flow rate of the base feed stream (FB), produced by the DA-GPC, A-GPC, and GPC schemes in tracking a series of set point changes in the pH of the solution. The DA-GPC and A-GPC schemes were activated at 2700 s.

Figure 9. Performance of the DA-GPC, A-GPC, and GPC schemes in tracking a series of set point changes in the pH of the solution. The DA-GPC and A-GPC schemes were activated at 2700 s.

following the previous set point test, while for the latter, white noise with zero mean and variance = 1 × 10−5 mol2/L2 was added to the process. Figure 8a shows the regulatory performance of the controller, where the DA-GPC was able to arrest the effect of CTF on CU within 300 s. Good manipulated variable profiles were also obtained. In the noisy environment test, Figure 8b shows that the DA-GPC responded satisfactorily. Note that in Figure 8b, a huge jump in the set point was also introduced. This revealed a test on the controller’s response when the constraint FRmin = 0 L s−1 was active. The DA-GPC readjusted itself quickly and actuator saturation was not prolonged. This was not the case for the AGPC and GPC (not shown), where frequent actuator saturations were observed within the same duration. All these showed that the DA-GPC scheme was not only adept in both servo and regulatory control of the nonlinear Van de Vusse CSTR, but that it was also reasonably robust toward noise addition and instances when a constraint was active. 4.2. Case Study II: pH Control of Nonlinear Neutralization Process Involving Strong Acid and Strong Base. We further consider the control of a pH neutralization process involving a strong acid (e.g., hydrochloric acid, HCl) and a

strong base (e.g., sodium hydroxide, NaOH) in a stirred tank, where the model equations are given as38 d[CA(t ) − C B(t )]/dt = [1/VN]{FA(t ) CAF(t ) − FB(t ) C BF − [FA(t ) + FB(t )] [CA(t ) − C B(t )]}

(49)

[C H+(t )]2 − [CA(t ) − C B(t )]C H+(t ) − K w = 0

(50)

pH(t ) = −log10[C H+(t )]

(51)

In eqs 49−51, CA is the concentration of the acid in the tank, CB is the concentration of base in the tank, CAF is the concentration of the acid feed stream, CBF is the concentration of the base feed stream, FA is the flow rate of the acid feed stream, FB is the flow rate of the base feed stream, VN is the volume of the tank content which is kept constant during operation, Kw is the water dissociation constant at the operating temperature, and CH+ is the activity of hydrogen ions. The controlled variable in this case is the pH of the solution in the tank, whereas FB is the manipulated variable (i.e., acid is titrated 9420



Article

Figure 12. Closed loop pH profile exhibited by the A-GPC and DAGPC schemes for different values of VFF-RLS design parameter σ. The corresponding integral absolute error (IAE) for each run is also given.

Figure 11. Temporal evolution of the estimated process gain (Kp), estimated FOPDT time constant (τp), and corresponding evolution of move suppression weight (R) for the DA-GPC scheme for the neutralization process.

The remaining degree of freedom, i.e., the VFF-RLS design parameter σ, was adjusted using the approach presented for the Van de Vusse CSTR example. For the GPC component, initialization of the internal model, the tuning parameters, and ξ was also done as described for the Van de Vusse CSTR example, while t∞ was retained as default. The constraints on the controller slew rates, i.e., ΔFBmin and ΔFBmax, were selected as ±5% of the difference between FBmax and FBmin. Figures 9 and 10 show the performance of the various control schemes in tracking consecutive set point changes in the pH of the solution. All simulation runs were initialized with the GPC controller and the VFF-RLS algorithm was activated at the onset of simulation, while the deployment of the DAGPC and A-GPC schemes were initiated at 2700 s to show the difference (if any) in the performance of the different adaptive controllers. From the results, the GPC scheme performed reasonably well up to 25 000 s, albeit some overshoots were observed in the response. This is expected as the GPC scheme was designed around the nominal pH 7 (i.e., FOPDT and SOPDTN models were identified in that region); thus moving the set point from a high gain region (i.e., pH ≈7) to a low gain region (i.e., pH ≈5) had little impact on the performance of the GPC scheme. However, this is not the case when moving the set point in the opposite direction, where a sharp and significant overshoot was observed after 25 000 s. For the AGPC scheme, unstable closed loop response was observed for 22 000 s and beyond, and a sudden fluctuation in the flow rate of the base feed stream was observed in the manipulated variable profile. This shows that while the A-GPC attempted to

with base), while FA and CAF are the possible sources of disturbance. The neutralization process described above poses a challenging process control problem because the gain of the process changes drastically as the pH of the solution increases from an acidic condition to the end point. For this simulation, the initial values and nominal operating conditions of the neutralization process, as well as the controller settings for the DA-GPC, A-GPC and GPC schemes are given in Table 2. Note that, for the VFF-RLS algorithm, the values of θ̂T0 , P0, σw, C, and λmin were retained as the default values used in the Van de Vusse CSTR example, while the selection of the values of D and ts were also as alluded to in the same example. The base case FOPDT model and SOPDTN model for the neutralization process were identified from a step change in FB from 199.5 to 200 mL min−1, which corresponds to pH values of 7.0 and 7.7, respectively: base case FOPDT model: pH(s) 90.326 × 103 = FB(s) [L/s] 291.1s + 1

(52)

base case SOPDTN model: pH(s) 90.538 × 103(255.21s + 1)e−1.3559s = FB(s) [L/s] 6.227 × 104s 2 + 553.5296s + 1

(53) 9421



Article

Figure 13. Performance of the DA-GPC scheme in (a) rejecting a 5% step increment in the nominal concentration of the acid feed stream (CAF) and the flow rate of the acid feed stream (FA), introduced separately at 30 000 s as a disturbance to the process, and (b) responding to white noise addition of zero mean and variance = 1 × 10−3. The pH of the solution is the controlled variable, while the flow rate of the base feed stream (FB) is the manipulated variable.

For the sake of completeness, the performance of the DAGPC scheme in regulating a 5% step increment in the nominal concentration of the acid feed stream (CAF) or the nominal flow rate of the acid feed stream (FA) is given in Figure 13a. The disturbances were introduced at 30 000 s in separate simulation runs following the previous set point test. The results showed that both disturbances were rejected within 500 s and the manipulated variable profiles exhibited were sufficiently smooth and tractable. In addition to the above, as demonstrated by Figure 13b, the proposed scheme showed a good degree of robustness toward white noise addition (with zero mean and variance = 1 × 10−3) and closed loop instability was not observed.

adapt to a low gain internal model within that duration, the move suppression weight which was tuned based on a high gain model (i.e., resulting in a sluggish controller for a low gain process) was simply inappropriate. With the dual adaptation capability, the DA-GPC scheme showed good closed loop response as well as a manipulated variable profile throughout the entire duration of simulation and none of the drawbacks exhibited by the GPC and A-GPC schemes were observed. Figure 11 shows the temporal evolution of the estimated values of Kp, τp, and R for the DA-GPC scheme. The values of Kp and τp, which determine how the value of R will evolve, changed drastically with time. Beyond 22 000 s, Kp changed from high gain to low gain, denoting a decrease in the process sensitivity (requiring a smaller R to avoid having a sluggish response), while τp increased with time, denoting a more sluggish process (requiring a smaller R to increase the speed of response). The combined effects of Kp and τp changed the value of R substantially, giving the DA-GPC an improved performance compared to the A-GPC scheme with a constant R = 5.133 × 109. Moreover, although sudden spikes in the parameter values were observed in Figure 11, the impact on the performance of the controller was small as the parameters were continuously updated, and the effect of possibly erroneous parameter estimates was not propagated in time. As σ is a critical parameter in determining the forgetting factor (λ) profile for the parameter estimator, it could be debated that, for such a highly nonlinear neutralization process, the poor performance of the A-GPC scheme shown thus far could be due to inappropriate selection of σ and had nothing to do with the choice of R. However, Figure 12 shows that, under different values of σ spanning 3 orders of magnitude, the AGPC scheme with constant R produced consistently poor and sluggish performance, whereas with proper values of R, the DAGPC scheme gave satisfactory performance without closed loop instability. To further ascertain that the poor performance of the A-GPC was due to inappropriate tuning and not to the misselection of σ, the simulation using the A-GPC controller was repeated with different ranges of σ and the poor results obtained (not shown) were consistent with those in Figure 12.

5. CONCLUSIONS Challenged by two benchmark case studies, the newly proposed DA-GPC scheme showed remarkable improvements to the existing A-GPC scheme. With a discrete second order CARIMA model structure, the DA-GPC scheme is capable of controlling complicated processes which exhibit numerator dynamics. Moreover, the automated modeling and tuning procedure enables the controller to deal with nonlinear processes of which conventional linear controllers are incapable of handling. In implementation, the DA-GPC scheme mimics a one degree of freedom controller, where only a single parameter, i.e., the VFF-RLS design parameter σ, needs to be adjusted while all other parameters can be set to their nominal default values. It should be noted, however, that the DA-GPC formulation does not cater to processes with varying dead times spanning a large range. Nevertheless, for processes with slight changes in the time delays, choosing a conservative but constant dead time would still make DA-GPC work in principle, albeit less nimbly. While the DA-GPC scheme is designed for SISO systems here as proof of concept, our future efforts will examine to what extent a decentralized DA-GPC scheme is able to handle multivariable systems, in particular the loop interactions. 9422



Article

Table 3. Procedures for Obtaining the Initial Guess Values for the FOPDT Dead Time θp,0 and the FOPDT Time Constant τp,0 for Different Ranges of Values of the SOPDTN Numerator Time Constant (τn) When ζSOPDTN ∈ [0.2, 5]a

a

For values outside this range, the search is not performed. All time related quantities have units of seconds. 9423



■

Article

APPENDIX: ALGORITHMS FOR DETERMINING THE INITIAL GUESS VALUES AS WELL AS FAULT DETECTION SCHEME FOR THE NEWTON−RAPHSON TECHNIQUE The Newton−Raphson (NR) technique in section 3.3 requires an initial guess vector [θp,0 τp,0] which is sufficiently near the

can be bracketed. On bracketing point S, a one-step regula falsi (RF) yields the guess value for θp,0. This is just one of the multiple possibilities to bracket point S. In the procedure described above, the crucial step is to find the minimum point L. Since χ is known, the time to reach point L is obtained by adding Q to χ. For the overdamped and critically damped cases, Q is found by solving for the root of d[ySOPDTN(t)]/dt = 0 analytically.32 The resulting expressions for the overdamped and critically damped cases are shown in eqs A.1 and A.2, respectively: Q = [τ1τ2/(τ1 − τ2)] ln[(1 − τn /τ2)/(1 − τn /τ1)]

(A.1)

Q = τ1τn /(τn − τ1)

(A.2)

For the case of underdamped behaviour, casting the complex poles of G3(s) in the form of τ1, τ2 ∈ C and τ2 = τ*1 , a similar expression can be obtained by substituting the complex τ1 and τ2 conjugates into eq A.1. Further simplification leads to eq A.3: Q = [(a 2 + b2)/2b] tan−1[−2XY /(X2 − Y 2)]

Figure 14. Simulated response curve ySOPDTN (t ) = 3−1{Gj(s)/s} within the DA-GPC structure for illustrating the search for the initial guess of FOPDT dead time (θp,0) and time constant (τp,0), where Kp is the model gain, τn < 0 is the SOPDTN numerator time constant, and χ is the SOPDTN dead time.

Table 4. Procedures To Trap Erroneous Newton−Raphson (NR) Solutions, Categorized by the Ranges of Values of the SOPDTN Numerator Time Constant (τn) range of τn τn < 0 τn > 0 τn = 0

(A.3)

where a = Re(τ1), b = Im(τ1), X = a + b − aτn, and Y = bτn. For τp,0, a good guess value is the interval between point S and point V, since the FOPDT time constant is defined as the time for the response to reach 63.2% of the final steady state value.32 By shifting the response curve downward by 0.632Kp, again point V will be the zero of the shifted curve, and the root bracketing technique for θp can be used to locate it. The challenge is to find the lower bracket of point V. An obvious choice is θp,0, or failing which, the lower bracket in the RF procedure that yielded θp,0. On finding the location of point V at t = κ, assign the guess value τp,0 = κ − θp,0. For τn > 0 (which is similar to a first order response except that overshoot can occur when τn is large), the search for the initial guess values of θp,0 and τp,0 is similar to that presented for τn < 0, except that θp,0 can be initialized directly as χ plus a small constant δ < 0.1 (in order to avoid multiplication by zero in the subsequent search procedure) and bracketing is only required for the search of τp,0. In summary, Table 3 gives the numerical procedures for obtaining the initial guess values for different ranges of τn. In addition to giving the NR good initial guess solutions for θp and τp as described above, an additional safeguard is used to trap erroneous NR solutions, as outlined in Table 4. Such solutions will be suboptimal as they do not truly satisfy the minimization required by eq 29. However, their existence for short periods of time, if at all, is not expected to materially affect the controller performance. 2

procedure to trap and salvage erroneous NR solutions If θp ∉ [t0,Q, t1,Q] or τp ≤ 0 Use multistep RF to find θp within [t0,Q, t1,Q] and κ within [T0,Q, T1,Q], then τp = κ − θp. If θp < χ or τp ≤ 0 Set θp = χ. Then, use multistep RF to find κ within [T0,Q, T1,Q], then τp = κ − θp. If θp < χ or τp ≤ 0 Depending on the value of ζSOPDTN: Set θp = θp,0 as calculated by eq A.15, A.16, or A.17 in Table 3. Set τp = τp,0 as calculated by eq A.32, A.33, or A.34 in Table 3.

intended solutions. Here, this initialization is based on three different ranges of the SOPDTN numerator time constant, namely τn = 0, τn < 0, and τn > 0. For τn = 0, the initial guess values for θp,0 and τp,0 are given by adapting the empirical approximations proposed by Panda et al.33 The remaining categories of τn < 0 and τn > 0 share similar search methodology in principle; therefore, we will only elaborate the ideas behind the search technique for τn < 0. Figure 14 represents a simulated response curve based on eq 28 for τn < 0. A good initial guess for θp,0 is point S (θp,0 ≈ χ + φ), which is a zero of the response curve. Thus, a root-finding approach can be used to locate it. Point S is bracketed by ySOPDTN (t ) = 3−1{Gj(s)/s} (where j = 1, 2, and 3 refers to the overdamped, critically damped, and underdamped cases, respectively) of different signs and always occurs after the response reaches a minimum at point L. Hence, by marching forward from point L with a sufficiently small interval, point S

■

2

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +6012-631 3536. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge the Faculty of Engineering from the University of Malaya as well as the Petroleum & Chemical Engineering Department from Sultan Qaboos University for providing full technical support for this work. In addition, we also thank the University of Malaya for the financial support in the form of a UMRG grant (RG134-11SUS) and a PPP grant (PG091-2012B), both of which H.K.Y. is the principal 9424



Article

investigator. We also express our special thanks to the Bright Sparks Unit from the University of Malaya for the financial support to Y.K.H.

Δ u k − 1 = vector of past slew rates (not including current ⃖ value) Δumin = lower limit for future slew rate constraints Δumax = upper limit for future slew rate constraints VR = volume of the Van de Vusse reactor content y = process output y = vector of future outputs (not including current value)

■

NOMENCLATURE a1, a2, b1, b2 = estimated model parameters C = VFF-RLS design parameter CA = concentration of the acid for the neutralization tank case study CAF = concentration of the acid feed stream for the neutralization tank case study CB = concentration of base for the neutralization tank case study CBF = concentration of the base feed stream for the neutralization tank case study CEIL[·] = nearest higher integer to [·] CH+ = activity of hydrogen ions CT = concentration of T for the Van de Vusse reactor case study CTF = concentration of T in the feed stream for the Van de Vusse reactor case study CU = concentration of U for the Van de Vusse reactor case study D = discrete dead time d = bias parameter D = diagonal matrix D* = base case discrete dead time calculated from θ*p e = absolute feedback error FA = flow rate of the acid feed stream for the neutralization tank case study FB = flow rate of the base feed stream for the neutralization tank case study FR = flow rate of T for the Van de Vusse reactor case study H = matrix associated with the future slew rates in the GPC prediction equation K = matrix associated with the past slew rates in the GPC prediction equation k1, k2, k3 = reaction rate constants for the Van de Vusse reaction scheme Kp = gain of the discrete time SOPDT model/SOPDTN model/FOPDT model Kw = water dissociation constant M = control horizon N1 = minimum prediction horizon N2 = maximum prediction horizon P = covariance matrix P0 = initial values for the covariance matrix Q = matrix associated with the past outputs in the GPC prediction equation R = move suppression weight R̅ = positive definite diagonal weighting matrices containing R r k − 1 = vector of future set points ⃗ t = physical quantity denoting time T(z−1) = GPC design polynomial t∞ = time for process responses to reach steady state and beyond ts = sampling time u = process input umin = lower limit for future input constraints umax = upper limit for future input constraints Δ u k − 1 = vector of future slew rates (including current value)

⃗ k−1

y

⃖ k−1

= vector of past outputs (including current value)

z1, z2 = poles of the discrete time SOPDT model Greek Symbols

■

γ = Kalman gain ε = prediction error ζSOPDTN = damping factor Φ = dead time of the Van de Vusse reactor θ̂ = vector of estimated model parameters θT0̂ = initial values for vector of estimated model parameters θp = dead time for the FOPDT model θp,0 = initial guess values for θp θ*p = base case FOPDT dead time λ = forgetting factor λmin = minimum forgetting factor ξ = threshold for e σ = VFF-RLS design parameter σw = process measurement noise τp = time constant for the FOPDT model τp,0 = initial guess values for τp τ*p = base case FOPDT time constant τ1, τ2 = denominator time constants for the SOPDTN model τn = numerator time constant for the SOPDTN model χ = dead time for the SOPDTN model χ* = base case SOPDTN dead time ψ = regressor vector

REFERENCES

(1) Clarke, D. W.; Mohtadi, C.; Tuffs, P. S. Generalized predictive controlPart I. The basic algorithm. Automatica 1987, 23 (2), 137− 148. (2) Clarke, D. W.; Mohtadi, C.; Tuffs, P. S. Generalized predictive controlPart II Extensions and interpretations. Automatica 1987, 23 (2), 149−160. (3) Seborg, D. E.; Edgar, T. F.; Shah, S. L. Adaptive control strategies for process control: A survey. AIChE J. 1986, 32 (6), 881−913. (4) Dougherty, D.; Cooper, D. A practical multiple model adaptive strategy for single-loop MPC. Control Eng. Pract. 2003, 11 (2), 141− 159. (5) Dougherty, D.; Cooper, D. A practical multiple model adaptive strategy for multivariable model predictive control. Control Eng. Pract. 2003, 11 (6), 649−664. (6) Aufderheide, B.; Bequette, B. W. A variably tuned multiple model predictive controller based on minimal process knowledge. In American Control Conference, Arlington, VA, June 25−27, 2001; American Automatic Control Council: Evanston, IL, USA, 2001; pp 3490−3495. (7) Diaz, C.; Dieu, P.; Feuillerat, C.; Lelong, P.; Salome, M. Adaptive predictive control of dissolved oxygen concentration in a laboratoryscale bioreactor. J. Biotechnol. 1995, 43 (1), 21−32. (8) Ho, Y. K.; Mjalli, F. S.; Yeoh, H. K. Generalized predictive control with dual adaptation. Chem. Eng. Sci. 2012, 84, 479−493. (9) Shridhar, R.; Cooper, D. J.; Tuning, A. Strategy for Unconstrained SISO Model Predictive Control. Ind. Eng. Chem. Res. 1997, 36 (3), 729−746. (10) Ho, Y. K.; Mjalli, F. S.; Yeoh, H. K. Multivariable adaptive predictive model based control of a biodiesel transesterification reactor. J. Appl. Sci. 2010, 10 (12), 1019−1027.

⃗

9425



Article

(33) Panda, R. C.; Yu, C. C.; Huang, H. P. PID tuning rules for SOPDT systems: Review and some new results. ISA Trans. 2004, 43, 283−295. (34) Skogestad, S. Simple analytic rules for model reduction and PID controller tuning. J. Process Control 2003, 13 (4), 291−309. (35) Skogestad, S. Corrections and comments from the author on the paper “Simple analytic rules for model reduction and PID controller tuning”, Journal of Process Control 13 (2003) 291−309. J. Process Control 2004, 14 (4), 465. (36) Chapra, S. C.; Canale, R. P. Numerical Methods for Engineers, 5th ed.; McGraw-Hill: New York, 2006. (37) Van de Vusse, J. G. Plug flow type reactor versus tank reactor. Chem. Eng. Sci. 1964, 19, 964. (38) Ali, E. pH control using PI control algorithms with automatic tuning method. Trans. Inst. Chem. Eng. 2001, 79 (Part A), 611−620. (39) Sistu, P. B.; Bequette, B. W. Model predictive control of processes with input multiplicities. Chem. Eng. Sci. 1995, 50 (6), 921− 936. (40) Aufderheide, B.; Bequette, B. W. Extension of dynamic matrix control to multiple models. Comput. Chem. Eng. 2003, 27, 1079−1096. (41) Kuure-Kinsey, M.; Bequette, B. W. Multiple Model Predictive Control of Nonlinear Systems. In Nonlinear Model Predictive Control; Magni, L., Raimondo, D., Allgöwer, F., Eds.; Springer: Berlin, 2009; Vol. 384, pp 153−165. (42) Aufderheide, B.; Bequette, B. W. A comparison of fundamental model-based and multiple model predictive control. Presented at the 40th IEEE Conference on Decision and Control, Orlando, FL, USA, 2001. (43) Balaguer, P.; Alfaro, V.; Arrieta, O. Second order inverse response process identification from transient step response. ISA Trans. 2011, 50, 231−238.

(11) Moon, S. M.; Clark, R. L.; Cole, D. G. The recursive generalized predictive feedback control: theory and experiments. J. Sound Vib. 2005, 279 (1−2), 171−199. (12) Moon, S. M.; Cole, D. G.; Clark, R. L. Real-time implementation of adaptive feedback and feedforward generalized predictive control algorithm. J. Sound Vib. 2006, 294 (1−2), 82−96. (13) Corrêa, N. A.; Corrêa, R. G.; Freire, J. T. Self-tuning control of egg drying in spouted bed using the GPC algorithm. Drying Technol.: Int. J. 2002, 20 (4), 813−828. (14) Al-Ghazzawi, A.; Ali, E.; Nouh, A.; Zafiriou, E. Online tuning strategy for model predictive controllers. J. Process Control 2001, 11, 265−284. (15) Han, K.; Zhao, J.; Qian, J. A novel robust tuning strategy for model predictive control. In 6th World Congress on Intelligent Control and Automation, Dalian, China, June 21−23, 2006; Institute of Electrical and Electronics Engineers: New York, 2006; pp 6406−6410. (16) Kawai, F.; Ito, H.; Nakazawa, C.; Matsui, T.; Fukuyama, Y.; Suzuki, R.; Aiyoshi, E. Automatic tuning for model predictive control: Can particle swarm optimization find a better parameter?. In 22nd IEEE International Symposium on Intelligent Control, Singapore, Oct 1−3, 2007; Institute of Electrical and Electronics Engineers: New York, 2007; pp 646−651. (17) Ali, E. In Automatic Tuning of Model Predictive Controllers Based on Fuzzy Logics. In Proceedings of IASTED Conference on Control and Applications, Banff, Canada, June 27−29, 2001; International Association of Science and Technology for Development: Calgary, Canada, 2001; pp 136−142. (18) van der Lee, J. H.; Svrcek, W. Y.; Young, B. R.; Tuning, A. Algorithm for Model Predictive Controllers Based on Genetic Algorithms and Fuzzy Decision Making. ISA Trans. 2008, 47, 53−59. (19) Liu, W.; Wang, G. Auto-tuning procedure for model predictive controller. Proceedings of the 2000 IEEE International Conference on Systems, Man, and Cybernetics, Nashville, TN, Oct 8−11, 2000; Institute of Electrical and Electronics Engineers: New York, 2000. (20) Li, S.; Du, G. In Online tuning scheme for Generalized Predictive Controller via simulation-optimization. IEEE International Conference on Fuzzy Systems, Honolulu, HI, USA, May 12−17, 2002, Institute of Electrical and Electronics Engineers: New York, 2002; pp 1381−1386. (21) Valencia-Palomoa, G.; Rossiter, J. A. Programmable logic controller implementation of an auto-tuned predictive control based on minimal plant information. ISA Trans. 2010, 50, 92−100. (22) Fortescue, T. R.; Kershenbaum, L. S.; Ydstie, B. E. Implementation of self-tuning regulators with variable forgetting factors. Automatica 1981, 17 (6), 831−835. (23) Cordero, A. O.; Mayne, D. Q. Deterministic convergence of a self-tuning regulator with variable forgetting factor. IEEE Proc. D: Control Theory Appl. 1981, 128 (1), 19−23. (24) Landau, I. D.; Lozano, R.; M’Saad, M.; Karimi, A. Adaptive Control: Algorithms, Analysis and Applications; Springer: Berlin, 2011. (25) Mikleš, J.; Fikar, M. Process Modelling, Identification and Control; Springer-Verlag: Berlin, 2007. (26) Shah, S. L.; Cluett, W. R. Recursive least squares based estimation schemes for self-tuning control. Can. J. Chem. Eng. 1991, 69 (1), 89−96. (27) Bierman, G. J. Measurement updating using the U-D factorization. Automatica 1976, 12 (4), 375−382. (28) Rossiter, J. A. Model-Based Predictive Control; CRC Press: Boca Raton, FL, 2003. (29) Camacho, E. F.; Bordons, C. Model Predictive Control; SpringerVerlag: Berlin, 1999. (30) Hinde, R. F., Jr.; Cooper, D. J. A pattern-based approach to excitation diagnostics for adaptive process control. Chem. Eng. Sci. 1994, 49 (9), 1403−1415. (31) McIntosh, A. R.; Shah, S. L.; Fisher, D. G. Analysis and tuning of adaptive generalized predictive control. Can. J. Chem. Eng. 1991, 69, 97−110. (32) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control, 2nd ed.; John Wiley & Sons, Inc.: New York, 2004. 9426


Generalized Predictive Control Algorithm with Real-Time

Recommend Documents