5006
Ind. Eng. Chem. Res. 2002, 41, 5006-5015
Frequency-Domain-Based, Control-Relevant Model Reduction for Nonlinear Plants C. Shreesha,†,‡ Ravindra D. Gudi,*,§ and P. S. V. Nataraj‡ Systems and Control Engineering Group, Department of Electrical Engineering, and Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India
Chemical plants exhibit nonlinear dynamics that need to be captured accurately for modelbased control. This paper focuses on a frequency-domain-based, control-relevant model reduction strategy for nonlinear process plants. The generalized frequency response function proposed for the nonlinear plants has been used to analyze the characteristics of nonlinear structures in the frequency domain, from a control-relevant viewpoint. Control-relevant model reduction is sought by minimizing the mismatch between a higher order, nonlinear proxy plant and candidate reduced order nonlinear model structures. The minimization is proposed to be done by prefiltering of identification data to shape the mismatch in a control-relevant manner. For these prefilters as well, candidate linear and nonlinear structures are evaluated and an algorithm to design the prefilters to reflect the servo and regulatory performance specifications is proposed. The model reduction methodology is used along with a nonlinear internal model control (NIMC) algorithm and is demonstrated on a simple simulation as well as on a representative nonlinear chemical reactor studied in the literature. 1. Introduction Most chemical process plants exhibit nonlinear dynamics. The modeling and control of these plants to achieve certain desired specifications is a challenging task. The conventional method of designing a controller for such a nonlinear plant is to find a local linear model which describes the plant dynamics in the neighborhood of the chosen operating point. A controller is designed based on this local linear model and is expected to give the specified performance around the chosen point of operation. However, because of the changes in operating conditions, plant parameters, and the effect of noise and disturbances, the operating point is susceptible to change and, thus, the system may operate unsatisfactorily. Hence, methodologies for nonlinear model estimation and model-based control design have been proposed for such nonlinear plants. Two issues need to be clearly addressed in the task of identifying models for such plants. First, the identified models need to be simple in structure, and hence the choice needs to be based on parsimony considerations. The second issue is the inevitable model-plant mismatch that comes into play because of the parsimony in the model structure. For linear plants, these issues have been elegantly addressed through the application of control-relevant identification by Rivera et al.,1-4 Shook et al.,5 and Kwok and Shah.6 Control-relevant identification can be performed in open- or closed-loop schemes. For linear systems, a comparison of open-loop verses closed-loop identification, when the system itself * To whom all correspondence should be addressed. Email:
[email protected]. Telephone: (91) (22) 576 7204. Fax: (91) (22) 572 6895. † Present address: Department of Electrical and Electronics Engineering, N.M.A.M. Institute of Technology, Post Nitte, Karkala, Udupi, Karnataka, India. ‡ Systems and Control Engineering Group, Department of Electrical Engineering. § Department of Chemical Engineering.
belongs to the model class, has been presented by Hjalmarsson et al.7 They show that the best control performance is achieved via identification with the intended controller in place in a closed loop. For nonlinear plants, structures such as the Hammerstein, the Weiner, and the Volterra series have been proposed for the identification. The Hammerstein and the Weiner structures typically emerge by a combination of linear dynamic and nonlinear static elements. The identification method for such models is discussed in refs 8-11. Using these models, controllers having linearlike characteristics can be synthesized using a gain scheduling method for the Hammerstein models and a variable transformation method for the Weiner models.12 The Volterra models capture the true plant dynamics as a set of impulse response coefficients associated with the appropriate order kernels.13-15 However, even for a second-order Volterra model with a finite lag, the number of impulse response coefficients required to represent the dynamics is very large. Various approaches for the estimation of second-order Volterra models are proposed in refs 16-20. An advantage of the Volterra model is that the nonlinear internal model control (NIMC) design method based on this requires only inversion of the linear kernel and, therefore, the Volterra model has emerged as an alternate to conventional modeling.14,18,21 The large number of coefficients required to represent even a second-order Volterra model still makes the controller design a tedious task, and the structure of the controller could be complex because of the presence of a large number of parameters. Therefore, model reduction into simpler model structures that would retain the relevant dynamics is desirable. When satisfaction of the specification on the control performance is explicitly considered to determine the relevant dynamics, the reduction strategy is termed control-relevant model reduction (CRMR). Zheng and Zafiriou22 have proposed a methodology for control-relevant identification of a second-order
10.1021/ie010946x CCC: $22.00 © 2002 American Chemical Society Published on Web 08/28/2002
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002 5007
Volterra model for a given nonlinear process. They have formulated the problem by modeling uncertainty due to the truncation error and kernel mismatch. The restricted complexity models such as the Hammerstein models can be approximated into Volterra series under certain conditions. Ling and Rivera21 have proposed a CRMR methodology for the second-order Volterra models to obtain a control-relevant Hammerstein model that is substantially simpler. They adapted the NIMC strategy, based on Volterra models, to develop an algorithm for CRMR. An alternate representation of nonlinear models is by a nonlinear autoregressive with moving average and exogenous input (NARMAX) structure, which yields a linear-in-parameter, identification methodology. Ling and Rivera12 have proposed a method for estimating higher order NARX model coefficients using the classical orthogonal least-squares method of structure selection and parameter estimation, proposed by Chen et al.23 Then, considering the resulting NARX model as a plant proxy, an equivalent Volterra model for the plant was estimated in terms of the NARX model coefficients following an algorithm proposed by Diaz and Desrochers.24 The CRMR of the resulting Volterra model was carried out in a manner similar to that proposed by Ling and Rivera.21 Frequency-domain-based analysis and synthesis of linear systems is a well-developed tool, and a similar approach for the nonlinear system analysis and synthesis is desirable. In this paper, a frequency-domainbased methodology for the CRMR of a class of nonlinear plants is proposed using the generalized frequency response functions of Peyton-Jones and Billings.25-27 In the proposed method, a higher order NARX model (designed as a plant proxy) is fitted to the appropriately obtained input/output open-loop data from the plant, following the classical orthogonal least-squares procedure of structure selection and parameter estimation, proposed by Chen et al.23 The suitability of simpler nonlinear structures to capture the relevant dynamics of the higher order plant is evaluated in the frequency domain. Specifically, we analyze the bilinear model, which has been studied quite extensively, as a candidate model structure for the reduced-order nonlinear model.15,28-30 The bilinear model is a special class of the NARX model, which is shown to capture the dynamics of many physical systems in a simple structure, especially in chemical engineering systems. Bilinear models can approximate the nonlinear processes to an accuracy restricted by the order of the model. They represent a mathematically tractable structure over the Volterra models for a nonlinear plant. Also, the bilinear model can obviously represent the dynamics of a nonlinear plant more accurately than a linear model. Hence, modeling and control of nonlinear systems in a bilinear framework have been extensively studied by Baheti et al.,31 Hsu and Mohler,32 Yeo and Williams,33,34 and Bartee and Georgakis.35 In this work, the CRMR and estimation is sought to be done by shaping the mismatch between the higher order plant proxy and the bilinear model structure in the frequency domain. Toward this objective, the controlrelevant prefiltering methodology proposed by Rivera et al.2 for linear systems is used here to prefilter the identification data. The design of the prefilter to reflect the specifications on the servo and regulatory control performance is done by appropriately minimizing the mismatch between the higher order plant proxy and the
bilinear model structure expressed as multiplicative uncertainty. The uncertainty is estimated in the frequency domain, using the generalized frequency response function for nonlinear systems as proposed by Peyton-Jones and Billings.25-27 The specification on the control performance and the nature of the set point/ disturbance signal that the system will be subjected to, therefore, define the control-relevant frequency band. For the prefilter as well, candidate linear and nonlinear structures are evaluated, and an algorithm for the synthesis of the control-relevant prefilter is proposed. The closed-loop performance of the nonlinear plant is evaluated using the bilinear model based NIMC scheme. The closed-loop simulation results and an analysis through the generalized frequency response are presented to validate the proposed methodology. The paper is structured as follows: In section 2, interpretation of the generalized frequency response functions defined for nonlinear systems by Peyton-Jones and Billings25-27 is introduced. The formulation of a CRMR problem is presented using the NIMC principles in section 3. The proposed CRMR methodology for a class of nonlinear systems is outlined in section 4, using the generalized frequency response function defined for the nonlinear systems. The closed-loop simulation results and the generalized frequency response characteristics of the estimated models, obtained by applying the proposed methodology to a simple illustrative example and to another representative example taken from the literature, are presented in section 5. 2. Generalized Frequency Response Functions The frequency domain description of linear systems is a well-established concept and has been used for the analysis and design of linear systems quite extensively. Describing functions represent a frequency domain representation of nonlinear systems which is derived from a quasi-linearization approach. The describing function is one-dimensional in frequency and is dependent on the input amplitude. This is defined for a given input type and, therefore, the describing function approach may not be applicable when the input type is different from the specific type for which it is defined or if the input contains two or more frequency components. Similarly, using the multidimensional Laplace transform and Fourier transform, nonlinear systems represented in the time domain by Volterra convolution integrals are mapped into multidimensional transfer functions and generalized frequency response functions, respectively.13,15 These multidimensional Volterra transfer functions incorporate the nonlinear frequency interaction through the multidimensional frequency components as its arguments and are not restricted to a particular input type. However, the multidimensionality in the representation makes the analysis and interpretation difficult and, therefore, use of these frequency response functions in the analysis of nonlinear systems is not quite useful for order beyond 2 or 3. The interpretation of the frequency space of increasing dimension and relating this higher dimension frequency space to input/output unidimensional spectra have been proposed by Peyton-Jones and Billings.26 They interpret the nonlinear multidimensional frequency response characteristics such as intermodulation, harmonic generation, and amplitude-dependent gain phase characteristics as being due to intrakernel and interkernel interferences, using a graphical analysis.27 Intrakernel
5008
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002
interference is nonlinear in frequency and describes the transfer of energy between spectral components within any of the Volterra kernels. The nonlinear, input amplitude-dependent, interkernel interference occurs between each of the Volterra kernels. The relationship between the multidimensional forms and the unidimensional system output can be described by these interference effects. Using the approach illustrated in ref 27, the multidimensional transfer functions are collapsed into one-dimensional form for any specific input type, in a specified frequency space, like an input frequency domain, an output frequency domain, or a constantdifference frequency domain. The evaluation of the one-dimensional generalized frequency response function for nonlinear systems having the NARX/bilinear structure can be carried out in the output frequency domain as follows. (i) Using the orthogonal least-squares method of Chen et al.,23 the NARX model of appropriate order and lag is estimated from the identification data. Using the recursive algorithm proposed by Peyton-Jones and Billings,25 the Volterra transfer functions H1(ejω1), H2(ejω1,ejω2), H3(ejω1,ejω2,ejω3), etc., are estimated, depending on the order of nonlinearity in the problem, in terms of the NARX model coefficients. (ii) The nth-order output frequency response can be evaluated by using the following expression: n
Yn(ejω1,...,ejωn) ) Hn(ejω1,...,ejωn)
U(ejω ) ∏ i)1
(1)
i
Here, U(ejω) represents the spectrum of the specified input. (iii) The one-dimensional frequency, as a function of which the generalized frequency response function has to be evaluated, can be defined in the output frequency domain as n
ω′n )
ωi w ωi ) ∑ i)1
{
ω′i - ω′i-1 1 < i e n ω′i i)1
(2)
(iv) The frequency response of the nth-order output can be defined in the one-dimensional output frequency domain as
Yn(ejω′n) )
1 (2π)n-1
∫-ππ‚‚‚∫-ππ Yn(ejω′ ,..., ej(ω′ -ω′ 1
n
n-1)
)
n-1
dω′1 ... dω′n-1 (3) (v) The total output response can be obtained by summing each of the nth-order output frequency responses, as follows: N
Y(ejω) )
Yn(ejω) ∑ n)1
(4)
This combination of (one-dimensional) nth-order output responses does not involve any further intermodulation or transfer of energy between frequencies. However, because the components in the output from different orders are vector quantities, the summation of eq 4 does introduce an amplitude-dependent interference effect between the kernels termed “interkernel interference”. (vi) If the input waveform is scaled by a constant factor Aw, then each of the nth-order outputs Yn(ejω),
being homogeneous, are scaled by the constant factor Anw. Therefore, eq 4 can be written as N
Y(Aw,ejω) )
Anw Yn(ejω) ∑ n)1
(5)
(vii) Hence, the generalized describing function, relating the input and output components at the same frequency, can be estimated as
N(Aw,ejω) )
Y(Aw,ejω) AwU(ejω)
(6)
Combining eqs 5 and 6 gives
N(Aw,ejω) )
N
1 jω
U(e )
jω An-1 ∑ w Yn(e ) n)1
(7)
(viii) Equations 2 and 7 together express the generalized frequency response function in the output frequency domain. The generalized frequency response function obtained following the above-specified procedure gives the onedimensional frequency response description of nonlinear systems. In the work done in this paper, the generalized frequency response function has been evaluated for the higher order NARX plant proxy and the reduced order bilinear model structures. Using these frequency response functions, the model-plant mismatch is estimated and an algorithm is proposed to estimate a control-relevant bilinear model which minimizes the mismatch in the chosen frequency range of interest, from the closed-loop performance viewpoint. This procedure is described in section 4 of the paper. 3. CRMR Problem Formulation As has been stated in the literature (see, for example, Van den Hof and Schrama36 and Rivera and Gaikwad4), a model exhibiting an excellent match with the true plant dynamics in an open-loop condition may not result in a good closed-loop control performance, using modelbased controllers. Thus, it is imperative that the closedloop performance requirements and the nature of set point/disturbance signals to which the closed-loop system will be subjected be incorporated in the estimation, to bias the estimation data. Rivera et al.1-4 proposed a prefilter-based control-relevant identification methodology for linear, single input single output (SISO)/multiinput multioutput (MIMO) systems, using appropriately obtained open-loop data. They have proposed the synthesis procedure of the prefilter required to be used for filtering the input/output data, so that the prefiltered data are biased to have the frequency content of importance from the closed-loop performance viewpoint. This has been achieved by equating the frequency domain expressions for minimization of the filtered prediction error with that of the control error minimization function. The prefilter expression thus obtained has been shown to explicitly include the sensitivity function, complimentary sensitivity function (defining the closedloop performance requirement), and the nature of the set point/disturbance signals that the closed-loop system will be subjected to. Also, the prefilter design is dependent on the model structure chosen and the control strategy to be adapted. The identification procedure has
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002 5009
where Ts is the sampling time and τcl is the specified rise time for the output response. The control error ec(t) is (cf. Figure 1)
ec(t) ) r(t) - y(t)
(10)
The objective function for the CRMR is to minimize the closed-loop control error using a control-relevant modelbased controller. Let J be the 2-norm of the control error defined as Figure 1. Single-loop nonlinear internal model control based system.
to be iterative because of the presence of the nominal model in the prefilter expression. Extension of these methodologies for the nonlinear systems is not straightforward because nonlinear systems are not amenable to being represented in the closed form in terms of transfer functions. For the nonlinear case as well, Ling and Rivera12,21 present an optimization-based approach to estimate the kernels of the reduced-order Volterra model. Their minimization procedure is decoupled with respect to each of the kernels. They propose the estimation of the linear kernel similar to the prefilter-based approach of linear systems. The successive higher order kernels are estimated by minimizing the control error contribution by the respective kernels. This minimization procedure makes use of only the subsequent lower order kernels, which are estimated in the minimization. In the following the CRMR problem is formulated for the nonlinear systems as a minimization of the weighted multiplicative error problem. It is shown that the weight function naturally influences the model reduction methodology in a control-relevant manner by explicitly incorporating the closed-loop performance requirements. In the formulation proposed here, we perform this minimization in the frequency domain using the analysis tools presented by Peyton-Jones and Billings25-27 for nonlinear systems. Furthermore, we restrict the analysis to bilinear structures for the model reduction; that is, we seek a reduced-order bilinear model. The CRMR problem for nonlinear systems is formulated as follows. Consider the NIMC-based single-loop system shown in Figure 1, where p(u,y) is the nonlinear plant to be controlled as per the desired specifications, p˜ (u,y˜ ) is the reduced-order model expressed in the bilinear model structure, q(e,u) is the bilinear model based nonlinear internal model controller, r(t) is the reference signal, ν(t) is the disturbance signal, e(t) is the controller input, and y(t) is the output from the system. The controller q(e,u) is synthesized as12
q(e,u) ) p˜ -1*Fimc(z)
J :) |ec(t)|2 ) lim
Nf∞
(12)
CRMR problem: Find min J′ ) min |ec(ejω)|2 (13) p˜ p˜ An expression for the control error represented in eq 10 can be obtained in the following manner. The output y(t) is expressed as
y(t) ) p(u,y) * q(e,u) * e(t) + ν(t)
(14)
Hence,
e(t) ) r(t) - y(t) + y˜ (t) ) r(t) - p(u,y) * q(e,u) * e(t) ν(t) + p˜ (u,y˜ ) * q(e,u) * e(t) (15) or
[1 + p(u,y) * q(e,u) - p˜ (u,y˜ ) * q(e,u)] * e(t) ) r(t) ν(t) (16) Using eqs 14 and 16, the control error ec(t) in eq 10 can be expressed as
ec(t) ) r(t) - p(u,y) * q(e,u) * [r(t) - ν(t)] - ν(t) (17) [1 + p(u,y) * q(e,u) - p˜ (u,y˜ ) * q(e,u)]
(
) 1-
)
[ (
p*q * [r(t) - ν(t)] 1 + p * q - p˜ * q
)
]
(18)
(1 + p * q - p˜ * q) - p * q * [r(t) - ν(t)] 1 + p * q - p˜ * q
(19)
1 - p˜ * q * [r(t) - ν(t)] 1 + (p - p˜ ) * q
(20)
)
) (1 - p˜ * q) * [1 + (p -p˜ ) * q]-1 * [r(t) - ν(t)] (21) Under the assumption that |(p - p˜ ) * q| , 1 (cf. remark 1), the Taylor series expansion of [1 + (p - p˜ ) * q]-1 in eq 21 gives
[1 + (p - p˜ ) * q]-1 ≈ [1 - (p - p˜ ) * q] (9)
(11)
Then, the CRMR problem in the frequency domain is stated as the problem of minimization of J′
where * represents an appropriate composition/convolution operation to obtain a composite representation of the dynamics and Fimc(z) is the linear prefilter required to achieve robustness and map the error e (controller input) into the output domain of the model (i.e., into the input range of the model inverse). As in earlier approaches,12,21 an IMC filter Fimc(z) is chosen based on the desired specifications as
z - e-Ts/τcl
1/2
J′ :) |ec(ejω)|2, -π e ω e π
)
Fimc(z) )
]
ec(t)2 ∑ N t)0
where N is the number of data points considered in the identification. In the frequency domain, J is equivalent to J′ defined as
(8)
(1 - e-Ts/τcl)z
[
1 t)N
and eq 21 can then be approximately written as
(22)
5010
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002
ec(t) ≈ (1 - p˜ * q) * [1 - (p - p˜ ) * q] * [r(t) - ν(t)] (23) Assuming all signals in eq 23 are square integrable and quasi-stationary, Parseval’s theorem37 can be applied to give
[2π1 ∫ |1 - p˜ (e π
||ec(ejω)||2 ≈
jω
-π
||
) * q(ejω) 2 1 - [p(ejω) -
||
| ]
p˜ (ejω)] * q(ejω) 2 r(ejω) - ν(ejω) 2 dω
1/2
(24)
p(ejω) - p˜ (ejω)
∫-ππ|˜ (ejω)|2|em(ejω)η˜ (ejω)|2|r(ejω) -
| )
ν(ejω) 2 dω
(25)
jω
p˜ (e )
where p(ejω) and p˜ (ejω) are the generalized frequency response functions of the plant and the reduced-order bilinear model, respectively. Now, the frequency response dynamics of the plant is estimated as that of a proxy plant pˆ obtained by fitting a higher order NARX/ NARMAX model to the identification data. The multiplicative uncertainty term em(ejω) can, therefore, be expressed in terms of this proxy plant as
em(ejω) )
CRMR problem: 1 Find min 2π p˜
(
Define the multiplicative error em(ejω) as
em(ejω) :)
tion (which is inevitable even when p˜ ) p). Therefore, the effect of this term on the control error minimization, and hence on the CRMR problem, is neglected. Therefore, the CRMR problem reduces to minimizing the contribution of only the second term on the right-hand side of the above inequality, which explicitly contains the model-plant mismatch description that affects the control error.2 Hence, the CRMR problem can be defined as
pˆ (ejω) - p˜ (ejω)
(26)
p˜ (ejω)
Further, from eq (8),
Fimc(ejω) ) p˜ (ejω) * q(ejω)
(27)
1/2
(32)
Minimization of the CRMR problem stated above minimizes the second component of inequality (31), through minimizing the model-plant mismatch in the desired frequency range of importance. Define a weight function W(ejω) as
|W(ejω)| :) |˜ (ejω) η˜ (ejω) [r(ejω) - ν(ejω)]|
(33)
where ˜ (ejω) and η˜(ejω) are the performance specifications for the overall closed-loop system and r(ejω) and ν(ejω) are the frequency domain descriptions of the set-point and disturbance signals for which the closed-loop system is designed. Because this weight function characterizes the desired closed-loop performance specifications, it can be considered as a control-relevant weight function. Then, the (reduced) CRMR problem in eq 32 can be written as
Using eqs 26 and 27 in eq 24 gives
|ec(ejω)|2 ≈
1 2π
(
∫
π
|
CRMR problem:
||
1 - Fimc(ejω) 2 1 -π
||
| )
em(ejω) Fimc(ejω) 2 r(ejω) - ν(ejω) 2 dω By the properties of IMC
˜ (ejω) ) 1 - Fimc(ejω)
1 Find min 2π p˜
(
1/2
(28)
synthesis,38
η˜ (ejω) ) Fimc(ejω) (29)
and
Using eq 29 in eq 28 gives
1 |ec(e )|2 ≈ 2π
(
jω
∫-π|˜ (e π
||
jω 2
jω
||
jω 2
jω
) 1 - em(e ) η˜ (e ) r(e ) -
| )
ν(ejω) 2 dω
1/2
(30)
By the triangular inequality,39 the above gives
|ec(ejω)|2 e
(2π1 ∫ |˜(e )| |r(e π
jω 2
-π
(2π1 ∫ |˜(e )| |e (e π
-π
jω 2
m
jω
jω
| )
) - ν(ejω) 2 dω
||
1/2
+
| )
) η˜ (ejω) 2 r(ejω) - ν(ejω) 2 dω
1/ 2
(31)
Now, the CRMR problem in eq 13 involves minimization of ||ec(ejω)||2, given by the above inequality. However, the first term on the right-hand side of the above inequality is based on the nominal properties of the closed-loop performance (assuming that p˜ ) p), and therefore its contribution to the control error is due to the performance degradation from perfect control. The IMC filter, which is necessary to achieve robustness in performance, contributes to this performance degrada-
∫-ππ|W(ejω) em(ejω)|2 dω)
1/2
(34)
Replacing the integral in eq 34 by an equivalent sum over the desired frequency range of importance from the closed-loop performance point of view gives
CRMR problem: 1 i)K W(ejωi) em(ejωi) Find min p˜ K i)1
(
∑|
|)
1/2
2
, -π e ωi e π (35)
where K is the number of discrete frequency points chosen. The minimization of eq 35 yields a reducedorder model that captures the plant dynamics in the frequency range, important from a closed-loop performance viewpoint. This is seen from the fact that the weighting term associated with the multiplicative uncertainty explicitly contains the closed-loop performance specifications, in terms of τcl, and set point/disturbance characteristics, in terms of r - ν. In the proposed approach, a prefilter for identification data is sought to be designed, via an explicit frequency domain minimization given by eq 35 in the space of prefilter parameters. For analysis, we consider filters with linear as well as nonlinear structures (for example, bilinear structure) that could aid in shaping the modelplant mismatch in the desired frequency range. A control-relevant bilinear model is fitted to the prefiltered data, and using the identified bilinear model, the closedloop performance is analyzed in an NIMC scheme.
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002 5011
Remark 1. The derivation of the CRMR problem is based on the assumption made in eq 21, i.e., that |(p p˜ ) * q| , 1 holds over the frequency range of importance from a closed-loop performance viewpoint. Now, from eqs 8, 25, and 29, we have q ) p˜ -1Fimc, η˜ ) Fimc, and em ) p - p˜ /p˜ . Hence,
(p - p˜ ) * q ) (p - p˜ ) * p˜ -1 * p˜ * q ) p - p˜ * p˜ * p˜ -1Fimc ) emFimc ) emη˜ p˜
(
)
Therefore, the assumption |(p - p˜ )q| , 1 is equivalent to the assumption |emη˜ | , 1. Rivera et al.2 show the assumption to be valid in the bandwidth of |˜ (r - ν)|, i.e., in the control-relevant frequency range. 4. Algorithm for CRMR In this section, the algorithm for control-relevant bilinear model identification proposed in section 3 is outlined. The closed-loop system is assumed to be subjected to step set point/disturbance signals. The first step is to synthesize the weight function. Using eq 29, the frequency domain representations of η˜ (ejω) and ˜ (ejω) are estimated. The frequency domain representation of the set point/disturbance signal r(ejω)/ν(ejω) is obtained based on the nature of the signal. Using these, the weight function is estimated as defined in eq 33 and is represented as a vector W(ejω). The weight function shapes the model-plant mismatch to be minimum in the desired band of frequency important from the closedloop performance perspective because it explicitly incorporates the closed-loop performance requirement and the nature of the set point/disturbance signal that the closed-loop system will be subjected to. The CRMR can then be performed by following the steps listed below. 1. The nonlinear plant to be modeled is perturbed about the given operating point by an appropriately designed perturbation signal, and the open-loop input/ output data are obtained. In the analysis following here, a multilevel pseudo random signal of appropriate variation and level, suitably designed from a priori knowledge of the plant, is used as the perturbation signal. 2. Using the input/output data set, various possible combinations of input, output, and product of input/ output having different lag elements are generated. Based on the classical orthogonal least-squares method proposed by Chen et al.,23 a NARX model of optimal order and lag is identified using this data set. This NARX model is considered as an equivalent or proxy to the true plant [pˆ (y,u) ≈ p(y,u)] in the model reduction process. 3. The generalized frequency response function of the identified NARX model is obtained by following the procedure given in refs 25-27, 40, and 41. The frequency response is obtained in the output frequency domain and is recorded as a vector pˆ (ejω). 4. The open-loop input/output data are appropriately prefiltered by a prefilter having a linear/bilinear structure, obtained by initialization or from previous iteration. Let the coefficients representing the prefilter structure be represented as a vector θpf. 5. Using the prefiltered data, a first-order bilinear model is identified by assuming an multiinput singleoutput (MISO) autoregressive with an exogenous (ARX) inputlike structure.
6. The generalized frequency response function of the bilinear model identified is obtained and recorded as a vector p˜ (ejω). 7. The multiplicative error between the NARX model and the bilinear model identified is estimated as specified in eq 26 and is recorded as a vector em. 8. The objective function of identification, i.e., the weighted multiplicative error, is estimated in the frequency domain as
J′ )
(
∑ |W(ejω ) em(ejω )|2 K i)1 1 i)K
i
i
)
1/2
, -π e ωi e π
(36)
9. An optimization algorithm is used which minimizes the weighted multiplicative error J′ to bring it within the specified tolerance limit in the space of prefilter parameters θpf. This will involve the iterations of steps 4-8 until the prefilter parameters θpf converges. 10. The identification data prefiltered by the prefilter synthesized in the last iteration when used to fit a bilinear model is expected to give a control-relevant bilinear model. Using this control-relevant bilinear model, a NIMC-based controller is synthesized in the following manner. The bilinear model p(yd,ud) identified is represented as
yd(k+1) ) amyd(k) + bmud(k) + cmyd(k) ud(k) (37) Equation 37 can be rearranged as
ud(k) )
yd(k+1) - amyd(k) bm + cmyd(k)
(38)
Equation 38 is a formal representation of the inverse of a bilinear system, which is required for the implementation of a NIMC strategy based closed-loop system. However, it is to be noted here that the controller (model inverse) output ud(k) is dependent on the next sample of the controller input, that is, e(k+1) (replacing output yd with error e in eq 38). When ud(k+1) is substituted for ud(k) on the left-hand side of eq 38, the estimation of inverse becomes causal. This is similar to partitioning of the linear model into all-pass and minimum-phase elements and writing the linear model inverse as the inverse of the all-pass element. The above procedure for arriving at a noncausal structure for the inverse of nonlinear models is similar to the one given by Morari and Zafiriou38 for linear models. Therefore, the model inverse can be defined as
ud(k+1) )
e(k+1) - ame(k) bm + cme(k)
(39)
where ud is the output of the controller and e is the input to the controller. This model inverse given by eq 39 is augmented by the IMC filter, Fimc(z), as defined in eq 9. The linear filter Fimc(z) is introduced for robustness and realizability. Also, when the error e does not belong to the output space of the model (which is input space for the model inverse), it is assumed that Fimc(z) will project e into the appropriate space before it is applied to the input of the model inverse. Using the aboveproposed bilinear, model-based, controller, the true plant performance is analyzed. 5. Simulation Results The control-relevant bilinear model reduction methodology proposed in this paper is validated by applica-
5012
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002
tion to representative nonlinear plants. In this section the results obtained from simulation studies are presented. 5.1. Illustrative Example 1. Consider a system described by the following difference equation:
y(k+3) ) 0.6y(k+2) + 0.5y(k+1) - 0.3y(k) + 2u(k+ 2) + 0.5u(k+1) + 0.6y(k+2) u(k+2) + 0.2y(k+1) u(k+1) (40) This system is already in NARX structure. The objective here is to identify a first-order bilinear model in a control-relevant manner for this plant. The performance requirement from the closed-loop system is specified to have an 8 s rise time. The sampling time chosen is of 1 s. This plant is perturbed by a three-level pseudo random signal generated with a shift register length of 7, having signal levels of (+1, 0, -1) in the presence of independent and identically distributed (iid) random noise of variance 0.1 acting at the output, and the input/ output data are recorded. Using these data, a new set of data is generated with a product of input and output vectors, which is used as a second input in fitting the bilinear element. A MISO ARX routine is used to fit the bilinear model. Using this set of data, a bilinear model is identified which can be written as
yd(k+1) ) 0.9029yd(k) + 2.0691ud(k) + 0.5902yd(k) ud(k) (41) A linear structure for the prefilter is assumed and is initialized as
L1(z) )
0.2z-1 + 0.1z-2 + 0.1z-3 + 0.1z-4 1 - 0.3z-1 - 0.1z-2 - 0.1z-3 - 0.1z-4
(42)
The optimization routine minimizing the weighted multiplicative error in frequency domain resulted in a linear prefilter given by
L1o ) 0.2740z-1 - 0.0794z-2 + 0.3275z-3 + 0.0283z-4 1 - 0.0259z-1 - 0.1475z-2 + 0.2033z-3 - 0.1268z-4 (43)
yf(k+1) ) 0.4366yf(k) - 2.9423y(k) - 0.028yf(k) y(k) (44) Similarly, the input and the product of input/output are also filtered. The control-relevant bilinear model identified using the prefiltered data set can be written as
yc(k+1) ) 0.8579yc(k) + 1.8393uc(k) + 0.1502yc(k) uc(k) (45) The frequency response characteristics of the given NARX plant, bilinear model identified without prefiltering, the control-relevant bilinear model identified with linear prefiltering, and the closed-loop performance specification function [Fimc(z)] are shown in Figure 2. A similar plot comparing the frequency response characteristics of the control-relevant bilinear model identified with a prefilter having bilinear structure is shown in Figure 3. The closed-loop responses obtained using a direct identified bilinear model and the two controlrelevant bilinear models using NIMC strategy are shown in Figure 4, with a set point of unity and a random noise of zero mean and variance 0.01 acting at the output. From Figures 2 and 3, it can be seen that the frequency response characteristics of the two controlrelevant models exhibit minimum mismatch with that of the true plant in the bandwidth of the closed-loop specification function. However, the linear prefilterbased control-relevant bilinear model matches the true plant in a slightly higher band of frequency. From Figure 4, it can be seen that the closed-loop response obtained from the linear prefilter-based control-relevant bilinear model exhibits marginally higher settling time and rise time, while the response based on the controlrelevant bilinear model identified using a prefilter having bilinear structure is superior as per specification. However, the direct identified nominal model-based system performance is very poor compared to those obtained from the two control-relevant models. 5.2. Continuous Stirred Tank Reactor (CSTR) Problem. As a second illustration, the CSTR described in refs 12, 18, and 21 is considered. The description of the reactor is given as follows:
x˘ 1 ) 10(6 - x1) - 2.4568x1xx2 x˘ 2 ) 80u - 10.1022x2
The time taken for the linear prefilter coefficients to converge was 278.64 units. This identification data prefiltered by L1o resulted in a control-relevant bilinear model expressed as
x˘ 3 ) 0.002412x1xx2 + 0.11219x2 - 10x3
ycl(k+1) ) 0.7653ycl(k) + 2.8354ucl(k) + 1.1163ycl(k) ucl(k)
y ) x4/x3
A similar analysis is carried out with a prefilter in bilinear structure. The prefilter coefficients are initialized as af ) 0.3, bf ) -0.1, and cf ) -0.05. Using the algorithm proposed in section 4, the prefilter in bilinear structure is synthesized to have coefficients of [0.4366, -2.9423, -0.028]. The time taken for the convergence of the coefficients of the prefilter having bilinear structure was 85.41 units. The prefilter equation, to obtain a filtered signal [yf(k)] using the prefilter having bilinear structure, for the signal [y(k)] is given by
x˘ 4 ) 245.978x1xx2 - 10x4 (46)
Identification data are obtained by perturbing the plant by a five-level pseudo random signal, of levels [-0.0002, -0.0001, 0, 0.0001, 0.0002], generated with a shift register of length 5. A random noise of zero mean and variance 10-8 was also present during data collection. Nominal operating points for the system are given as x10 ) 5.506 77, x20 ) 0.132 906, x30 ) 0.001 975 2, x40 ) 49.3818, u0 ) 0.016 78 m3 h-1, and y0 ) 25 000.5 NAMW. State variables xi, input u, and output y are first normalized using zi ) (xi - xi0)/xi0, ν ) (u - u0)/u0, and w ) (y - y0)/y0. Using this normalized data, various
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002 5013
Figure 2. Generalized frequency response function comparison: solid line, true plant; dashed line, nominal bilinear model; dots, control-relevant bilinear model based on linear prefiltering; dashdotted line, nominal complementary sensitivity function.
Figure 4. Closed-loop performance comparison: solid line, using a nominal bilinear model; dots, using a control-relevant bilinear model estimated using a linear prefilter; dashed line, using a control-relevant bilinear model estimated using a prefilter having bilinear structure.
With this prefilter, the identification data are prefiltered and a bilinear model is identified. The model-plant mismatch is then estimated, and it is weighted by the designed weight function, based on the specification and step set point/disturbance signal. The linear prefilter synthesized by minimizing the weighted model-plant mismatch is given as
L2o(z) )
-0.0658z-1 - 0.4256z-2 1 - 1.9151z-1 - 0.9977z-2
The time taken for convergence of the linear prefilter coefficients was 826.13 units. The control-relevant bilinear model identified using the prefilter L2o can be written as Figure 3. Generalized frequency response function comparison: solid line, true plant; dashed line, nominal bilinear model; dots, control-relevant bilinear model obtained using a prefilter having bilinear structure; dash-dotted line, nominal complementary sensitivity function.
combinations of product of input/output and the two are obtained with different lags. Using this normalized input/output data, a higher order NARX model is identified using the classical orthogonal least-squares method proposed by Chen et al.23 The NARX model identified is considered as the plant proxy in further analysis. The identified NARX model is represented as follows:
y(k+2) ) 1.235y(k+1) - 0.382y(k) - 0.0348u(k+ 1) - 0.0253u(k) - 0.0305u(k+1) u(k+1) 0.0266u(k) u(k) (47) The closed-loop system is specified to have a rise time of 0.3 h. The sampling time chosen for this problem is 0.05 h. To obtain a reduced-order model of this plant proxy, the method proposed in section 4 is adapted. A linear prefilter having the following structure is initialized.
L2(z) )
-0.3z-1 - 0.4z-2 1 - 2z-1 + z-2
ycl(k+1) ) 0.8541ycl(k) - 0.0909ucl(k) 0.0011ycl(k) ucl(k) The control-relevant bilinear model identification is also carried out with a prefilter having bilinear structure. The prefilter in bilinear structure is initialized, with coefficients af ) -0.4, bf ) -0.2, and cf ) -0.1, having the filter structure yf(k) ) afyf(k-1) + bfy(k-1) + cfy(k-1) yf(k-1), where yf is the prefiltered signal of y obtained using a prefilter having bilinear structure. The optimal prefilter coefficients obtained are af ) 0.6045, bf ) -0.2786, and cf ) -0.0539. In this case the time taken for convergence of the coefficients of the prefilter having bilinear structure was 229.09 units. With these prefilter coefficients, the normalized input, output, and their product are prefiltered. A controlrelevant first-order bilinear model is fitted to the prefiltered data and is obtained as
yc(k+1) ) 0.8894yc(k) - 0.0581uc(k) 0.0291yc(k) uc(k) (48) A nominal model is identified by fitting to the identification data without prefiltering and can be written as
yd(k+1) ) 0.9261yd(k) - 0.0341ud(k) 0.4057yd(k) ud(k) (49)
5014
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002
Figure 5. Generalized frequency response comparison: solid line, true plant; dashed line, nominal bilinear model; dots, controlrelevant bilinear model estimated using a linear prefilter; dashdotted line, closed-loop specification function.
Figure 6. Generalized frequency response comparison: solid line, true plant; dashed line, nominal bilinear model; dots, controlrelevant bilinear model estimated using a prefilter having bilinear structure; dash-dotted line, closed-loop specification function.
The generalized frequency response functions of the plant proxy (equivalent to the true plant), bilinear models obtained with linear prefilter and without prefiltering are plotted in Figure 5, along with the frequency response characteristics of the desired closedloop behavior. Also, a similar plot comparing the frequency domain characteristics of the control-relevant model estimated using a prefilter having bilinear structure is shown in Figure 6. From Figures 5 and 6, it can be shown that the control-relevant bilinear models exhibit minimum mismatch with the NARX plant proxy in the frequency range of importance, compared to the mismatch of the direct identified bilinear model with the plant proxy. It can also be shown that the linear prefilter-based control-relevant model shows a better match with the frequency domain characteristics of the plant proxy than that obtained with the control-relevant model identified with a prefilter having bilinear structure. The closed-loop responses obtained using these two control-relevant bilinear model-based systems are compared with that obtained from the direct identified bilinear model-based system. These are shown in Figure 7 (normalized output). In performing the closed-loop simulation, a set point of 0.1 (normalized output) for the output is specified and zero mean random noise of
Figure 7. Closed-loop tracking performance comparison: solid line, based on a nominal bilinear model; dots, based on a controlrelevant bilinear model estimated using a linear prefilter; dashed line, based on a control-relevant bilinear model estimated using a prefilter having bilinear structure.
Figure 8. Closed-loop regulation performance comparison: solid line, based on a nominal bilinear model; dots, based on a controlrelevant bilinear model obtained using a linear prefilter; dashed line, based on a control-relevant model obtained using a prefilter having bilinear structure.
variance 10-8 is considered to be acting at the output of the system. From Figure 7, it can be shown that the bilinear control-relevant model identified using a prefilter in bilinear structure exhibits a marginally superior performance than that obtained from the controlrelevant bilinear model identified using a linear prefilter. The step disturbance regulation performance obtained is shown in Figure 8. From Figure 8 also, it may be seen that the control-relevant model identified using a prefilter in bilinear structure results in a marginally superior performance than that obtained using the control-relevant model identified based on a linear prefilter. Even in this case, the direct identified nominal model-based system performance is shown to be poor. 6. Conclusions A prefilter-based approach for the control-relevant bilinear model identification is proposed. A prefilter having linear/bilinear structure is synthesized to minimize the model-plant mismatch, expressed in terms of the generalized frequency response functions of the plant (or its proxy) and the bilinear model. The model-
Ind. Eng. Chem. Res., Vol. 41, No. 20, 2002 5015
plant mismatch is weighted by the closed-loop performance requirement to bias the identification data from a closed-loop performance point of view. Using the bilinear, model-based, NIMC strategy, the closed-loop simulation results are obtained, which validate the proposed approach. A prefilter in bilinear structure having only three parameters is found to achieve the desired minimization compared to the large number of parameters required for a linear prefilter. Also, with respect to computation time, it is observed that the prefilter synthesis in bilinear structure is much faster than the linear prefilter synthesis. The identification and control methodology proposed here is applicable, in the presence of iid white noise added at the output, during identification data collection and also during the closed-loop simulation. However, issues related to controlrelevant identification in the presence of linear/nonlinear noise affecting the output are potential areas of future investigation. Acknowledgment C.S. expresses his gratitude to the management of NMAM Institute of Technology, Nitte, Udupi, Karnataka, India, for sponsoring him in the Ph.D. program at IIT Bombay, India, and the financial assistance. Literature Cited (1) Rivera, D. E. Control Relevant Parameter Estimation: A Systematic Procedure for Pre-filter Design. Proceedings of the 1991 American Control Conference; American Automatic Control Council (AACC): Evanston, IL, 1991; pp 237-241. (2) Rivera, D. E.; Pollard, J. F.; Garcia, C. E. Control-Relevant Pre-filtering: A Systematic Design Approach and Case Study. IEEE Trans. Autom. Control 1992, 37 (7), 964-974. (3) Rivera, D. E.; Gaikwad, S. V. Modeling for Control Design in Combined Feedback/Feed forward Control. Proceedings of the American Control Conference; American Automatic Control Council (AACC): Evanston, IL, 1992; pp 1445-1446. (4) Rivera, D. E.; Gaikwad, S. V. Systematic Techniques for Determining Modeling Requirements for SISO and MIMO Feedback Control. J. Process Control 1995, 5 (4), 213-224. (5) Shook, D. S.; Mohtadi, C.; Shah, S. L. A Control-Relevant Identification Strategy for GPC. IEEE Trans. Autom. Control 1992, 37 (7), 975-980. (6) Kwok, K. Y.; Shah, S. L. Long-Range Predictive Control with a Terminal Matching Condition. Chem. Eng. Sci. 1994, 49 (9), 1287-1300. (7) Hjalmarsson, H.; Gevers, M.; Bruyne, F. D. For Model Based Control Design, Closed Loop Identification Gives Better Performance. Automatica 1996, 32, 1659-1673. (8) Narendra, K. S.; Gallman, P. G. An Iterative Method for Identification of the Nonlinear Systems Using the Hammerstein Models. IEEE Trans. Autom. Control 1966, 12, 546. (9) Eskinat, E.; Johnson, S. H. Use of Hammerstein Models in Identification of Nonlinear Systems. AIChE J. 1991, 37 (2), 255268. (10) Pearson, R. K. Nonlinear Input/Output modeling. J. Process Control 1995, 5 (4), 197-211. (11) Bai, E. W. An Optimal Two-Stage Identification Algorithm for Hammerstein-Weiner Nonlinear Systems. Automatica 1998, 34 (3), 333-338. (12) Ling, W. M.; Rivera, D. E. A Methodology for ControlRelevant Nonlinear System Identification using Restricted Complexity Models. J. Process Control 2001, 11 (2), 209-222. (13) Schetzen, M. Volterra and Wiener Theories of Non-Linear Systems; Wiley: New York, 1980. (14) Rugh, W. J. Non-Linear System Theory; The John Hopkins University Press: Baltimore, MD, 1981. (15) Mathews, V. J.; Sicuranza, G. L. Polynomial Signal Processing; John Willey & Sons: New York, 2000.
(16) Barker, H. A.; Pradisthayon, T. Performance of Antisymmetric Pseudo Random Signals in the Measurement of 2nd Order Volterra Kernels by Cross Correlation. IEE Proc. 1972, 119 (3), 353-362. (17) Barker, H. A.; Davy, R. W. Measurement of Second-Order Volterra Kernels using Pseudo Random Ternary Signals. Int. J. Control 1978, 27 (2), 277-291. (18) Doyle, F. J., III; Ogunaike, B. A.; Pearson, R. K. Nonlinear Model-based Control using Second-Order Volterra Models. Automatica 1995, 31 (5), 697-714. (19) Pearson, R. K.; Ogunaike, B. A.; Doyle, F. J., III. Identification of Structurally Constrained Second-Order Volterra Models. IEEE Trans. Signal Process. 1996, 44 (11), 2837-2846. (20) Maner, B. R.; Doyle, F. J., III; Ogunaike, B. A.; Pearson, R. K. Nonlinear Model Predictive Control of a Simulated Multi Variable Polymerization Reactor using Second-Order Volterra Models. Automatica 1996, 32, 1285. (21) Ling, W. M.; Rivera, D. E. Control Relevant Model Reduction of Volterra Series Models. J. Process Control 1998, 8 (2), 7988. (22) Zheng, Q.; Zafiriou, E. Control-Relevant Identification of Volterra series Models. Proceedings of American Control Conference; American Automatic Control Council (AACC): Evanston, IL, 1994; p 2050. (23) Chen, S.; Billings, S. A.; Luo, W. Orthogonal Least Squares Methods and Their Application to Non-Linear System Identification. Int. J. Control 1989, 1873-1896. (24) Diaz, H.; Desrochers, A. A. Modeling of Nonlinear Discrete Time Systems from Input/Output Data. Automatica 1988, 24, 629. (25) Peyton-Jones, J. C.; Billings, S. A. Recursive Algorithm for Computing the Frequency Response of a Class of Non-Linear Difference Equation Models. Int. J. Control 1989, 50 (5), 19251940. (26) Peyton-Jones, J. C.; Billings, S. A. Interpretation of NonLinear Frequency Response Functions. Int. J. Control 1990, 52 (2), 319-346. (27) Peyton-Jones, J. C.; Billings, S. A. Describing Functions, Volterra Series, and the Analysis of Nonlinear Systems in the Frequency Domain. Int. J. Control 1991, 53 (4), 871-887. (28) Espana, L. Reduced Order Bilinear Models for Distillation Columns. Automatica 1978, 14, 345. (29) Mohler, R. R. Nonlinear Systems Volume 1, Dynamics and Control; Prentice-Hall: Englewood Cliffs, NJ, 1991. (30) Mohler, R. R. Nonlinear Systems Volume 2, Application to Bilinear Control; Prentice-Hall: Englewood Cliffs, NJ, 1991. (31) Baheti, R. S.; Mohler, R. R.; Spang, H. A. Second-Order Correlation Methods for Bilinear System Identification. IEEE Trans. Autom. Control 1980, 25 (6), 1141-1146. (32) Hsu, C. S.; Mohler, R. R. On the Inverse of a Special Class of Bilinear Systems. J. Dyn. Syst., Meas., Control 1981, 102, 103105. (33) Yeo, Y. K.; Williams, D. C. Adaptive Bilinear Model Predictive Control. Proc. Am. Control Conf. 1986, 1455-1461. (34) Yeo, Y. K.; Williams, D. C. Bilinear Model Predictive Control. Ind. Eng. Chem. Process Des. Dev. 1987, 26, 2267-2274. (35) Bartee, J. F.; Georgakis, C. Identification and Control of Bilinear Systems. Proc. Am. Control Conf. 1992, 2576-2580. (36) Van den Hof, P. M. J.; Schrama, R. J. P. Identification and ControlsClosed-loop Issues. Automatica 1995, 31 (12), 1751-1770. (37) Ljung, L. System Identification: Theory for the User; Prentice-Hall: Englewood Cliffs, NJ, 1987. (38) Morari, M.; Zafiriou, E. Robust Process Control; PrenticeHall: Englewood Cliffs: NJ, 1989. (39) Skogestad, S.; Postlethwaite, I. Multi Variable Feedback Control: Analysis and Design; John Wiley & Sons: New York, 1996. (40) Nataraj, P. S. V.; Date, P.; Umrani, A. Robust Feedback Synthesis for Nonlinear Integro Differential Equation Models Using Generalized Describing Functions. Automatica 1997, 33 (5), 959-962. (41) Kotecha, K. V. Design of Robust Controller for Nonlinear Plant Using Generalized Describing Function. M. Tech. Dissertation, IIT Bombay, Mumbai, India, 2000.
Received for review November 26, 2001 Revised manuscript received May 23, 2002 Accepted July 23, 2002 IE010946X
