8966
Ind. Eng. Chem. Res. 2009, 48, 8966–8979
PROCESS DESIGN AND CONTROL Development of ARX Models for Predictive Control Using Fractional Order and Orthonormal Basis Filter Parametrization Muddu Madakyaru,† Anuj Narang,‡,§ and Sachin C. Patwardhan*,‡ Systems and Control Engineering and Department of Chemical Engineering, Indian Institute of Technology, Bombay, Powai, Mumbai 400076, India
Among various industrial vendors of model predictive control (MPC) schemes, ARX appears to be a popular choice of model structure while models for predictive control schemes are being developed. These models are, however, nonparsimonious in the number of model parameters. As a consequence, the length of the data required to keep the variance errors low while using the conventional ARX structure is significantly large. Thus, if it is possible to reparametrize the ARX model such that fewer parameters are required at the identification stage, then it is possible to reduce the length of identification data and thereby reduce the cost involved in model identification exercise. In this work, we explore the possibility of reparametrizing ARX models using the fractional-order differential operators (FO-ARX) and orthonormal basis filters (OBF-ARX). We also propose a novel approach for identification of time delay matrix from multivariate data using ARX and OBF-ARX models. The efficacy of the proposed modeling technique is demonstrated by conducting simulation studies on the benchmark Shell control problem. Analysis of the simulation results reveals that, when compared with the conventional high-order ARX structure, the FO-ARX and the OBF-ARX are better model parametrizations when the data length is less. In particular, OBF-ARX parametrization is able to estimate the time delay matrix from multivariate data quite accurately. The experimental studies establish the feasibility of using the proposed FO-ARX and OBF-ARX models for formulating MPC schemes. 1. Introduction In the past two decades, model predictive control (MPC) technology has become a touchstone tool for industrial advanced process control in refining and petrochemical industries and is now beginning to attract interest from other process industries as well.1 The need to solve complex control problems in the process industries has led to significant research activity in the field of model predictive control over the past two decades. Given a dynamic model of the system under consideration, the issues such as controller design to ensure stability and robustness of MPC schemes or efficient online computations have been very well researched in the process control literature. The key component of any MPC scheme is the dynamic model used for carrying out online predictions. The development of the dynamic model is of paramount importance for achieving good closedloop performance. The model developed at the beginning of a project may require maintenance when the operating conditions become significantly different from the design conditions. Thus, if it is desired to achieve uniform operating performance over a wide operating range using a linear MPC scheme, then it often becomes necessary to reidentify the model parameter. However, model identification for MPC is a costly affair as it involves perturbing the plant to generate data relevant for model identification, which results in the loss of production during the perturbation time. One major concern while reidentifying models is the length of data required for parameter estimation, which translates to the time required for conducting identification * To whom correspondence should be addressed. Tel.: +91-2225767211. Fax: +91-22-25726895. E-mail:
[email protected]. † Systems and Control Engineering. ‡ Department of Chemical Engineering. § Current address: Department of Chemical and Materials Engineering, University of Alberta, Canada T6G2V4.
experiments. Thus, it is necessary to develop a modeling scheme that can generate good-quality models using relatively short data length so that the cost involved in model identification can be reduced. When compared to the controller design aspects, relatively less literature is available on cost-effective model identification approaches for developing/maintaining models in an MPC scheme. Industrial applications of model predictive control rely mostly on linear empirical models obtained by employing time series analysis approaches.2 There are several methods reported in the literature for developing stochastic time series models3,4 Shook and Shah5 have proposed an approach that modifies the conventional time series model identification strategy to suit the specific needs of predictive control. Shouche et al.6 have developed a simultaneous MPC and identification (MPCI) approach that proposes to identify MPC relevant ARX type models online while minimally disturbing the process. McArthur and Zhan7 have recently developed a multistage approach for identification of MPC relevant models under closed-loop conditions based on time series analysis. There are other system identification methods, such as subspace identification method,8 that have come up recently and have been employed to develop models suitable for MPC. Lakshminarayanan et al.9 have presented an approach based on projection on partial least squares for the development of multivariable time series models. Biao and co-workers10-12 have established the feasibility of using subspace approaches for developing state space models for MPC schemes. They have developed open-loop as well as closed-loop identification schemes. The development of a conventional time series model using the prediction error method, except for the ARX structure, requires a nonlinear optimization problem to be solved. On the other hand, the main
10.1021/ie8009439 CCC: $40.75 2009 American Chemical Society Published on Web 09/16/2009
Ind. Eng. Chem. Res., Vol. 48, No. 19, 2009
advantage of the subspace approach is that the parameter identification is carried out using projections. However, the number of parameters required to be identified by this approach is, in general, significantly large when compared to the conventional time series models such as ARMAX or BJ. As a consequence, the data set required to generate good parameter estimates is significantly large. Most MPC relevant model identification approaches discussed above have been developed with the motivation of generating a deterministic component of the model that is better suited for long-range predictions. These approaches typically generate the model parameter estimates by minimizing some norm of either the filtered or the multistep prediction errors. As a consequence, the noise model obtained in this process may not be a good representation of the true unmeasured disturbance behavior. Recently, Patwardhan et al.13 have shown that the regulatory performance of an MPC formulation can be significantly improved if the unmeasured disturbance model identified from the operating data by minimizing one-step prediction errors is directly used in the controller formulation. These improvements can be attributed to the fact that the noise model can be used for predicting the effects of past unmeasured disturbance on the future output behavior, thereby improving the prediction ability of the model. Thus, from the viewpoint of achieving better regulatory performance, the noise model developed together with the deterministic model by minimizing unfiltered one step ahead prediction errors can also be viewed as an MPC releVant model. Finite impulse response (FIR) and autoregressive with exogenous input (ARX) model structures are widely used for identifying control relevant multivariable models from plant data.2 The main advantage of ARX structure is that the resulting model parameter estimation problem can be solved analytically. Another advantage of ARX structure over FIR structure is that ARX facilitates estimation of the noise model simultaneously with the deterministic model. To capture the noise part really well, the model order is typically chosen quite high.14 Moreover, unlike the FIR model, the ARX model structure can be used for identification of both unstable and stable systems. These advantages of the ARX model have made it popular among various industrial vendors of MPC.2 It may be noted that the variance errors in parameter estimates are directly proportional to the number of model parameters and inversely proportional to the length of data. As a consequence, the length of the data required to keep the variance errors low while using the conventional ARX structure is significantly large. Also, it is known that using a model with more degrees of freedom (or parameters) is not always better as it can result in the modeling of nonexistent dynamics and noise characteristics. Thus, if it is possible to reparametrize the ARX model such that fewer parameters are required at the identification stage, then it is possible to reduce the length of identification data. A direct consequence of reduction in the data length is that the loss of production during the model identification exercise is reduced and, effectively, the cost involved in model identification exercise is reduced. In this context, the following two model parametrization techniques appear promising: • fractional-order (FO) differential equations;15 • generalized orthonormal basis filters (GOBF).16 The fractional calculus (FC) based on generalized fractionalorder integrodifferential operators is a very old topic in mathematics.18 Since its inception, the theory of fractional calculus developed primarily as a pure theoretical field of mathematics has been useful only to mathematicians. However,
8967
in the past few decades many authors have pointed out that using the fractional system theory can actually describe the properties of various real-world physical systems in a better way. Podlubny and co-workers have contributed a significant amount of work aimed at this subject.19,20 Hartley and Lorenzo21 and Poinot and Trigeassou22 have discussed in detail the application of fractional system theory. A few of the applications of fractional system theory are long distributed lines, electrochemical processes, dielectric polarization, viscoelastic materials, and in control theory. The fractional-order derivative of a variable can be viewed as a limit of an infinite series involving integer-order derivatives. Consequently, a fractional-order model can be viewed as a parsimonious representation of an infiniteorder ODE model. We propose to exploit this feature of a fractional order model to develop a high-order ARX model, which is parsimonious in parameters to be identified. It may be noted that, in the control literature, fractional-order representation has been mainly used for the development of models with output error (OE) structure. The possibility of using the fractional-order representation for unmeasured disturbance modeling has not been exploited. In recent years there has been growing interest in the use of orthogonal basis filters (OBF), which form a complete set or a orthonormal basis in the space of all stable transfer functions, for representing process dynamics.16 The orthogonal filter approximations provide a simple and elegant method of representing open-loop stable systems. In fact, these models can be looked upon as compact representations of convolution type models, which have been widely used in MPC schemes. Heuberger et al.16 have shown how OBF parametrization can significantly reduce the number of parameters to be estimated when compared to an FIR model. In MPC literature, however, the OBF parametrization has mostly been employed to develop models with output error structure.23,24 The possibility of using GOBF parametrization of the ARX model has not been exploited in the literature. It may be noted that an ARX model, in fact, can be viewed as a canonical parametrization of a observer and is required to be stable even if the deterministic dynamics is unstable. The OBF-based parametrization provides a direct method to ensure stability of the observer at the identification stage. This work is aimed at the development of reparametrized ARX type models that are suitable for developing linear predictive control schemes. We explore the feasibility of reparametrizing ARX models using the fractional-order differential operators (FO-ARX) and the orthonormal basis filters (OBF-ARX). To begin with, we develop nested optimizationbased schemes for identification of parameters of FO-ARX and OBF-ARX models using input-output perturbation data. We also propose a novel approach for identification of time delay matrix from multivariate data using ARX and OBF-ARX models. The efficacy of the proposed modeling technique is demonstrated by conducting simulation studies on the benchmark Shell control problem. We demonstrate that the quality of the conventional ARX model diminishes as data length is reduced. On the other hand, the proposed FO-ARX and OBFARX parametrizations generate good-quality models even with the reduced data length. To establish the feasibility of using the proposed FO-ARX and OBF-ARX models for formulating MPC schemes, experimental studies are also conducted on a benchmark heater-mixer setup.25 The latter part of the paper is organized in the following sequence. Section 2 presents the rationale behind the proposed model parametrizations and the model parameter estimation
8968
Ind. Eng. Chem. Res., Vol. 48, No. 19, 2009
schemes. Sections 3 and 4 present results of simulation and experimental studies, respectively. The main conclusions reached through the analysis of these results are presented in the last section. 2. Reparametrization of ARX Models In the model identification exercise typically carried out before developing a predictive control scheme, the problem at hand is to construct a suitable time series using a set of sampled sequence of input and output vectors ΣN ) {y(k), u(k)} k ) 1, 2, ..., N
y(k) ) G(z-1) u(k) + H(z-1) e(k)
where y ∈ R and u ∈ R . Given input and output data set ΣN collected from a plant by injecting input perturbations in open loop, the problem of identifying a linear time series model can be stated as finding a linear operator φ[...] r
noise properties using the same set of poles, the resulting model may be unnecessarily high in order. While developing time series models from input output data, the first task is to select model structure and model order. In practice, the true model order is rarely known. Also, the length of data available for model identification is finite, and the data contain unmeasured disturbances and measurement noise. As a result of these factors, the two basic types of errors can occur during parameter estimation exercise, namely, the variance error and the bias error.3 For the general model structure, represented as (6)
m
y(k) ) φ[u(k - 1), ..., u(1), y(k - 1), ..., y(1), θ] + e(k) (1) such that a suitable norm of model residuals {e(k)} (k ) 1, ..., N) is minimized with respect to (wrt) parameter vector θ. The steps involved in model development are (a) planning and injecting perturbations in the plant such that the relevant modes of the plant are excited and (b) selecting a suitable model structure and identifying model parameters from the perturbation data collected. It is important to note that the two steps are not quite independent and the choice of model structure has a bearing on the duration of the perturbation exercise for model identification. To understand their interdependence, consider a general form of single-input-single-output (SISO) time series model y(k) ) G(z-1) u(k) + H(z-1) e(k)
(2)
where {e(k)} is a zero mean white noise sequence. For the purpose of parameter estimation, this model is rearranged in one-step predictor form as follows: yˆ(k|k - 1) ) H(q, θ)-1 G(q, θ) u(k) + (1 - H(q, θ)-1)y(k) ) Wu(q, θ) u(k) + Wy(q, θ) y(k) (3) The model parameters θ are identified by minimizing an objective function of the form
it has been shown that the variances of transfer function estimates as a function of frequency behave according to the following relationship:3 Φ (eiω) ˆ (eiω)] = n V cov[G N Φ (eiω)
(7)
u
where n represents the number of model parameters, N ˆ (...) represents the estimated represents the data length, and, G deterministic model. To keep the variance errors under check, it is desirable to keep the signal-to-noise ratio (Φu(eiω)/ΦV(eiω)) high. However, the signal-to-noise ratio is governed by considerations such as tolerable perturbations in the outputs during identification tests. Once the input excitation signals are selected, the signal-to-noise ratio is fixed. The other factor that influences the variance error is the ratio n/N, which should be kept small. To keep variance errors in check, we have the option to manipulate either the number of parameters(n) or the data length (N) in such a way that the ratio n/N remains small. The use of higher order ARX models for identification implies that a large number of parameters needs to be estimated, and, consequently, we require a large data set to keep the n/N ratio small. Thus, if we can explore a new way of parametrization of the ARX models with a lower number of parameters, then we can work with smaller data lengths while keeping the variance errors in check. The ARX model given by eq 5 can be expressed as follows: n
y(z) ) z-d
∑bz i
i)1
n
-i
u(z) -
∑az i
-i
y(z) + e(z)
(8)
i)1
N
J)
∑
1 [y(k) - yˆ(k|k - 1)]2 N k)1
(4)
This approach is called prediction error method, and details can be found in Ljung3 and Sodderstrom and Stoica.4 Different parametrizations of G(q) and H(q) result in different model structures, such as ARX, ARMAX, and Box-Jenkins. A significant advantage of ARX model over ARMAX and Box-Jenkins models is that the resulting predictor eq 3 assumes the form n
yˆ(k|k - 1) )
n
∑ b u(k - d - i) - ∑ a y(k - i) i
i)1
i
(5)
i)1
which is linear in parameters θ. Thus, unlike the other two model structures, the resulting parameter estimation problem can be solved analytically using linear least-squares approach. However, since the ARX model describes both the system dynamics and
If it is desired to reparametrize the above model, then we can think of two possible approaches. (1) The first is to parametrize the parameter sets {ai} and {bi} using some alternate representation: The discrete wavelet transform based compression method proposed by Nikolaou and Vuthandam17 can be viewed as a strategy belonging to the first approach. Discrete wavelet transforms can provide a compact and parsimonious reparametrization of the parameter sets {ai} and {bi}. In this work, we achieve this by parametrizing ARX models using fractional-order differential equations. (2) The frequency domain operator, z-i, has a short memory, which results in nonparsimonius representation. Parsimony could be achieved by using operators with longer memory. In this work, we propose to achieve this using the orthonormal basis filters (OBF). The main advantage of selecting OBF is that these filters have longer memory when compared to the operator z-i. Details of the proposed reparametrization schemes are presented in the subsequent subsections.
Ind. Eng. Chem. Res., Vol. 48, No. 19, 2009
2.1. Parametrization Using FO Models. Fractional calculus is a generalization of integration and differentiation to noninteger-order fundamental operator aDtλ, where a and t are the limits of the operation and λ is the fractional integration order. The continuous integrodifferential operator is defined as
{
dλ dtλ λ aDt ) 1
∫
t
a
R(λ) > 0 R(λ) ) 0 (dτ)-λ R(λ) < 0
}
(9)
When the real part of λ is negative, it is in fact a fractional integral operator which then performs integration over the limits of the operation, and if the real part of λ is positive, then the operation is differentiation. So, we can refer to fractional derivative or fractional integral indifferently. The two most popular definitions used for the general fractional differential are the Gru¨nwald-Letnikov (GL) discrete form of the definition and the Riemann-Liouville (RL) definition.18 The GL definition for function f(t) is given as Dλ λεC
f(t) ) lim hf0
1 hλ
∑ [(-1) (λi ) f(t - ih)] ∞
i
(10)
i)0
such that λ > 0. Here,
()
λ Γ(λ + 1) ) i Γ(i + 1) Γ(λ - i + 1)
(11)
Γ(...) is the well-known Euler’s Gamma function and h is the finite sampling step. For convenience, Laplace domain notion is usually used to describe the fractional integrodifferential operation as follows: L{0Dλt f(t)} ) sλF(s) where λ ∈ R. Here it is assumed that all initial values are zero. The general SISO FO system representation in continuous time domain is shown in (12). R1
[τ1D
+ τ2D
R2
Rr
0
