Automatica 40 (2004) 1891 – 1903 www.elsevier.com/locate/automatica
Two-dimensional frequency analysis for unconstrained model predictive control of cross-directional processes夡 Junqiang Fana,∗ , Gregory E. Stewarta , Guy A. Dumontb a Honeywell Automation and Control Solutions, 500 Brooksbank Avenue, North Vancouver, BC, Canada, V7J 3S4 b Department of Electrical and Computer Engineering, University of British Columbia, 2356 Main Mall, Vancouver, BC, Canada V6T 1Z4
Received 27 March 2003; received in revised form 19 May 2004; accepted 24 June 2004
Abstract This paper presents consistent criteria for evaluating the selection of tuning parameters for an industrial model predictive control of large-scale cross-directional (CD) processes using a two-dimensional (temporal and spatial) frequency analysis technique. The concept of rectangular circulant matrices (RCMs) and their properties are presented. It is shown that large-scale CD processes can be approximated as RCMs and then diagonalized by complex Fourier matrices, allowing analysis in terms of a family of SISO transfer functions across the spatial frequencies. Familiar concepts from control engineering such as bandwidth and stability margin are extended into the twodimensional frequency domain, providing intuitive measures of closed-loop performance and robustness. 䉷 2004 Elsevier Ltd. All rights reserved. Keywords: Frequency domain; Matrix methods; Model predictive control; Model uncertainty; Paper machine; Performance analysis; Robust control; Sensitivity function; Stability analysis; Two-dimensional systems
1. Introduction A typical paper machine is shown in Fig. 1. For profile control, it is equipped with a number of cross-directional (CD) actuator arrays and scanning sensor(s). As illustrated, the direction of sheet travel is known as the machine direction (MD) while the direction perpendicular to the sheet travel is known as the cross direction (CD). The MD represents the temporal dimension while the CD represents the spatial dimension. The most important controlled properties of paper are weight, moisture, and caliper (thickness). The quality of the paper sheet is defined in terms of the variability of these properties. Details describing the functioning of actuator arrays (such as slice lip, etc. in Fig. 1) may be found in Cutshall (1991). The control objective of these actuators is to 夡 This paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor Tor Arne Johansen under the direction of Editor Frank Allgöwer. ∗ Corresponding author. Tel.: +1-604-982-3634; fax: +1-604-980-0120. E-mail address:
[email protected] (J. Fan).
0005-1098/$ - see front matter 䉷 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2004.07.002
minimize the variability of the measured properties of the paper sheet (Heath & Wellstead, 1995). In industrial paper machine CD processes, the measurements generally outnumber the actuators (typically 3–10 times). The resulting rectangular process model transfer matrices are usually transformed into square plants by first filtering the high-resolution outputs with a low-pass antialiasing spatial filter and then downsampling the outputs to the resolution of the inputs. Because the actuators in an array are identically constructed and evenly distributed across the sheet, the responses to these actuators are modeled as identical except for those near the edges of the sheet. In other words, these systems are spatially invariant except near the spatial boundaries (Duncan & Bryant, 1997; Laughlin, Morari, & Braatz, 1993; Mijanovic, 2004; Stewart, Gorinevsky, & Dumont, 2003a). Based on a spatially invariant approximation, the spatial frequency concept for web processes was thoroughly addressed by Duncan (1989). Spatial frequencies were also proposed for application in CD control for web forming processes by Laughlin et al. (1993). The spatial bandwidth
1892
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
Process response [gsm]
Fig. 1. Wide view of the paper machine showing typical positions of the various actuator arrays and scanning sensor(s). (Artwork courtesy of Honeywell).
0.01 0.005 50
0 0.005 0.01
100 0
50
100
150
200
250
300 150
Scanner measurement locations Actuator position [microns]
2 200 1
-1 300 -2
(a)
250
0
10
20
30
40
Actuator number
50
0 20 40 60
60
(b)
Fig. 2. (a) Nominal steady-state response of basis weight to slice lip actuators; (b) non-zero elements of the spatial interaction matrix G0 in (2).
of CD control systems for web processes was investigated by Duncan and Bryant (1997), where it was illustrated that the spatial bandwidth of the whole array of actuators depends on the spatial frequency response of a single actuator. A two-dimensional frequency domain loop shaping design technique was developed in Stewart, Gorinevsky, and Dumont (2003b) and applied to CD controller design in Stewart, Gorinevsky, and Dumont (2003a). More recently, a spatial frequency anti-windup strategy for cross-directional control problems was provided by Rojas, Goodwin, Desbens, and Johnston (2002a) and Rojas, Goodwin, and Johnston (2002b). If the responses of the actuators near the sheet edges (see Fig. 2a for example) are approximated as being circularly shifted from the central actuator’s response1 and the measurement array is processed to be the same resolution as the inputs using the method mentioned above, then these square plants are defined in terms of circulant matri1 There exist certain plastic film processes (Featherstone, VanAntwerp, & Braatz, 2000; Laughlin et al., 1993) that have this circularly shifting property as they are extruded in a tube.
ces (Hovd, Braatz, & Skogestad, 1997; Hovd & Skogestad, 1994; Laughlin et al., 1993; Mijanovic, 2004; Stewart et al., 2003b), which can be diagonalized by real Fourier matrices (Davis, 1994; Gray, 2002). However, the use of highresolution measurements is desirable from the papermakers’ point of view as paper quality is evaluated at the highest resolution available. In addition, it is generally difficult to square the process model when considering CD control with multiple actuator arrays and multiple controlled properties (Backström, Gheorghe, Stewart, & Vyse, 2001). In this paper, the rectangular CD process model matrix will be approximated as spatially invariant. Rectangular circulant matrices (RCMs) are proposed analogous to the square circulant matrices used in Laughlin et al. (1993) and Stewart et al. (2003b). As explained below (in Section 4.1), the proposed class of RCMs is related through a permutation to the rectangular matrices of square circulant blocks used in Chan, Nagy, and Plemmons (1992, 1993) and Davis (1994). The RCMs as presented here may also be stacked together to construct the class of g-circulant square matrices described in Davis (1994). Regarding the control problem at hand, the RCMs have a similar advantage to square circulant matrices as they can be transformed into diagonal blocks by multiplication with Fourier matrices (Davis, 1994; Hovd & Skogestad, 1994), with each diagonal element containing a SISO transfer function across the spatial frequencies. Using twodimensional (temporal and spatial) frequency analysis, this paper defines a two-dimensional closed-loop bandwidth and a robust stability margin for CD control systems. Model predictive control (MPC) has been very successful in the process industries (Qin & Badgwell, 1996). However, until recently (Backström et al., 2001), MPC has been rarely used in industrial CD control except a simplified optimization version (Chen & Wilhelm, 1986). The main reasons being: • Size of control problem: One single actuator array and one controlled property CD system could involve up to 300 identical actuators, 2000 measurements, and 1200 hard input constraints. • Perceived complexity: CD process outputs are affected by both temporal and spatial responses; a wide range of possible responses commonly occur on a working paper machine. • Robustness: CD process models are often severely illconditioned and are generally accompanied by model uncertainty. These issues have been considered for some time in the CD control research community. The large-scale problem may be solved by reducing the model size (Fan & Dumont, 2001; Haznedar & Arkun, 2002; Heath, 1996; Kristinsson & Dumont, 1996; Rigopoulos, 1999) or finding a fast quadratic programming (QP) algorithm (Bartlett, Biegler, Backström, & Gopal, 2002; Campbell, Rawlings, & Rao,
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
1998; VanAntwerp & Braatz, 2000) or using a linear programming (LP) solver (Dave, Doyle III, & Pekny, 1999). Although there is significant literature about robust MPC for handling model uncertainty, robust constrained MPC may not yet be implementable in real large-scale industrial processes due to the intensive on-line computational time required (Bemporad & Morari, 1999; Mayne, Rawlings, Rao, & Scokaert, 2000). Refer to (Duncan, 1995; Featherstone et al., 2000; Laughlin et al., 1993; Stewart et al. 2003a, b; VanAntwerp, Featherstone, & Braatz, 2001) and references therein for other robust controller structures related to CD processes. The new state of the art in industrial CD control is the model predictive controller described in Backström et al. (2001), Backström, Henderson, and Stewart (2002) and Bartlett et al. (2002), which uses a realistic rectangular process model identified from input–output data by commercial identification software (Gorinevsky & Gheorghe, 2003; Gorinevsky & Heaven, 2001). At the time of writing, this CD-MPC has been used in more than 10 real paper machines and shows advantages over traditional CD controllers which are restricted to have block-decentralized structures. The advantages are due to the coordination of multiple actuator arrays. However, there is currently no consistent method for evaluating the selection of the optimization weights in the model predictive controller to achieve the desired performance. Therefore, in this paper, we mainly focus on analysis of the stability and the performance of CD-MPC in the presence of inevitable model uncertainty. The contributions of this article are: (1) the introduction of the concept of the rectangular circulant matrix (RCM) and its properties; (2) the application of these properties to openloop and closed-loop transfer matrices of CD systems in order to decouple these large-scale MIMO transfer matrices into a family of SISO transfer functions across the spatial frequencies by using the two-dimensional frequency analysis method; (3) consistent criteria are given for evaluating the selection of optimization weights of the CD-MPC for closedloop CD systems.
1893
process, and h(z) is the Z-transform of the first order with deadtime temporal response of the process. For illustration purposes throughout this article, we will utilize an example model representing a typical system of a slice lip actuator array used to control the basis weight profile of a newsprint paper machine. Fig. 2a illustrates the nominal steady-state response of the basis weight profile to slice lip actuators. The spatial interaction matrix G0 in (2) is a band diagonal matrix as illustrated in Fig. 2b. This band diagonal and nearly spatially invariant structure of the spatial interaction matrix G0 in (2) is very important and is extensively used in the subsequent sections of this paper. We also assume that the true process response belongs to a set of possible response models, which is described by Gp (z) ∈ : ={G(z) + G (z) : (G (ei )) < G(z) ∞ }, (3) where ∈ (0, 1) is a positive scalar number which limits the perturbed transfer matrix Gp (z) to a neighborhood of the nominal model G(z) in (1), denotes the maximum singular value, z=ei and the symbol denotes the temporal frequency. 2.2. Industrial CD model predictive controller This article considers a recently developed industrial MPC controller for CD processes that has been gaining widespread application in paper mills. This industrial controller is developed by making use of the sparse structure of G0 in (2) (Backström et al., 2001) and choosing an efficient optimization algorithm (Bartlett et al., 2002). In this industrial MPC, the quadratic programming problem is proposed as follows: Hp (y(k ˆ + j ) − ysp (k + j ))T Q1 (y(k ˆ + j) min u(k) j =1
−ysp (k + j )) +
H u −1
[u(k + j )T Q2 u(k + j )
j =0
2. Problem statement
+u(k + j )T Q3 u(k + j )]
(4)
2.1. Process model The standard paper machine CD process model is used in this work and given by y(z) = G(z)u(z) + d(z), G(z) = G0 h(z),
h(z) =
(1) z−Td , 1 − az−1
(2)
where y(z) ∈ Cm , u(z) ∈ Cn , and d(z) ∈ Cm are the Z-transforms of the measurement profile, the actuator setpoint profile and the disturbances, respectively. The constant matrix G0 ∈ Rm×n describes the spatial response of the
subject to (1) and the following constraints on the actuator array:
u(k)b − Υ u(k − 1),
(5)
where yˆ ∈ Rm is the predicted output profile, ysp ∈ Rm are the setpoints for controlled sheet property, Hp and Hu are the prediction horizon and the control horizon, respectively, Q1 ∈ Rm×m and Q2 ∈ Rn×n are diagonal weighting matrices for the controlled sheet property y and the changes in the actuator array u(k)(=u(k) − u(k − 1)), respectively, and Q3 ∈ Rn×n is the weighting matrix which is either a
1894
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
diagonal matrix or given by Q3 = B T B, −1 1 1 −2 0 1 . .. B = .. . . .. . . . . .. .. . 0 ···
Gp(z)
(6) 0 1 −2 .. . .. . .. .
··· ··· .. .. . . .. .. . . .. .. . . .. . −2 .. . 1 ··· 0
··· .. . .. . .. .
0 .. . .. . .. .
, 0 1 −1
ysp
Kr
+
K(z)
d(z)
∆G (z)
u(z)
G(z)
+ +
+ +
y(z)
+
(7)
Fig. 3. The closed-loop system block diagram with the unconstrained MPC in (4).
where is a positive scalar, B ∈ Rn×n is known as the “bending moment matrix” (Bartlett et al., 2002; VanAntwerp & Braatz, 2000); Q1 , Q2 , and Q3 are positive semi-definite, , Υ, and b are constraint matrices and vector, respectively, derived from consideration of the physical limitations on the actuators and the details may be found in Backström et al. (2001, 2002). The design degrees of freedom available in the industrial model predictive controller (4) are the weighting matrices Q1 , Q2 , and Q3 , the prediction horizon Hp , and the control horizon Hu . These parameters should be designed such that:
ing their unconstrained behavior (for example loop shaping, pole placement, internal model control, etc.). Indeed many practical CD control systems, if well tuned, will not have active constraints. If constraints (5) are neglected, then the system is linear and we may write down a closed-form solution for these operations (Maciejowski, 2001). Fig. 3 illustrates the block diagram structure for such a system. The derivation of the constant prefilter matrix Kr ∈ Rm×m and the complex transfer matrix K(z) ∈ Cn×m from (1)–(2) and (4) is standard (Maciejowski, 2001) and the details may be found in Fan (2003). The linear closed-loop system in Fig. 3 is robustly stable for all plants Gp (z) in (3) if it is nominally stable and
···
1 −2 1
• the closed-loop control system given by (4)–(5) and (1)–(2) is stable for all Gp (z) ∈ in (3); • the estimated standard deviation of the controlled sheet property profile at steady state is kept as small as possible. 3. Closed-loop transfer functions Strictly speaking, the problem as described above in Section 2.2 would require analyzing the robust stability of a closed-loop control system containing a nonlinear optimization. While such a solution would be beneficial (Kothare, Balakrishnan, & Morari, 1996; Kouvaritakis, Rossiter, & Schuurmans, 2000), it is quite challenging for CD processes due to their large-scale characteristics (Bemporad & Morari, 1999; Mayne et al., 2000). In this section, we propose a technique based upon a linear approximation to the problem. We proceed as follows: (i) We first ignore constraints (5) such that the closed-loop system given by (1)–(2) and (4) is linear. (ii) We then compute equivalent transfer matrices of the CD-MPC in (4). (iii) We compute the closed-loop transfer matrices and analyze the performance and robustness of the system. (iv) Finally, we re-introduce constraints (5) for implementation. While it is generally recognized that the performance of a control system is affected by its actuator constraints, the practical tuning of controllers often concentrates on design-
Tud (z)G (z) ∞ < 1 ⇐ (Tud (ei )) <
1 , (G (ei ))
Tud (z) = K(z)[I − G(z)K(z)]−1 ,
∀,
(8) (9)
where G (z) is the unstructured uncertainty in (3) and result (8) follows from the small gain theorem (Skogestad & Postlethwaite, 1996), Tud (z) ∈ Cn×m is the nominal closedloop transfer matrix connecting the actuator setpoints u(z) to the output disturbance d(z). As CD control is largely a regulation problem with only occasional setpoint changes, the performance criterion is a function of its disturbance rejection properties. The nominal transfer matrix linking the measured profile y(z) with the output disturbances d(z) is obtained from Fig. 3 as Tyd (z) = [I − G(z)K(z)]−1 ,
(10)
where Tyd (z) ∈ Cm×m is also known as the sensitivity function (Skogestad & Postlethwaite, 1996). For stable plants like CD processes, properly choosing the prediction and the control horizon Hp and Hu in (4) can guarantee that Tud (z) and Tyd (z) are stable (García, Prett, & Morari, 1989; Mayne et al., 2000). The two transfer matrices Tud (z) in (9) and Tyd (z) in (10) will be used to analyze the behavior of the linear closed-loop system illustrated in Fig. 3. However, as Tud (z) ∈ Cn×m and Tyd (z) ∈ Cm×m have very large dimensions (n300, m2000), the analysis of their properties can be quite complicated. The next two sections introduce a framework and present a two-dimensional frequency domain analysis approach for solving this problem.
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
4. Rectangular circulant matrices In this section, we first propose the concept of a rectangular circulant matrix and then show that the large-scale closedloop transfer matrices Tyd (z) in (10) and Tud (z) in (9) are an approximate square circulant matrix (SCM) and an approximate rectangular circulant matrix (RCM), respectively. 4.1. Rectangular circulant matrices A circulant matrix is a particular square matrix whose definition and properties may be found in Davis (1994) and Gray (2002). Here, we propose the definition of a rectangular circulant matrix (RCM): Definition 1. A matrix N ∈ Cm×n is a rectangular circulant matrix if it satisfies the following conditions: (1) If m = na n, then column j = column j − 1 circularly down shifted na times, j = 2, 3, . . . , n; (2) If n = na m, then row k = row k − 1 circularly right shifted na times, k = 2, 3, . . . , m; for a positive integer na . The class of RCMs defined above is related through a permutation to the class of rectangular matrices R composed of na square circulant blocks C1 , . . . , Cna ∈ Cn×n described in Chan et al. (1992, 1993) and Davis (1994). That is, R = P N for m = na n (and R = N P for n = na m) where P ∈ Rmax(m,n)×max(m,n) is a permutation matrix (Fan, 2003). Also, the class of g-circulant matrices defined in Davis (1994) may be constructed from RCMs as presented above in Definition 1. That is, for m = na n, the square matrix Nc : =[N, . . . , N] constructed by stacking together na RCMs is an na -circulant matrix by Definition 5.1.1 in Davis (1994) (the n = na m case follows by a straightforward extension). The motivation for selecting the RCM structure as in Definition 1 for the CD control problem at hand is presented in Section 4.2. Rectangular circulant matrices (RCMs) are analogous to square circulant matrices (SCMs) with respect to several properties (see (Fan, 2003) for the proofs of these properties): (1) If N1 , N2 ∈ Cm×n are both RCMs, then N3 = N1 + N2 is an RCM. (2) If N1 ∈ Cm1 ×n1 and N2 ∈ Cm1 ×n2 are both RCMs, then N4 = N1T N2 ∈ Cn1 ×n2 is an RCM for n1 = na n2 or n2 = na n1 . (3) If N1 ∈ Cm×n is an RCM with m = n and N1 is invertible, then N5 = N1−1 is an RCM. = (4) The matrix N ∈ Cm×n is an RCM if and only if N H Fm N F n is a matrix of diagonal blocks: T = [N1 , . . . , Nna ] N na ] [N1 , . . . , N
if m = na n, if n = na m,
(11)
1895
where Fm and Fn are the Fourier matrices,2 H denotes k ∈ Cmin(m,n)×min(m,n) is a conjugate transpose, each N diagonal matrix for k = 1, . . . , na . = Fm N F H (5) If N ∈ Cm×n is an RCM, then N n has the diagonal block structure of (11), and ) j (N )=j (N n a (j + (k − 1)n, j )|2 |N k=1 = na |N (j, j + (k − 1)m)|2 k=1
j = 1, . . . , min(m, n),
if m = na n, if n = na m,
k = 1, . . . , na .
(12)
Note that the singular values j (N ) in (12) are not necessarily sorted in descending order as usually presented, but rather ordered according to the frequency content of the corresponding singular vector. (6) If N ∈ Cm×n is an RCM and nd ∈ Cq (q = max(m, n)) is a vector which represents the diagonal elements of 1 , . . . , N na in (11), N 0 N1 0 · · · 0 0 N2 · · · , nd : =diag(N), N= .. .. .. ... . . . na 0 0 ··· N (13) then, the elements of nd occur in conjugate pairs nd (j + 1) = nd (q − j + 1)H ,
j = 1, . . . , k,
(14)
where k = q/2 − 1 for even q and k = (q − 1)/2 for odd q. Note that nd (1) is the DC component, and nd (q/2) for even q is the frequency component at = 2 which is not conjugate with nd (1). Fig. 4 shows non-zero elements of a SCM Ns ∈ C320×320 , an RCM Nr ∈ C320×64 and the non-zero diagonal elements s and N r . of N 4.2. Frequency domain interpretation Given the spatial signals y ∈ Cm and u ∈ Cn , then the formulas y = Fm y
and u = Fn u,
(15)
where multiplication by Fourier matrices Fm and Fn is equivalent to performing the standard m-point and n-point discrete Fourier transforms, respectively. The usual interpretation is that this is a decomposition of the (spatial) signals into spatial harmonic functions with m (respectively n) equally spaced frequencies between = 0 and 2 (Oppenheim, Schafer, & Buck, 1999). 2 The (j, k)th element of the Fourier matrix of order n is F (j, k) = n
√
−(j −1)(k−1) / n with j, k = 1, . . . , n and = e2 i/n , i2 = −1.
1896
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
spatial frequencies are the spatial Nyquist frequency of the measurements and the actuators respectively defined by
0
50
50
100
100
150
150
200
200
250
250
300
300
y
N =
0 100
100
200
200
300 0
300 100 200 300
(a)
100
200
300
(b)
(c)
20
40
60
0
20 40 60
(d)
Fig. 4. (a) Non-zero elements of a SCM Ns ; (b) Non-zeros of the diagonal s ; (c) Non-zero elements of an RCM Nr ; (d) Non-zeros of elements of N r . the diagonal elements of N
We will consider the case m=na n in Definition 1 in some detail (the case n = na m is a straightforward extension). If u, with =Fm N F H y =N u, then y = N y and u in (15) and N n. From (11) we find that for an input signal u of a single spatial frequency u (j ), 1 for i = j, u(i) = (16) 0 otherwise,
u (j ) = 2
j −1 , n
j = 1, 2, . . . , n,
(17)
is composed of na spatial frequency components, 2 j − 1 y (j + (k − 1)n) = + (k − 1) , na n for k = 1, . . . , na .
uN =
1 . 2Xu
f(j ) = 0
for jb j m − jb .
B =
1 y (jb ). 2 X y
where X denotes the physical spacing of the elements, and typical units for the spatial frequency are cycles/m. For CD control systems in general the physical spacing of actuators and measurements Xu > Xy due to the fact that there are typically several times more measurements than actuators distributed across the same spatial domain. Two important
(23)
Theorem 1. If N ∈ Cm×n is an RCM, y = N u, and the spatial bandwidth in (23) is limited to B uN (the spatial Nyquist frequency of input u in (21)), then (i, j ) N u(j ), i = j + k, y (i) = (24) 0 otherwise j = 1, . . . , f loor(n/2) + 1, j = f loor(n/2) + 2, . . . , n.
Proof. See Fan (2003).
Each of the discrete spatial frequencies u or y may be converted into engineering units with the usual formula 1 0 , 2 X , (20) = 1 2 X (2 − ), < < 2 ,
(22)
The discrete indices in (22) correspond to a physical spatial bandwidth B from (20),
(18)
(19)
(21)
Although the general case in (16)–(19) will show the presence of resonant frequencies, practical CD control systems are typically dominated by the response at a single output spatial frequency due to the process response being “bandlimited” in the majority of cases. Denote the j th column of N ∈ Cm×n as nj ∈ Cm , and define the Fourier transform of nj as nj = Fm nj . Note that due to the circular nature of the columns of an RCM, then | nj | = | nk | for any j, k=1, . . . , n and is denoted as f=| nj | for any j =1, . . . , n, where f ∈ Rm contains the magnitude of the spatial frequency components of the columns of N. A bandlimited fis one for which there exists an index jb < f loor(m/2) such that3
where 0, k= m − n,
then the output y (j + (k − 1)n) + (k − 1)n, j ) N(j u(j ), k = 1, . . . , na , = 0 otherwise
1 , 2Xy
(25)
This result indicates that when the column vectors of N are bandlimited to the spatial Nyquist frequency of inputs uN , the matrix operator N can be replaced by a family of scalar equations, one for each discrete spatial frequency u (j )with j = 1, 2, . . . , n in (17). Corollary 1. If N ∈ Cm×n is an RCM and B uN in (21) and (23), then (i, j )|, j (N ) = |N
(26)
where i = j for j = 1, . . . , f loor(n/2) + 1, or i = j + m − n for j = f loor(n/2) + 2, . . . , n. Proof. See Fan (2003).
3 f loor(x) is a Matlab command for rounding the elements of x to the nearest integer towards minus infinity.
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
1897
Based on the RCM’s properties and Theorem 1, the following theorem may be presented:
for downsampling from the high resolution m1 to the true measurement resolution m. Then, the model matrix G0 in (2) is given by
Theorem 2. For the closed-loop system in Fig. 3, if
G0 = G2 · G1 ∈ Rm×n .
m×n
(1) G0 ∈ R in (2) is an RCM by Definition 1, (2) the columns of G0 are bandlimited with B uN in (21) and (23), and (3) the weighting matrices Q1 ∈ Rm×m , Q2 , Q3 ∈ Rn×n in (4) are RCMs, then, the open-loop plant model G(z) ∈ Cm×n in (2), the closed-loop transfer matrices Tud (z) ∈ Cn×m in (9) and Tyd (z) ∈ Cm×m in (10) are also RCMs. Proof. See Fan (2003).
4.3. Approximate rectangular circulant matrices Theorem 2 shows that the RCM structure would be preserved for each of the transfer matrices important to the problem at hand. For practical CD control systems, there are two violations to the assumptions in Theorem 2: (1) The static matrices G0 in (2) and Q3 in (6) are not exact RCMs due to edge effects (see Fan, 2003; Mijanovic, 2004). (2) The dimensions of the m × n matrix G0 in (2) do not satisfy m = na n for integer na . The column vectors of the spatial interaction matrix G0 in (2) are circularly shifted except the column vectors near the edges (see Figs. 2a and b). Comparing Fig. 2b with Fig. 4c (which illustrates a true RCM), these two figures are exactly same except for the upper right corner and the lower left corner. If Q3 = B T B with the bending moment matrix B in (7), it is not an exact RCM due to several edge (near its 4 corners) elements in the first and last rows and columns. The difference between an RCM and its approximation arising from the above spatial edge problems (of G0 and Q3 ) has been investigated in Fan (2003). The spatial edge effects on the RCM approximations of closed-loop transfer matrices Tyd (z) and Tud (z) are minor. As for the second violation, a technique analogous to that used to change the sampling frequency of a continuous signal in multirate signal processing (Vaidyanathan, 1992) may be used to obtain a useful result by upsampling and downsampling in turn. First, create a new spatial interaction matrix G1 ∈ Rm1 ×n where m1 = m · n by upsampling the continuous spatial response (i.e., the columns of G0 ) at the corresponding high resolution. Then, create a second matrix
1 0 G2 = ... ∈R
0
m×m1
··· 0 0 ··· 0 1 . . .. . .. .. 0 ··· 0 0
0 0 .. .
n
··· 0 ··· 0 ··· 0 ··· 0 . . . .. . .. . . .. 0 ··· 0 ··· 1
0 0 .. .
n
··· 0 · · · 0 . .. . .. 0 ··· 0
0 0 .. .
n
(27)
(28)
The downsampling matrix G2 in (27) is an RCM by Definition 1, and we will approximate the upsampled model matrix G1 ∈ Rm1 ×n with an RCM through the use of periodic boundary conditions. However, G0 = G2 · G1 itself is not an RCM due to the fact that the dimensions of G2 ∈ Rm×m1 and G1 ∈ Rm1 ×n do not satisfy Property 2 of the RCM. 2 = Fm G2 FmH and From the RCM’s property 4, both G 1 1 = Fm1 G1 FnH are each a matrix of diagonal blocks as G illustrated in (11), and further, the column vectors of G1 satisfy the condition B uN of Theorem 1. Therefore, 0 = G 1 ∈ Cm×n has the following characteristic: 2 G G 0 (i, j ) = 0 for i = j + k, (29) G = 0 otherwise, where k is given in (25), but in this case m is not restricted to be an integer multiple of n. Comparing (29) with (24), G0 ∈ Cm×n with a non-integer na = m/n will have a properly analogous property to that of the RCM N in Theorem 1. Furthermore, the closed-loop transfer matrix Tud (z) will have a similar property as that of the RCM N in Theorem 1 and the sensitivity function Tyd (z) will be a SCM (Fan, 2003). 5. Two-dimensional frequency domain analysis In this section, we will first examine the closed-loop system behavior of the control system in Fig. 3, and then discuss how the tuning parameters affect the closed-loop system’s robustness and performance. 5.1. Transfer functions in the two-dimensional frequency domain From the previous section, we know that the open-loop plant model G(z) in (1), the closed-loop transfer functions Tyd (z) and Tud (z) are approximate RCMs. Their corresponding representations in the two-dimensional frequency domain can then be obtained by y
y
{ g (1 , z), . . . , g (n , z)} = Diag(Fm G(z)FnH ),
(30)
y y tyd (m , z)} = diag(Fm Tyd (z)FmH ), { tyd (1 , z), . . . ,
(31)
{ tud (u1 , z), . . . , tud (un , z)} = Diag(Fn Tud (z)FmH ),
(32)
y
where j and uj are spatial frequencies with engineering units in (20), ‘diag(T )’ denotes the operation of extracting the diagonal elements of T, ‘Diag(T )’ denotes the operation of obtaining the following elements of a rectangular matrix4 4 It is a straightforward extension for the case T ∈ Cm×n , such as G(z) in (2).
1898
J. Fan et al. / Automatica 40 (2004) 1891 – 1903
2
0 0 10 -2 10
-4
10
-6
10
-8
10
0
1
2
3
4
5
iω
40
u j
0.5 0 0 10 -2 10
temporal frequency spatial frequency [cycles/second] [cycles/metre]
(a)
1
ud
y
0.02
60
1.5
|t (ν ,e )|
iω
|tyd(νj ,e )|
0.04
y j
iω
|g(ν ,e )|
0.06
-4
10
-6
10
-8
10
0
1
2
3
4
5
20
-2
10
-6
10
-8
10
0
1
2
3
4
temporal frequency spatial frequency [cycles/metre] [cycles/second]
temporal frequency spatial frequency [cycles/second] [cycles/metre]
(b)
-4
10
(c)
Fig. 5. The two-dimensional frequency representations of (a) the nominal plant model G(z), (b) the sensitivity function Tyd (z), and (c) the closed-loop transfer matrix Tud (z).
T ∈ Cn×m : Diag(T )={T (1, 1), . . . , T (k, k), T (k + 1, k + 1 + m − n), . . . , T (n, m)},
(33)
where k = n/2 for even n or k = (n + 1)/2 for odd n. y Fig. 5a shows the nominal model’s representation g (j , z) in (30) in the two-dimensional frequency domain. Fig. 5a y illustrates that the process model g (j , z) in (30) has large gain for low spatial and temporal frequencies, and small gain for high spatial and temporal frequencies. This highfrequency gain roll-off is a familiar feature of physical processes and limits the controllability of disturbances d(z) in (1) to the low spatial and temporal frequencies (Duncan & Bryant, 1997; Featherstone et al., 2000; Stewart et al., 2003b). The ill-conditioning of G0 in (2) can be interpreted using (26) in Corollary 1 as gain roll-off for high spatial frequencies. Figs. 5b and c show the two-dimensional frequency repy resentations tyd (j , z) in (31) and tud (uj , z) in (32) of the closed-loop transfer matrices Tyd (z) and Tud (z) in (10) and (9), respectively. Here again the nominal model G(z) in (1) of Section 2.1 is used, and horizons Hp = 8 and Hu = 1, and weights Q1 = 4I, Q2 = 10−4 I , and Q3 = 0.002I in (4). The two-dimensional frequency sensitivity function y tyd (j , z) in (31) may be used for calculating the control system’s bandwidth. The bandwidth in the two-dimensional frequency domain can be similarly defined as traditionally done in the temporal frequency domain (Skogestad & Postlethwaite, 1996): y
Definition 2. The closed-loop bandwidth Bw (j , ) y of CD processes is the contour in (j , ) such that y y y | tyd (j , ei )| < √1 for all (j , ) enclosed by Bw (j , ). 2
The standard robust stability condition in (8) may be analogously transferred into the two-dimensional frequency domain. If Tud (z) in (9) is an RCM, then, based on (26) in
Corollary 1, the maximum singular value of Tud (z) is equal to the maximum value of | tud (uj , z)| in (32). Therefore, the robust stability condition (8) for the closed-loop RCM system may be rewritten as | tud (uj , ei )| < 1/(G (ei )),
∀uj , .
(34)
This allows us to define a robust stability margin for unstructured additive model uncertainty in (3) as
: =1/ max sup | tud (uj , ei )|, j
(35)
for uj ∈ {u1 , u2 , . . . , un } and ∈ [0, N ], where N is the temporal Nyquist frequency (N = ). This robust stability margin can be explained as follows: given the nominal process model G(z) in (1) and the controller K(z) in Fig. 3, for any stable unstructured additive uncertainty G (z) in (3), if (G (ei )) < for all , then the closed-loop system described in Fig. 3 is stable. 5.2. Role of MPC tuning parameters Typically, the design of the horizons Hp and Hu in (4) is based on the dynamics of the process (Alamir & Bornard, 1995; Mayne et al., 2000; Soeterboek, 1992) and is not explored further here. In practice, Hu > 1 is rarely used as it leads to heavier on-line computational burden for solving the QP problem in (4) (Bartlett et al., 2002). The number of degrees of freedom in the selection of tuning parameters in (4) can be reduced by noting that without loss of generality, Q1 in (4) may be fixed and the closed-loop performance modified via Q2 and Q3 (Soeterboek, 1992). Typically in MPC, the weighting matrix Q2 is considered to affect the dynamics of the closed-loop system while Q3 affects the slower behavior and steady-state performance. Here we explore the justification of this assumption. Let us revisit the sensitivity function Tyd (z) in (10) in order to evaluate the weights’ influence on the performance. The controller K(z) in Fig. 3 and (10) can be obtained after
J. Fan et al. / Automatica 40 (2004) 1891 – 1903 2 spatial frequency [cycles/metre]
4.5 Q =0.0001I 2
Q2=0.001I Q2=0.01I
iω
|tyd(0,e )|
1.5
1
0.5
0
0
0.01
(a)
0.02
2
Q =0.01I
3
2
2 1.5 1
|tyd(νyj ,eiω)|
