- Email: [email protected]
Long-range predictive control with a terminal matching condition
Pergamon
Chnntcof E#a&eering Scimce, Vol.
49. No. 9, p,,. 1287%13W, ,994 Copyright 8 1994 Elsevicr Science Ltd Printed in Gsal Britain. All rights rcsorvcd OOIX-2509,Y.t 56.00 + 0.00
LONG-RANGE
PREDICTIVE CONTROL WITH A TERMINAL MATCHING CONDITION KUN-YU
Department
KWOK
and SIRISH L. SHAH’
of Chemical Engineering, University of Alberta, Edmonton,
Alberta, Canada T6G 2G6
(First receiued 27 December 1991; accepted in revisedfirm 6 September 1993) Abstract-The long-range predictive control strategy is modified by augmenting the generalized predictive control law with a terminal matching condition. This condition is formulated as a weighting term on the output steady-state error. The final control law only minimizes the square of prediction errors over a small future prediction horizon and at steady state. It is applicable to either transfer function or convolution modeling. For transfer function models, the modified approach reduces the computational load of an adaptive controller by mimicking a large prediction horizon with a smaller one plus a weighting at the steady state. For convolution models, only a few initial step response coefficients and the steady-state gain are required to formulate the new long-range predictive controller. The nonadaptive closed-loop analysis reveals that this new steady-state error weighting term (i) reduces the additional order otherwise required due to inclusion of an integrator for offset elimination, and (ii) provides the stabilizing effect of a large predictive control horizon. A closed-loop adaptive control system based on this modified control law has been successfully applied to the control of mean arterial blood pressure. 1. INTRODUCTION
for minimum variance (bistriim and Wittenmark, 1973) and generalized minimum variance control (Clarke and Gawthrop, 1975) were aimed at providing tight control of processes by minimizing the variance of the process output predictions. A key factor in the successful use of these methods is correct prior selections of the model order, n, and time delay, k, such that accurate k-step-ahead predictions are made. Poor selection of k or a varying time delay affects the control performance or even leads to instability. Improper selections of model order may result in similar problems although overspecified model orders may not be as sensitive as the uncertainties in time delay. Recent developments in long-range predictive control (LRPC) which attempt to minimize the sum of prediction errors over a user-specified time horizon provide an alternative solution to the abovementioned drawbacks (Richalet et al., 1978; Cutler and Ramaker, 1980; Clarke and Mohtadi, 1989). Since the strength of the LRPC algorithm lies in the longrange predictions of process outputs and control actions, these two predictions which have a direct influence on stability and dynamic response in any longrange predictive controller (Garcia and Morari, 1982) naturally become the tuning parameters. The selection of the output prediction horizon is less influential than the choice of the time delay or model order. However, the long-range predictive nature necessarily makes the computation load of an adaptive controller very heavy for each new control action to be calculated. On the other hand, a relatively short output The
early
methods
prediction horizon, of which the minimum variance is the extreme case, may result in vigorous control actions or even unstable control. Moreover, the longrange control algorithms such as MAC (Richalet et
‘Author
to whom correspondence
should be addressed.
al., 1978), DMC (Cutler and Ramaker, 1980), and MOCCA (Sripada and Fisher, 1985). which employ a convolution model as the internal model (Garcia et al., 1989) to encompass most of the process dynamics, require as many as 50 step or impulse response coefficients, a size not economically sound for adaptive control. This paper proposes a modified approach to the LRPC scheme which introduces an output terminal condition in the minimization of a multistep prediction cost function. The terminal matching condition is applied in the form of weighting on the square of the error between the predicted process steady-state value and the desired steady-state target or setpoint. This weighting combined with a relatively short output prediction horizon, in the context of adaptive and nonadaptive control, is able to resolve the trade-off between economic computation and long-range robustness. The concept of steady-state error weighting is depicted in Fig. 1. An ordinary LRPC algorithm considers process output predictions over a large horizon from nl to nlonp.The modified approach proposed here would consider a relatively shorter horizon from nl to nsho,,plus the term at steady state so as to ensure the terminal matching of process output to the setpoint. Different forms of terminal conditions have appeared in the literature. However, most of the applications have occurred in the context of an infinite LQ controller with terminal state constraints for the improvement of stability (Chesna and Ydstie, 1988; Mayne and Michalska, 1991; Rawlings and Muske, 1992), and in a self-tuning controller with a cautious control element during a poor identification period (Papadoulis et al., 1987), but the concept of a terminal matching condition formulated as a weighted square of the steady-state error and the relative weight of this term as an on-line parameter has not been investigated in the literature. This concept proposed in this
KUN-YU KWOKand
1288
SIRISH
L. SHAH
(4) There are recommendations of two performance tuning strategies for LRPC with steadystate error weighting: one for fast response, the other for robust response. (5) The successful experimental evaluation of this (adaptive) control strategy in the control of mean arterial blood pressure.
Fig. 1. Long-range predictive control weighting.
with
steady-state error
paper minimizes a finite-time objective function augmented by a weighted infinite-time element (for the terminal condition), and is therefore functionally different from infinite-horizon LQ strategies, i.e. the proposed controller is simultaneously able to accommodate finite-time and infinite-time objectives. In this study, the main contributions, in the context of adaptive and nonadaptive control of processes, are:
(1) The incorporation
of a weighted steady-state term as a tuning parameter into the performance criteria of LRPC. Philosophically speaking, the idea of a terminal matching condition is intuitively appealing and easy to relate to by plant engineers. The concept of terminal matching can be readily explained when one considers the analogy of a car trip from point A to destination B. Whilst the performance of the vehicle over a finite and continuously moving window of say 200 m is important, equally important is the final destination B. The steady-state weighting gives us the capability of incorporating this condition in the design philosophy. (2) The incorporation of such weighting in the optimization criterion obviously requires steadystate output predictions. In order to reduce the heavy computational load as used in conventional steady-state prediction techniques, technical results presented in Section 2 provide a method to make steady-state predictions with fewer iterations. The method allows steadystate output predictions to be computed directly from the ARIMAX model without iterating to a large output prediction horizon. The same method can also apply to the convolution model in that a complete set of step response coefficients is not required. (3) The extra closed-loop model order incurred as a result of the integral model is eliminated by using the steady-state error weighting.
This paper first presents the procedure for on-line model-based finite-horizon and steady-state predictions of the process output responses based on ARIMAX or convolution models. Then the control law with the steady-state error weighting and its closed-loop analysis are discussed in the context of nonadaptive control, followed by the performance tuning strategy. The last section of this paper documents the results of the application of this modified LRPC strategy in the control of mean arterial blood pressure. 2.
PROCESS MODELING
Most model predictive control techniques are based on the optimization of a quadratic objective function involving the error between the setpoint and the predicted outputs. The distinguishing feature among different techniques is the model type and structure for making the predictions. The most commonly used model types are discrete-time transfer functions describing the input-output relationship, or discrete-time convolution models in the form of step or impulse response coefficients. The use of both model types for finite-horizon plus steady-state predictions is considered in the following sections. 2.1. Transfer junction model The process model is assumed to be an autoregressive integrated moving-average model with an exogenous input (ARIMAX):
A(q-‘)y(t) where
= &q-‘)u(t
- 1) + T(q+$
(1)
A=l-q-’ B^= Bq-‘+’ 4 -’ = backward shift operator.
Having the integrator in the noise model also represents a class of nonstationary dynamics such as random steps at random times, i.e. random-walk or Brownian motion type disturbances. Note that this integral mode serves as an internal model of the (steptype) disturbances and setpoint changes (Francis and Wonham, 1976). To derive ajth step-ahead predictor from eq. (i), the following Diophantine identity is considered: 7”= EjA^A + q-jFj.
(2)
Multiplying eq. (1) by AEjq’ gives E,AAy(t
+ j) = E,l?Au(t + j - 1) + EJTC; (t + j). (3)
Long-range predictivecontrol with a terminalmatchingcondition
1289
Substituting eq. (2) into eq. (3) results in the following jth step-ahead model:
and substituted into eq. (5), the result is
Ty(t + j) = E,BAu(t + j - 1) + F,y(t) -+ E,T<(t + j).
TAA
(4) Future control actions are separated from the known ones by substituting the following identity into eq. (4): (5)
Eji% = G,T + Hjq-j
where the elements of Gj now correspond to the first 1 open-loop step response coefficients. The optimal jth step-ahead predictor for j 2 1 becomes j -
(T - q-‘Fj)6
B q-‘Fja ,-,=G,++ AA
(6)
asj-bco.
qj($-G,)=++-%&
(7)
TAA l? -=Gj+q-j@+g) AA
(11)
For sufficiently large j, both Fj and Hi converge to F, and H,, so eq. (12) yields eq. (7). a (determination of F’, and H,): Recall that polynomials Fj+ 1 and Hj+ 1 can be defined as Remark
T=
Ej+,A^A + q-j-‘Fj+,
Ej+~~=Gj+~T+q-j-‘Hj+~.
04)
ejA^A = Fj - q-lFj,I ejB = gjT + q-‘Hj+
(15) 1 - H,
Proof: Rearranging eq. (2) gives = Ej + 9.
(8)
The denominator polynomial A^A has one root on the unit circle and nn roots inside the unit circle (na is the order of A^).According to the Taylor seriesexpansions of functions around the point on the unit circle for j approaching infinity, the tailing coefficients in the quotient polynomial Ej converge to the constant value T (l)//l^( 1) with order approaching infinity. Hence, the remainder q-jFj does not vanish but converges to a finite-length polynomial q-“F,, where Fj --+ F,
as+
CO.
E, =
?ZA
- H,
(18)
where T!I) es=A^o g.=-.
&l) A(l)
Therefore, F, is found from eq. (17) as F, = e,A^. Equation (18) becomes H,A = g,T - e,b.
When Ej in eq. (8) is rearranged as T - q-‘F,
(17)
ejg = e,6 = g,T + q-‘H,
z
(16)
where F and H are of degree na and max(nb, nt) - 1, respectively, and ej and g, are the last coefficients in Ej+ , and Gj+ 1. For sufficiently large j, F, H, and the tailing coefficients in E and G converge so that ejiA=e,AA=F,-qq’F,
T
(13)
Subtracting eqs (2) and (5) from eqs (13) and (14), respectively, gives
Lemma 1: Given two polynomials T and A with their roots all inside the unit circle. For sufJiciently large j, the
polynomials Fj in eq. (2) and Hj in eq. (5) conuerge to two Jinite constantpolynomials (de&d QSF, and H,), and the combination ofthe two identities in eqs (2)and (5) becomes
(IO)
Equation (10) indicates that the second Diophantine identity in eq. (5) also has one root on the unit circle and the rest inside the unit circle. Therefore, the same conclusion applies to eq. (5), i.e. H, converges to a finitelength constant polynomial as j approaches infmity. The last part of the lemma is verified by rearranging eq. (10) as follows:
$(f+jlt)=GjAu(t+j-l)+HjAuf(t-l)+Fjy,-(t) where uJ and y, are values filtered by l/T. Two methods of calculating the optimal steady-state prediction are avaiiable. The first method is simply by iterating eq. (6) for large values of j until $(t + jl t) converges. Each iteration requires the solutions of solving the Diophantine identities 2 and 5 with order j + degree (Fj, H,). As the prediction horizon jincreases, the computational load for the solutions also increases at every sampling interval. This tedious method is of course not recommended because the heavy computational load would defeat the purpose of our modified control strategy. Instead, a new and simpler method is described below to calculate the optimal steady-state output prediction at every sampling time. The following technical result in the form of a lemma is required for this simple prediction method.
=Gj+!+_
(19)
H, can be computed by comparing the coefficients in eq. (19). Thus, (9)
H, = h,,,, + h,, ,q-’
+
* - . + hS,nhqpmk
1290
KUN-YU KWOKand SIRISHL.
where nh = max(nb - 1, nt h‘.i = esjig,
1)
bj - g* ,=gl
ti*
q
Now applying the result in lemma 1 to the jth step ARIMAX model in eq. (6) yields (assuming unity time delay for simplicity)
= lim [G,Au(t tj j-02
- 1) + HjAu,(t
-
j-m
1)
-
y(t) = i giq-‘Au(t i- 1
lim [g,_,Atr(t) g,_,,.Au(t
+ gj_,Au(t
+ 1) +
lim 9(t +jl
a, = i=$+ 1 giq-‘Au(t
- 2) + goAu(t + j - l)] - k) + F,y,(t).
(20)
+ 1 - k)
i=
= m I t) = lim i-g_,Au(t)
+ g,._*Au(t + 1) +
i-m
+ IIU-
.. .
1)] + H,AU,(t - k)
+ F,y/(r) = [g%Au(t) + gSAu(r + 1) +
.. .
+ g,Au(t + nu - l)] + H,Aur(t - k) + FSY,0) = g,[Au(t) + Au(t + 1) +
-. .
+ Au@ + nu - 1)] + H,Au,(t - k) + Fa,-(t). (21) As a result, the optimal steady-state predictor based on an ARIMAX model is given as $(si t) = gs f Au(t +j i=i
1) + H,Au/ (t - 1) + FSyr(t). (22)
The first term on the right-hand side of eq. (22) is the current and future forced response, whereas the remaining terms corresponding to the impact of past
(truncation error).
T is usually chosen so that the truncation error is less than 1% of the step-input steady-state gain. Depending on the sampling rate, the number of coefficients may be as many as 50 to encompass 99% of the step response dynamics. The process model becomes
j(t) = f: g;q-‘Au(t
t)
+ gj-..Au(t
+ 1 - k) + ae + a,
T is the model horizon
For practical implementation, the projected control increments after a control horizon of nu are assumed to be zero (Cutler and Ramaker, 1980). Therefore, as gi converges to gs for increasing i, eq. (20) reduces to
j-m
(23)
where
. -
+ nu - 1) + * - ’
+ g,Au(t tj + H,Auf(t
+ 1 - k) + ass
In practice, the dynamic model is commonly assumed by a finite number of coefficients, i.e.
1)
FsY_r@)
j-.n
-t-
y(t) = 5 giq-‘Au(t 1=1
hi = Agi.
= lim [GjAu(t + j - l)] + H,Au,(t
=
2.2. Conuolution model The convolution model structure is usually represented by a step response model:
where the step response coefficients, gi, tend to g. when i tends to infinity and a,, is a reference value. It is readily turned into an impulse response model by differencing the step response such that
+ ~jY.f(t)l
+
SHAH
inputs and outputs on the prediction of the infinite time value of the process output.
+ 1 - k) + a_.
1
(24)
Its advantage over a parametric model is the capability to model unusual dynamic behavior which cannot be well represented by a reduced-order model. Moreover, many chemical processes are inherently infiniteorder and therefore a simple first- or second-order transfer function may introduce structural mismatch that is unfavorable to predictive control (Shook et al., 1992). However, on-line identification of such a large number of parameters is a difficult task. The process model in eq. (24) can be arranged into different predictive representations such as recursive form or state-space formulation (Garcia and Morari, 1982; Lim, 1988; Li et al., 1989). However, all of them assume a linear superposition of previous inputs to the total output such that (for unity time delay, i.e. k= I) jJ(t + jlt) = 5
giAu(t tj
- i)
i=l T +
,=F+, g&G
+j
- i) + ass,
(25)
The output comprises three terms: future control effects, past control effects, and a reference plus residual term. Equation (25) shows that even with steady-state weighting a complete step response sequence of up to
Long-range
T intervals is required regardless of the choice of output prediction, j. However, an approximation alternative, borrowed from Auslander et al. (1978), can be used so that a finite number of initial step response coefficients plus the steady-state gain are sufficient for long-range predictions. This approximation recognizes that most process responses are dominated by a single time constant at least in the low-frequency systems, the lowregion; i.e., for overdamped frequency dynamics can be well approximated by a first-order model. Consider that a unit step response from a stable second- or higher-order process is often characterized by an S-shaped curve, the low-frequency portion can be approximated by an exponential function from the nth point to the steady state. In terms of impulse response, the approximation is written as G,=
hlq-’
+ h2q-’
+ /&q-n + h,+Iq-n+’ =h,q-‘+h2q-2+
3. CONTROL
ALGORITHM
3.1. Derivation The design of a long-range, model-predictive controller is based on the predicted behavior of the process over the prediction horizon. The control objective is to minimize the errors between the trajectory of predicted values and setpoint over a horizon. The following quadratic function defines the performance index to be minimized to achieve the control objective:
r,(j) CP(t + j) - w(t +iW
.
+
...
to be A^= (1 - pq- ‘) such that the ARIMAX modeling conforms to an output error prediction type which is implicitly assumed in DMC type control algorithms that use the convolution model type.
JW, G n4 4 7) =
. . + hn_-lq-“+l
+
1291
predictive control with a terminal matching condition
+ F
hnq-”
+h,_lq-m+l+T
d(j)[Au(t
+j-
l)]’
ci -
1) - WI*
j=l
1 -pq(26) where p. being the rate of exponential decay, is a discrete characterization of the dominant time constant. If the steady-state gain of the process is known, for example, from steady-state design equations (Chesna and Ydstie, 1986), p is determined so that the model gain in eq. (26) is equal to the process gain:
h, LL - J$‘r;
hi’
(27)
The model in eq. (26) can be rearranged into a transfer function from G =b;q-LJ;b;q-2+ ... P 1 - pq-1
+b:q-” (28)
where b; = hI fori=2-n.
bi=hi-phi-1
A stable overdamped model is guaranteed as long as i= 1 It has been shown that the modeling error that stems from the exponential decay approximation can be made arbitrarily small by increasing the sampling period for a fixed n coefficients or vice versa (Takahashi et al., 1975). By transforming a “ limited” convolution model (a few step response coefficients and gain) into a transfer function form, the ARTMAX modeling procedure in Section 2.1 can be applied to make future predictions_ In doing so, the user-defined polynomial T is selected
5 r(j)CPM~ j=l
(29) where l
. l
G,(I) = 4s p=l-
+
l l l l
nl is the minimum output prediction horizon n2 is the maximum output prediction horizon nu is the control horizon such that Au@ + j) = 0, Vj >, nu y,(j) is an output weighting sequence d(j) is a control weighting sequence v(j) is a steady-state error weighting sequence s denotes a value at steady state.
The first two terms on the right-hand side form the standard generalized predictive control (GPC) objective. The last term corresponds to the additional terms penalizing the squares of steady-state predicted error. The summation from 1 to nu is required because, at each control interval j, where j = 1-nu, the steady-state prediction involves the sum of the first j consecutive control actions. It has been suggested that, for robust control, n2 should encompass the whole rise time of a process (Clarke et al., 1987), which often requires as many as 50 terms. Now with steadystate error weighting in place, n2 can be reduced to cover only the high-frequency dynamics. The process predicted values are required in eq. (29), from j = nl-n2 for j(t + j) and j = 1-nu for 9(slt + j - 1). To derive a control law from eq. (29), the predictive equations for the ARIMAX model in eqs (6) and (22) are used. The ARIMAX model can be expressed in the following compact vector/matrix equation (assuming unity delay): Y=GU+f
(30)
KUN-YU
1292
KWOK and SIRISH L. SHAH
GPC control law. Therefore, eq. (33) is referred to in short as “GPC with y-weighting”. For simplicity, the output weighting is assumed to be unity; thus, r, reduces to an identity matrix in the sequel.
where Y = [$(i + nl) P(t + nl + 1). -.
iL-1 9n1
..*
g, 0
0
*”
81
0 ..
Bo
-.
G:
9(t + n2)]’
.‘.
0
.-*
0
-.
90
-.-
CL2-
-
_..
clnl--I
9n2-2
II”
U = [Au(r) Au@ + 1) _. . Au@ + nu - l)lr f=[f(t+nl)f(t+nl+l) f(t +i) = H,Au,(t
...
f(t-tn2)]’
- 1) + Fj_Y,(t)
and Y, = GJJ + f,
(31)
where Y, = [J+Ic) j(slt
. Jvslt + nu - 1)lT
+ 1) .
rg. 0
3.2 Control law in linear form Although the control law in eq. (33) calculates the future control action from r = 1 to nu, only the current control action, i.e. Au(t), is implemented. At the next sample, the entire long-range minimization is repeated once again to calculate the next control action. This receding-horizon policy produces better control performance as the control law makes use of updated plant information. It has been shown that, in the absence of disturbance or model change, the control sequences in successive times are the same (Clarke and Mohtadi, 1989). Because of the receding-horizon policy, only the first element in vector U of eq. (33) is required. To enable closed-loop analysis of this GPC control law, the receding-horizon control law is expressed in a linear, difference equation form (Lambert, 1987; Mclntosh et al., 1991). A general linear structure of a controller with a setpoint and measurement variable is represented by
01
‘..
RAu(c) = VW(r) - Q(r) where R,V, and S are polynomials shift operator. Now define a vector
f&u x 1) = [l 1 . * . 13’CH.h, it - 1)+ F,Y, (dl. The equivalent of eq. (29) using vector/matrix notation is J = [Y - W]VJY
(32)
[GrG
W = [w(t
+ nl)
. - * y(m)]
w(r + nl + 1)
W,(nuxl)=
[1 1
..
. . .
w(r + n2)]r
+ A + Gfl-G,]
+ Gfr(W,
-
fJ1
Au(t) =
l]‘w(s).
- ‘[GTl-,(W
1~~~
w(t)
After substituting eqs. (30) and (31) into eq. (32), the minimization of the cost function yields the following control equation: U = [GTl-,G
(35)
.*a
+ (34)
+A2
+ A + GTl-G,] - ‘G:l-.
Then the current control action can be written as follows (assuming all future setpoints are the same as the present one),
A = diag [A(l) K(2) . . . A(m)]
r = diag [y(l) y(2)
a.21
...
where #Iicorresponds to the elements of the first row vector in the matrix
where r,. = diag Cry(nl) y,@l + 1) . . . y&2)3
in the backward
which denotes the first row of the matrix [GrG A + G:T&]-‘Gr, and a variable a,=Al+&+l
- W] + LJTAU
Wm-CY, - w.1
+ CY, -
Cm,, a,l+i
(34)
- f)
’ . . an2 a,]
- F.1 Y/N -F,i+,y#
44 x
an1 +1
- H., Au,+ - 1) -Kr+lAu,-(t-
1)
: - F,ay&)
w(t)
- H,,Au,(t
- 1)
- H,Au,(t i
- 1) I
- F,y, I (El
I w(t)
(33)
where GTT,G and GTTG. are both matrices of dimension nu x nu. The dynamic matrix, G, contains all the step response coefficients arranged in a lower triangular structure. The term G,TrG, is always of fun rank. It guarantees a nonsingular inverse even if the matrix G’T,,G is ill-conditioned due to large time delays or too short an output prediction horizon. When r is set to zero, eq. (33) reduces to the basic
(37) and, in compact Au(t) = w(r) 2
*=“I
- Aq(t
form, ai - yr(r)
- 1) 2
g F,a, I=“, ff{Kj
i=n, +
w(r)a. . , L-
v,(r)F.a. _, , _ - - Au&,_ - 1)H.a.. .__
(38)
Long-range predictivecontrol with a terminalmatchingcondition By comparing the terms in eq. (38) with eq. (34), the polynomials in eq. (34) can be shown to be as follows: R = Tfq-1
(
Hi& + H&.
z
i=nl
(39)
>
v=T(i$l~i+as)
WI
“2
S=
x
F;ai + FIas
(41)
i=nl
where the degrees are 6R = max(nb, nt), 6Y = nt, and SS = na. 4.
NONADAPTIVE
CLOSED-LOOP
ANALYSIS
OF
y-WEIGHTING
4.1. Closed-loop transfer function The equations for a general linear controller in eqs (39)-(41) can be combined with the ARIMAX model in eq. (1) to form a closed-loop system model. Assuming T to be 1, the designed closed-loop transfer function models for y(t) and u(t) are given by
Y(l) = u(t) =
ivq-‘w(t)
+ R&c(t)
RA^A •c q-‘lb _&w(t)
-
Sx(t)
RA^A + q- ‘6s
4.2. Elimination of the integrator A-’ with y-weighting Recall from Section 2.1 that an integrated process model is assumed for GPC strategy. The presence of the integrator, A-‘, in the noise model implicitly increases the overall closed-loop model order by one (McIntosh, 1988). In this section, y-weighting is shown to eliminate the additional model order caused by the integrator and retain offset-free performance at steady-state, so that the overall closed-loop model order using the ARIMAX model with y-weighting is the same as that for the open-loop model. In order to prove the model order reduction by using the steady-state error weighting, a minimum variance control tuning configuration is assumed for simplicity in this section, i.e. nl = 1 n2 = k
For any value of weighting on the steady-state error, co
the closed-loop characteristic polynomial can be derived from eqs (39) and (41) so that
(43)
where x(t) represents a nonstationary or integrated white noise sequence, and is defined as x(t) = &(t)/A. When there is a model-plant mismatch, particularly if the gain of the model is different from the actual gain of the plant due to inexact estimation of polynomials A and B, the presence of the steady-state error weighting term does not cause any offset in the closed-loop system. This property can be seen by analyzing the closed-loop transfer function at steady state. At steady state (i.e. q = l), the closed-loop transfer function in eq. (42) is reduced to
(4)
nu = 1.
o
1293
1 gh- 1Hk + g3.y
R=l+q-
s_
1
(gf-1
+ g:Y)
+gk-&+gsFsy
C&l + dY) The resulting characteristic polynomial ilrom eq. (42) becomes wMyss = RA^A + q-IL& = gk-IIA^A(gh-l
+ g-lHk)
+ q-‘B^F,]
+ gay[A^A(gs + q-‘H,)
+ q-‘BF,].
(45)
Recall from eq. (7) that
V(1)
Y(f) = -w(t) S(1)
and, forj = k and T=
where
1,
"2 V(l)=
C Cri+ a, i-n1
S(1) =
F alF1(I) + a,F,(l). i=“l
(46)
q’(&-gk-l)=H,+$-
Substituting eq. (46) into eq. (45) results in the follow-
ing characteristic polynomial:
Recall from eq. (2), for T = 1, that 1 = E,A^A + q-‘Fj. At steady state, it becomes 1 = FJl) = F,(l). Therefore, V(1) = S(1) =
nZ C oli+ a, i=ml
confirms that y(t) = w(t).
2 %?MY.S
=
8,
+
--&y‘x
(47)
The striking feature of eq. (47) is that the maximum order of the closed-loop model is na (assuming that nb is less than na as is usually the case). Unlike ordinary weighting on control increments in generalized minimum variance control, weighting on the output steady-state error does not incur an additional order in the closed-loop performance. The results in this section are summarized by a comparison of the following three cases:
KUN-YU KWOK and SIRISHL. SHAH
1294
Case 1: ARIMAX model ahead control increments:
with weighted one-step-
Model: Ay(t) = Bu(t - k) + 5B) A Control law: J, = [y(t + k) - w(t + k)]’ + [g.Au(t)]‘. The resulting closed-loop transfer function is Y(t) = Bk- 1Bq- l w(t) + (gk- l&B + s:NW Bk- iB + gf
AA
Although it provides ofieet-free servo control, an extra order in the denominator is incurred by the integrated mode, causing a larger phase shift in frequency domain which requires a lower gain and hence a more sluggish response. Case 2: ARMAX model ahead positional control:
4.3. Stability with y-weighting For any open-loop stable process, the controller with y-weighting tends to be mean-level control as y approaches infinity. This property makes y-weighting similar to a large n2 value in GPC. As n2 or y increases, the closed-loop poIes migrate towards the open-loop pole positions. Therefore, y-weighting is able to provide a stabilizing effect even if there is model-plant mismatch. The stabilizing ability of yweighting is formally stated in the following lemma. Lemma 2: For any open-loop stable process under long-range predictive control with a steady-state weighting term, y, ana’ a stable estimated model f&(q-‘). there exists a finite ualue, ym, such that y > yrn > 0 guarantees a stable closed-loop system. Proof: Although the following proof, for simplicity, is for the ease of nu = 1, similar proofs can be deduced for any size of the control horizon. Recall that the closed-loop characteristic equation from eq. (42) is Ws = RAA+
with weighted one-step-
q-l&.
Substituting eqs (39) and (41) into the closed-loop equation gives
Model:
Hjaj + H,tx,
Ay(t) = Eu(t - k) + e(t)
_
AA
Control law: 32 = [y(t + k) - w(t + k)]’ + [g,u(t)]‘.
\
j=nl
/
In the absence of the integrated mode, the modeling of the noise is assumed to be always stationary. The closed-loop transfer function is y(t)
=
(48)
gk-lBq-%r) + (gk- iJ-hB + g:)<(t) gk- 18 +
g,2A
The characteristic polynomial does not have an extra order. However, this approach is well known for having an oflet problem at steady state.
Then applying eqs (7) and (12) to eq. (48) turns eq. (48) into WS=~A+/iAq-‘[j~,($-Gj)qhj
Case 3: ARIMAX model with one-step-ahead control and steady-state error weighting:
+ (2)asli-ca]. As y increases to infinity, aj diminishes to zero and a, approaches I/g,. Then eq. (49) can be rearranged to the following format:
Model: Ay(t) = J3u(t - k) + 3 A
where
Control law: 33 = CY(~ + k) - w(t + WI2 + Y(Y.., - w.)’ Transfer function: y(t)
=
(gk- 1
+
YB.)4-1B”‘(t)+ (gk-lEkB + w,q-lHs + g:)<(t) a- 18 + w:A
The introduction of weighting on the steady-state error instead of on the control increments eliminates both the extra order in case 1 and the offset problem in case 2.
+ _ .)gj-1 = a,q-’
+&q-2
+ gsy + _. .
1
1
1295
Long-range predictivecontrol with a terminalmatchingcondition The coefficients Si diminish as y increases. According to Liu and Gertler (1987), there is a limiting value 6, > 0 so that the closed loop is guaranteed to be stable for 1Sil -c 6,. Therefore, a limiting value Y,,, 0 exists such that eq. (48) is stable (Kwok, 1992).
puted as follows:
Y 5. PERFORMANCE
TUNING
OF y-WEIGHTING
The important tuning parameters for GPC steady-state error weighting are listed below: l l l l
with
5.1. Strategy # 1: deadbeat approach Clarke and Mohtadi have shown that the equivalent of a stable deadbeat controller for an observable and controllable linear system is obtained by the following selections: n23
2N-
1, A=0
(50)
where N is the order of a state-space model (Mohtadi, 1987). When an input/output model is considered, the configuration of a deadbeat controller requires the following parameter settings: nu = na + 1,
nl = nb + 1,
y tr(G:G,)
= tr(mGrG)
g,2nu(nu+ 1) 2
= tr(mGTG)
2m tr(GTG) ’ =(nu
+ 1)nugf’
(52)
are roughly the sum of deadbeat output response
Since GPC without y-weighting can be reduced to many of the well-known control methods by different selections of the first three tuning parameters, two new tuning strategies are recommended by incorporating y as an “active” tuning parameter into two common control methods: deadbeat and mean-level control. The reason why it is called “active” is that, while keeping the other tuning parameters constant, only y is adjusted to obtain the desired response. Depending on the performance required, the first of the following approaches can be used for underdamped processes to achieve fast response, whereas the second can be used for overdamped systems to increase robustness.
nl = N,
= tr(mGrG)
Using the “formula” in eq. (52), the control increments
minimum output horizon nl maximum output horizon n2 control horizon nu steady-state error weighting y.
nu = N,
tr(GTrG,)
n2 = nu + nl -
1. (51)
The settings specified above are, therefore, fixed at the commissioning stage. Then y is adjusted to tailor the speed of response. In principle, the value of y can vary from zero to infinity. A very large value of y only makes the controller behave more like a mean-level controller. In order to find a y-weighting relative to the deadbeat control action, one can make the magnitude of GFIG. of the control law in eq. (33) “approximately equal” to mGTG, where m is the degree of detuning from deadbeat to mean level. Then the new control action is roughly l/(m + 1) times the original deadbeat action. If the magnitude of the matrix inverse is approximated by the trace of GrG as in McIntosh’s tuning strategy (McIntosh et al., 1991), y can be com-
weighted by l/(m + 1) plus the mean-level response weighted by m/(m + 1). The initial value of y can be chosen by setting m = 1 so that the control performance is about the average of the two responses for commissioning. Because the formulation of y-weighting already has g 3 as a factor in the inverse portion of eq. (33) the closed-loop performance is independent of any gain changes. This is an advantage of y over the use of L-weighting, which often has to be adjusted relative to gz by scaling. Another advantage is that the slowest control response resulting from a large value of y is limited to that of mean-level control; whereas increasing k can indefinitely decrease the response. 5.2. Strategy #2: mean-level approach For typical industrial plant models, a large value of n2 and a value of nu of 1 are usually sufficient to produce satisfactory control performance for openloop stable plants. The control equation in eq. (33) with nu = 1 and without any weighting can be written simply as
It is obvious from the equation that nonminimum phase model effect should be avoided by choosing a large n2 such that C;znl gi_i is strictly positive. Although a default output horizon window of n2 - nl + 1 = 10 has been found to produce robust response, a large n2 corresponding to a 90% rise time of the process is more desirable especially when commissioning a new controller. Depending on the sampling time, of n2 capturing the rise time can vary from 20 to 50, adding a lot of computational load to the controller. When y-weighting is used in conjunction with n2, it becomes the choice for performance tuning. Then the value of n2 can be reduced considerably to capture only the high-frequency dynamics which usually appear during the first “time constant” period. This is similar to the guideline suggested by Maurath et al. (1985) recommending n2 equivalent to the number of sampling periods required for the process open-loop step response to reach 50% of its final value. This is achieved by making n2 approximately equal to the number of sampling units for one overall process time constant. Therefore, the configuation for commission-
1296
KUN-YU
KWOK and SIRISH L. SHAH
ing a process is proposed to be nl = k,
n2 = 5 + k,
nu = 1,
/
=a
(54) where k is the process time delay T is the dominant or overall time constant process l T=is the sampling time 6 Elf., gi-1 + y > 0 ensures that the sign controller is compatible with the sign steady-state gain (to avoid nonminimum effect). l
.
of the
of the of the phase
For example, if the sampling period is chosen by taking one-fifth of a time constant, n2 will be 5 + k, where k is the units of time delay. Starting with a large value of y. the process response can be adjusted by gradually decreasing y until satisfactory results are observed. The control law for this approach is given below:
Fig. 2.
order over- or underdamped processes are discretized by the same factor of the individual time constants, i.e. zl/rl, = ~2/72,, the resulting step response models are exactly equal. Therefore, an approximate value of y and ns which will provide an effect similar to a larger n2 vlaue (with y = 0), can be computed from eq. (57) when nl and T/Z~are given. Figure 2 shows the trajectories of y corresponding to different n2 and ns. For this particular plot, nl is one and the sampling time is taken as one-tenth of the time constant, which is considered to be an appropriate discretization interval for most industrial processes.
Au = elf.,
si-
*Cw(t+ 4 -f( t + 41 + BJCM) -fMl
where ns denotes the reduced output prediction horizon in order to distinguish it from n2 in eq. (53). One can obtain a feel for the range of admissible y values by equating eqs (53) and (55) and finding the values of y corresponding to different n2s. Without any loss of generality, the setpoints in both equations are assumed zero. After some rearrangement, the resulting relationship between y and n2 is given as
Cz,,
!J-lfCt
Y=
i)C~lm,+l ST-1-C;1,+1
+
gi-lf(f
9j=
S*Sj
where the prime denotes a quantity corresponding to the case where the gain is unity. With this substitution, all gain terms in eq. (56) will be cancelled out. What remains is a relation between y and a combination of the model response coefficients with a standardized unity gain:
y, =
+
OCE,, CT?-1
kZC;2nlgi--lf(t+i)-kkS,C;2n,g:-1
It should be noted that the value of y in eq. (56) is dependent neither on the process gain nor on the process time constant (71, but only on the choice of nl, n2, ns, and the ratio r/~,. In order to show that y is independent of the process gain, the gain term is factored out from all step response coefficients, i.e.
CL
g:-lf(t
+
i,~;‘,+, C;f,,
Trajectories of y for different values of a2 for the special case of au = 1.
.
(56)
4. ADAPTIVE
@f-d2 - c12w+I g:- If’@ + Ok;“=., Cd- 3
gf-If’@
+ i) -_KClf.,
It should be noted again that the choice of y is dependent on the discretization factor, 7/7,, but not on 5. The reason is that when two equal-gain, equal-
GPC WITH y-WEIGHTING
One of the areas where closed-loop control has currently received much attention is in the control of biomedical systems (Linkens, 1984; Vozeh and Steimer, 1985; Linkens and Hacisalihzade, 1990). In particular, continuous infusion of a vasodilator for induced hypotension is often required for cardiac patients during and after surgery. The patient’s response to drug intervention is analogous to a complex and nonlinear chemical process subject to unmeasurable disturbance. The infusion of a vasodilator by a com-
(sI_ 1)’
.
(57)
Long-range
predictive control with a terminal matching condition
puter-controlled drug delivery pump (acting as the final control element) is the manipulative variable. Maintaining the mean arterial blood pressure (MAP) (which is the controlled variable) on the setpoint by the drug infusion is the control objective. Adaptive GPC has heen previously applied in this area by Kwok et al. (1991). The present investigation differs from the earlier work because it involves the use of GPC with steady-state error weighting and a control-relevant, long-range identification strategy for the control of MAP on a physiological system. This identification strategy was developed by Shook et al. (1991), designed to provide an optimal trajectory of output predictions instead of a single optimal prediction. Therefore, the identification provides the best match for a trajectory of predictions, which, in turn, is utilized by GfC with y-weighting. The initial tuning of the controller was done according to tuning strategy #2 in the previous section. Details of the theory, system description, and software package are available in Kwok et al. (1991), Shook et al. (1991), and Kwok (1992). 6.1. Experiment This study was approved by the Health Sciences Animal Care Committee of the University of Alberta. Experiments were performed, in an ethics-approved manner, on six healthy anesthetized mongrel dogs of either gender weighing between 16 and 18 kg. The infusion rate of sodium nitroprusside (SNP), which is a potent hypotensive agent, was used as the manipulative variable. Each run began with a controller initialization sequence which consisted of the infusion of SNP (7-15pgkg-‘min-‘) that produced a small hypotensive response. This also constituted the initial probing sequence for the on-line estimation scheme. Thereafter, the performance of the controller
was evaluated for setpoint changes by requesting a computer-controlled setpoint change of - 30 mmHg from the original baseline. Controller performance was also tested during a series of unpredictable disturbances produced by the administration of norepinephrine (NE) infusion to simulate a hypertensive disturbance and another SNP infusion to simulate a hypotensive disturbance. 6.2. Results For brevity, only one series of runs is presented. Detailed results are found in Kwok et al. (1992). For this particular run, the concentration of SNP was 200 pgml-‘. Figure 3 shows the initial probing effect and the subsequent servo response. Since the subject was highly sensitive to SNP infusion, the initial probing of 40 ml h- 1 already caused a large drop in MAP. During computer control, the MAP controller performed favorably because MAP was well within the acceptable boundary of f 5 mmHg around the setpoint. Two maintained disturbances were made consecutively, one in the absence and the other in the presence of closed-loop computer control. Deviations from baseline MAP under these two conditions were compared. In order to ensure that the disturbances in MAP were made from equivalent baseline pressures, MAP was lowered to similar levels by either a constant infusion of SNP (absence of closed-loop control) or by the controller (presence of closed-loop control). Figures 4 and 5 show the MAP trajectories in response to hypertensive and hypotensive disturbances, respectively. These results clearly quantify the performance of our closed-loop controller based on GPC with steady-state error weighting and the long-range identification strategy.
I
100’
a
011 0
1297
2
4
6
8
10
12
14
2
4
6
8
10
12
14
Tie. Fig. 3. MAP
minutes
response to the initial probing signal and setpoint change
I
KUN-YU KWOK and SIRISH L. SHA%I
1298
Without
closed-loop
control
+-
1-+ With closed-loop
control
601 0
10
20
30
40
50
60
I 70
405 0
10
20
30
‘IO
50
60
70
Time. minutes Fig. 4. MAP response
Without
to norepinephrine
closed-loop
control
+
causinghigh blood pressure.
(NE) infusion
1+ With closed-loop
control
90 85 _________________________________-________________________________ ____------_-
6.5 0
5
10
15
“0
5
IO
15
20
25
30
35
40
20
25
30
35
40
Time, minutes Fig. 5. MAP response
to another
nitroprusside
(SNP) infusion
7. CONCLUSIONS
A terminal matching condition is introduced into LRPC as a weighting on the square of the error between the predicted process steady-state value and the setpoint. The final control law only minimizes the square of the prediction error over a small future prediction horizon and the weighted terminal condition. A computational scheme is also derived for both parametric and convolution models so that, for the
causing
further blood pressure
drop.
ARIMAX model, the steady-state value is calculated without iterating the Diophantine identities and, for the convolution model, a complete set of step response coefficients is no longer required. The latter model type only relies on a few dynamic response coefficients plus the steady-state gain to make future predictions. The closed-loop analysis also shows that the weighted square of steady-state error in the control law not only provides a stabilizing effect as a large
Long-range predictivecontrol with a terminalmatchingcondition prediction horizon but also eliminates the extra closed-loop model order incurred as a result of the integral model. Two tuning strategies based on deadbeat and mean-level concept are provided for using the terminal matching condition with generalized predictive control, which has been successfully applied to the control of mean arterial pressure.
NOTATION
a,
Atq-‘1 bf
ml--‘) ctq-‘1 et Ej(q- '1 ftt+ i) f f,
Fj(q- ‘1 Si
Gj(q- ‘1 G,tq G
‘1
G hi h R;;q-
‘)
J k nl n2 na ns nb nt nu P --I 4 R(q- ‘)
:(4-l) T(q_‘) T u (0 U
approximation error in convolution model reference value in convolution model process model denominator polynomial ith coefficient of the numerator polynomial in G, process model numerator polynomial noise model polynomial ith coefficient of E,(q-I) quotient polynomial from jth-step-ahead Diophantine identity TJA ith open-loop process output prediction assuming no future control action vector off@ + i), where i is from nl to n2 vector off(s) remainder polynomial from jth stepahead Diophantine identity T/A ith coefficient of G,(q- ‘) quotient polynomial from jth-step-ahead Diophantine identity E,B/T transfer function of any process matrix of step response elements lower triangular matrix containing only &s ith impulse response coefficient ith coefficient of H,(q- ‘) remainder polynomial from jth-stepahead Diophantine identity EjBfA performance index time delay in the number of sampling intervals minimum output prediction horizon maximum output prediction horizon order of A(q- ‘) maximum output prediction horizon less than n2 order of B(q ‘) order of T(q- ‘) control horizon first-order coefficient in the denominator of G, backward shift operator polynomial in the linear form of GPC control law value at steady state polynomial in the linear form of GPC control law observer design polynomial model horizon for convolution model control input vector of incremental control inputs
v(q-‘) w(t) W x(r) v(t) Y
1299
polynomial in the linear form of GPC control law setpoint vector of setpoint trajectory noise input process output vector of output prediction values
Greek symbols ith element of the first row of ai [GTG + A + GfrG,] - ‘GT sum of pi from nl to n2 ith element of the first row of ; [GTG + A + G,TrG,] - ‘G:‘lcontrol steady-state weighting Y a critical steady-state weighting which Ya makes closed-loop control stable Idiagonal matrix containing control steady-state weightings 6 degree of the polynomial 4 A
1 A 6(t) Z r.
fth coefficient of L@ differencing operator
(= 1 - q-l) control weighting diagonal matrix containing control weightings uncorrelated random sequence of zero mean overall time constant sampling time
Superscripts
T
estimated value value corresponding to unity gain transposition
Subscripts /”
MV
S Others v a
d-step time delay signal filtered by l/T(q-‘) minimum variance steady state
characteristic equation a portion of %?s REFERENCES
K. J. and Wittenmark, B., 1973, On self tuning regulators.Automatica9, 185-199. Auslander,D. M., Takahashi,Y. and Tomizuka. M.. 1978. Direct diaital ntocesscontrol:nracticeand algorithmsfor micropro&& application.P&. IEEE 66, r99-208. Chesna, S. A. and Ydstie, B. E., 1986, Self tuning and adaptive control of chemical processes, in Adoptive and Learning Systems: Theory and Applications (Edited by Astrcm,
K. S. Narendra), pp. 119-130. Plenum Press, New York. Chesna. S. A. and Ydstie. B. E.. 1988. Controlline nonminimum chase svstems with nredictive controllers. in Workshop ‘on Ado&e Control’ Strategies for Industrial Use (Edited by S. L. Shah and G. A. Dumont), __ PD. 204-231. Lodge Kananaskis, Alberta. Clarke, D. W. and Gawthrop. P. J., 1975, Self-tuning controller. IEE Pmt. D 122.929-934. Clarke, D. W. and Mohtadi, C., 1989, Properties of generalized predictive control. Automatica 25, 859-875.
1300
KUN-YU KWOK and SIRISH L. SHAH
Clarke, D. W., Mohtadi, C. and Tuffs, P. S.. 1987. Gentralised predictive control-Part I. The basic algorithm. Automatica 23, 137-148. Cutler, C. R. and Ramaker, B. L., 1980, Dynamic matrix control-a computer control algorithm, in Proc. JACC, San Francisco. Francis. B. A. and Wonham. W. M.. 1976. The internal model principle of control thkory. Au&n&a 12457467. Garcia, C. E. and Morari, M., 1982, Internal model control. 1. A unifying review and some new results. Ind. Engng Chem Process Des. Den. 21, 308-323. Garcia, C. E., Prett, D. M. and Morari, M., 1989, Model medictive control: theorv and oracti-a survey. Auton&a 25, 335-348. Kwok, K. Y., 1992, Long-range adaptive predictive control: extensions and applications. Ph.D. thdsis, University of Alberta. Kwok, K. Y., Mutha, R K., Shah, S. L., Clanachan, A. S. and Finegan, B. A., 1991, Constrained long-range adaptive control of arterial blood pressure.. J. predictive Adaptive Conrrol Signal Processing 5, 363-374. Kwok; K. Y., Shah, s. L., Clanachan, A. S. and Finegan, B. A., 1992, Evaluation of a long-range adaptive predictive controller for computerized drug delivery systems, in Proceedings o/the IFAC International Symposium on Adaptive ControI and Signal Processing, Grenoble, France. Lambert, E. P., 1987, Process control applications of long range prediction. D.Phil. thesis, University of Oxford. Li, S., Lim, K. Y. and Fisher, D. G., 1989, A state space formulation for model predictive control. A.I.Ck.E. J. 35, 241-249. Lim, K. Y., 1988, Multivariable optimal constrained control aleorithm IMOCCAI. M.Sc. thesis. University of Alberta. Linkins, D. k., 1984, Computer codtrol in biimedicine, in Real-Time Computer ControL (Edited by S. Bennett and D. A. Linkens),Chap. 14. P. geregrin&London. Linkens, D. A. and Hacisalihzadc, S. S., 1990, Computer control systems and pharmalogical drug administration: a survey. .I. Med. Engng Technol. M(2), 41-54. Liu, K. and Gertler, J., 1987, On-line stabilization of adaptive controllers by detuning in a supervisory framework, in Proceedings of the American Control Confwence, Minneapolis, MN, pp. 194-200.
McIntosh, A. R., 1988, Performance and tuning of adaptive generalized predictive control. M.Sc. thesis, University of Alberta. McIntosh, A. R., Shah, S. L. and Fisher, D, G., 1991, Analysis and tuning of adaptive generalized predictive control. Can. J. them Engng 69,97-l 10. Maurath, P. R., Mellichamp, D. A. and Seborg, D. E., 1985, Predictive controller design for SISO systems, in Proceedings of the American Control Conference, Paper FP4, pp. 15461552. Mayne, D. Q. and Michalska. H., 1991, Model predictive control of nonlinear systems, in Proceedings ofthe American Control Conference, FA8, pp. 2343-2348. Mohtadi, C., 1987, Advanced self-tuning algorithms. D.Phil. thesis, University of Oxford. PaDadoulis. A. V.. Tsilixiannis. C. A. and Svoronos. S. A.. i987, A cautious self-&ning’controller for chemical pro: cesses. A.1.Ch.E. J. 33, 40149. Rawlings, J. B. and Muske, K. R., 1992, The stability of constrained receding horizon control. IEEE Trans. Automat. Control (submitted). Richalet, J., Rault, A., Testud, J. L. and Papon, J., 1978, Model predictive heuristic control: applications to industrial processes. Automatica 14, 413428. Shook, D. S., Mohtadi, C. and Shah, S. L., 1991, Identification for long range predictive control. IEE Proc. D 138,75-84. Shook, D. S., Mohtadi, C. and Shah, S. L., 1992, A controlrelevant identification strategy for GPC. ZEEE Trans. Automat. Control 37, 975-980. Sripada, N. R. and Fisher, D. G., 1985, Multivariable optimal constrained control algorithm (MOCCA): Part 1. Formulation and aDDlication. in Proceedinas of the Internatonal Conferench bn Industrial Process VMideling and Control, Vol. 1, Hangzhou, China. Takahashi, Y., Tomizuka, M. and Auslander, D. M., 1975, Simple discrete control of industrial processes (finite time setting control algorithm for single-loop digital controller). Trans. ASME, pp. 354-361. Vozeh, S. and Steimer, J. L., 1985, Feedback control methods for drug dosage optimisation. Clin. Pkarmacokin. 10, 457476.
Chnntcof E#a&eering Scimce, Vol.
49. No. 9, p,,. 1287%13W, ,994 Copyright 8 1994 Elsevicr Science Ltd Printed in Gsal Britain. All rights rcsorvcd OOIX-2509,Y.t 56.00 + 0.00
LONG-RANGE
PREDICTIVE CONTROL WITH A TERMINAL MATCHING CONDITION KUN-YU
Department
KWOK
and SIRISH L. SHAH’
of Chemical Engineering, University of Alberta, Edmonton,
Alberta, Canada T6G 2G6
(First receiued 27 December 1991; accepted in revisedfirm 6 September 1993) Abstract-The long-range predictive control strategy is modified by augmenting the generalized predictive control law with a terminal matching condition. This condition is formulated as a weighting term on the output steady-state error. The final control law only minimizes the square of prediction errors over a small future prediction horizon and at steady state. It is applicable to either transfer function or convolution modeling. For transfer function models, the modified approach reduces the computational load of an adaptive controller by mimicking a large prediction horizon with a smaller one plus a weighting at the steady state. For convolution models, only a few initial step response coefficients and the steady-state gain are required to formulate the new long-range predictive controller. The nonadaptive closed-loop analysis reveals that this new steady-state error weighting term (i) reduces the additional order otherwise required due to inclusion of an integrator for offset elimination, and (ii) provides the stabilizing effect of a large predictive control horizon. A closed-loop adaptive control system based on this modified control law has been successfully applied to the control of mean arterial blood pressure. 1. INTRODUCTION
for minimum variance (bistriim and Wittenmark, 1973) and generalized minimum variance control (Clarke and Gawthrop, 1975) were aimed at providing tight control of processes by minimizing the variance of the process output predictions. A key factor in the successful use of these methods is correct prior selections of the model order, n, and time delay, k, such that accurate k-step-ahead predictions are made. Poor selection of k or a varying time delay affects the control performance or even leads to instability. Improper selections of model order may result in similar problems although overspecified model orders may not be as sensitive as the uncertainties in time delay. Recent developments in long-range predictive control (LRPC) which attempt to minimize the sum of prediction errors over a user-specified time horizon provide an alternative solution to the abovementioned drawbacks (Richalet et al., 1978; Cutler and Ramaker, 1980; Clarke and Mohtadi, 1989). Since the strength of the LRPC algorithm lies in the longrange predictions of process outputs and control actions, these two predictions which have a direct influence on stability and dynamic response in any longrange predictive controller (Garcia and Morari, 1982) naturally become the tuning parameters. The selection of the output prediction horizon is less influential than the choice of the time delay or model order. However, the long-range predictive nature necessarily makes the computation load of an adaptive controller very heavy for each new control action to be calculated. On the other hand, a relatively short output The
early
methods
prediction horizon, of which the minimum variance is the extreme case, may result in vigorous control actions or even unstable control. Moreover, the longrange control algorithms such as MAC (Richalet et
‘Author
to whom correspondence
should be addressed.
al., 1978), DMC (Cutler and Ramaker, 1980), and MOCCA (Sripada and Fisher, 1985). which employ a convolution model as the internal model (Garcia et al., 1989) to encompass most of the process dynamics, require as many as 50 step or impulse response coefficients, a size not economically sound for adaptive control. This paper proposes a modified approach to the LRPC scheme which introduces an output terminal condition in the minimization of a multistep prediction cost function. The terminal matching condition is applied in the form of weighting on the square of the error between the predicted process steady-state value and the desired steady-state target or setpoint. This weighting combined with a relatively short output prediction horizon, in the context of adaptive and nonadaptive control, is able to resolve the trade-off between economic computation and long-range robustness. The concept of steady-state error weighting is depicted in Fig. 1. An ordinary LRPC algorithm considers process output predictions over a large horizon from nl to nlonp.The modified approach proposed here would consider a relatively shorter horizon from nl to nsho,,plus the term at steady state so as to ensure the terminal matching of process output to the setpoint. Different forms of terminal conditions have appeared in the literature. However, most of the applications have occurred in the context of an infinite LQ controller with terminal state constraints for the improvement of stability (Chesna and Ydstie, 1988; Mayne and Michalska, 1991; Rawlings and Muske, 1992), and in a self-tuning controller with a cautious control element during a poor identification period (Papadoulis et al., 1987), but the concept of a terminal matching condition formulated as a weighted square of the steady-state error and the relative weight of this term as an on-line parameter has not been investigated in the literature. This concept proposed in this
KUN-YU KWOKand
1288
SIRISH
L. SHAH
(4) There are recommendations of two performance tuning strategies for LRPC with steadystate error weighting: one for fast response, the other for robust response. (5) The successful experimental evaluation of this (adaptive) control strategy in the control of mean arterial blood pressure.
Fig. 1. Long-range predictive control weighting.
with
steady-state error
paper minimizes a finite-time objective function augmented by a weighted infinite-time element (for the terminal condition), and is therefore functionally different from infinite-horizon LQ strategies, i.e. the proposed controller is simultaneously able to accommodate finite-time and infinite-time objectives. In this study, the main contributions, in the context of adaptive and nonadaptive control of processes, are:
(1) The incorporation
of a weighted steady-state term as a tuning parameter into the performance criteria of LRPC. Philosophically speaking, the idea of a terminal matching condition is intuitively appealing and easy to relate to by plant engineers. The concept of terminal matching can be readily explained when one considers the analogy of a car trip from point A to destination B. Whilst the performance of the vehicle over a finite and continuously moving window of say 200 m is important, equally important is the final destination B. The steady-state weighting gives us the capability of incorporating this condition in the design philosophy. (2) The incorporation of such weighting in the optimization criterion obviously requires steadystate output predictions. In order to reduce the heavy computational load as used in conventional steady-state prediction techniques, technical results presented in Section 2 provide a method to make steady-state predictions with fewer iterations. The method allows steadystate output predictions to be computed directly from the ARIMAX model without iterating to a large output prediction horizon. The same method can also apply to the convolution model in that a complete set of step response coefficients is not required. (3) The extra closed-loop model order incurred as a result of the integral model is eliminated by using the steady-state error weighting.
This paper first presents the procedure for on-line model-based finite-horizon and steady-state predictions of the process output responses based on ARIMAX or convolution models. Then the control law with the steady-state error weighting and its closed-loop analysis are discussed in the context of nonadaptive control, followed by the performance tuning strategy. The last section of this paper documents the results of the application of this modified LRPC strategy in the control of mean arterial blood pressure. 2.
PROCESS MODELING
Most model predictive control techniques are based on the optimization of a quadratic objective function involving the error between the setpoint and the predicted outputs. The distinguishing feature among different techniques is the model type and structure for making the predictions. The most commonly used model types are discrete-time transfer functions describing the input-output relationship, or discrete-time convolution models in the form of step or impulse response coefficients. The use of both model types for finite-horizon plus steady-state predictions is considered in the following sections. 2.1. Transfer junction model The process model is assumed to be an autoregressive integrated moving-average model with an exogenous input (ARIMAX):
A(q-‘)y(t) where
= &q-‘)u(t
- 1) + T(q+$
(1)
A=l-q-’ B^= Bq-‘+’ 4 -’ = backward shift operator.
Having the integrator in the noise model also represents a class of nonstationary dynamics such as random steps at random times, i.e. random-walk or Brownian motion type disturbances. Note that this integral mode serves as an internal model of the (steptype) disturbances and setpoint changes (Francis and Wonham, 1976). To derive ajth step-ahead predictor from eq. (i), the following Diophantine identity is considered: 7”= EjA^A + q-jFj.
(2)
Multiplying eq. (1) by AEjq’ gives E,AAy(t
+ j) = E,l?Au(t + j - 1) + EJTC; (t + j). (3)
Long-range predictivecontrol with a terminalmatchingcondition
1289
Substituting eq. (2) into eq. (3) results in the following jth step-ahead model:
and substituted into eq. (5), the result is
Ty(t + j) = E,BAu(t + j - 1) + F,y(t) -+ E,T<(t + j).
TAA
(4) Future control actions are separated from the known ones by substituting the following identity into eq. (4): (5)
Eji% = G,T + Hjq-j
where the elements of Gj now correspond to the first 1 open-loop step response coefficients. The optimal jth step-ahead predictor for j 2 1 becomes j -
(T - q-‘Fj)6
B q-‘Fja ,-,=G,++ AA
(6)
asj-bco.
qj($-G,)=++-%&
(7)
TAA l? -=Gj+q-j@+g) AA
(11)
For sufficiently large j, both Fj and Hi converge to F, and H,, so eq. (12) yields eq. (7). a (determination of F’, and H,): Recall that polynomials Fj+ 1 and Hj+ 1 can be defined as Remark
T=
Ej+,A^A + q-j-‘Fj+,
Ej+~~=Gj+~T+q-j-‘Hj+~.
04)
ejA^A = Fj - q-lFj,I ejB = gjT + q-‘Hj+
(15) 1 - H,
Proof: Rearranging eq. (2) gives = Ej + 9.
(8)
The denominator polynomial A^A has one root on the unit circle and nn roots inside the unit circle (na is the order of A^).According to the Taylor seriesexpansions of functions around the point on the unit circle for j approaching infinity, the tailing coefficients in the quotient polynomial Ej converge to the constant value T (l)//l^( 1) with order approaching infinity. Hence, the remainder q-jFj does not vanish but converges to a finite-length polynomial q-“F,, where Fj --+ F,
as+
CO.
E, =
?ZA
- H,
(18)
where T!I) es=A^o g.=-.
&l) A(l)
Therefore, F, is found from eq. (17) as F, = e,A^. Equation (18) becomes H,A = g,T - e,b.
When Ej in eq. (8) is rearranged as T - q-‘F,
(17)
ejg = e,6 = g,T + q-‘H,
z
(16)
where F and H are of degree na and max(nb, nt) - 1, respectively, and ej and g, are the last coefficients in Ej+ , and Gj+ 1. For sufficiently large j, F, H, and the tailing coefficients in E and G converge so that ejiA=e,AA=F,-qq’F,
T
(13)
Subtracting eqs (2) and (5) from eqs (13) and (14), respectively, gives
Lemma 1: Given two polynomials T and A with their roots all inside the unit circle. For sufJiciently large j, the
polynomials Fj in eq. (2) and Hj in eq. (5) conuerge to two Jinite constantpolynomials (de&d QSF, and H,), and the combination ofthe two identities in eqs (2)and (5) becomes
(IO)
Equation (10) indicates that the second Diophantine identity in eq. (5) also has one root on the unit circle and the rest inside the unit circle. Therefore, the same conclusion applies to eq. (5), i.e. H, converges to a finitelength constant polynomial as j approaches infmity. The last part of the lemma is verified by rearranging eq. (10) as follows:
$(f+jlt)=GjAu(t+j-l)+HjAuf(t-l)+Fjy,-(t) where uJ and y, are values filtered by l/T. Two methods of calculating the optimal steady-state prediction are avaiiable. The first method is simply by iterating eq. (6) for large values of j until $(t + jl t) converges. Each iteration requires the solutions of solving the Diophantine identities 2 and 5 with order j + degree (Fj, H,). As the prediction horizon jincreases, the computational load for the solutions also increases at every sampling interval. This tedious method is of course not recommended because the heavy computational load would defeat the purpose of our modified control strategy. Instead, a new and simpler method is described below to calculate the optimal steady-state output prediction at every sampling time. The following technical result in the form of a lemma is required for this simple prediction method.
=Gj+!+_
(19)
H, can be computed by comparing the coefficients in eq. (19). Thus, (9)
H, = h,,,, + h,, ,q-’
+
* - . + hS,nhqpmk
1290
KUN-YU KWOKand SIRISHL.
where nh = max(nb - 1, nt h‘.i = esjig,
1)
bj - g* ,=gl
ti*
q
Now applying the result in lemma 1 to the jth step ARIMAX model in eq. (6) yields (assuming unity time delay for simplicity)
= lim [G,Au(t tj j-02
- 1) + HjAu,(t
-
j-m
1)
-
y(t) = i giq-‘Au(t i- 1
lim [g,_,Atr(t) g,_,,.Au(t
+ gj_,Au(t
+ 1) +
lim 9(t +jl
a, = i=$+ 1 giq-‘Au(t
- 2) + goAu(t + j - l)] - k) + F,y,(t).
(20)
+ 1 - k)
i=
= m I t) = lim i-g_,Au(t)
+ g,._*Au(t + 1) +
i-m
+ IIU-
.. .
1)] + H,AU,(t - k)
+ F,y/(r) = [g%Au(t) + gSAu(r + 1) +
.. .
+ g,Au(t + nu - l)] + H,Aur(t - k) + FSY,0) = g,[Au(t) + Au(t + 1) +
-. .
+ Au@ + nu - 1)] + H,Au,(t - k) + Fa,-(t). (21) As a result, the optimal steady-state predictor based on an ARIMAX model is given as $(si t) = gs f Au(t +j i=i
1) + H,Au/ (t - 1) + FSyr(t). (22)
The first term on the right-hand side of eq. (22) is the current and future forced response, whereas the remaining terms corresponding to the impact of past
(truncation error).
T is usually chosen so that the truncation error is less than 1% of the step-input steady-state gain. Depending on the sampling rate, the number of coefficients may be as many as 50 to encompass 99% of the step response dynamics. The process model becomes
j(t) = f: g;q-‘Au(t
t)
+ gj-..Au(t
+ 1 - k) + ae + a,
T is the model horizon
For practical implementation, the projected control increments after a control horizon of nu are assumed to be zero (Cutler and Ramaker, 1980). Therefore, as gi converges to gs for increasing i, eq. (20) reduces to
j-m
(23)
where
. -
+ nu - 1) + * - ’
+ g,Au(t tj + H,Auf(t
+ 1 - k) + ass
In practice, the dynamic model is commonly assumed by a finite number of coefficients, i.e.
1)
FsY_r@)
j-.n
-t-
y(t) = 5 giq-‘Au(t 1=1
hi = Agi.
= lim [GjAu(t + j - l)] + H,Au,(t
=
2.2. Conuolution model The convolution model structure is usually represented by a step response model:
where the step response coefficients, gi, tend to g. when i tends to infinity and a,, is a reference value. It is readily turned into an impulse response model by differencing the step response such that
+ ~jY.f(t)l
+
SHAH
inputs and outputs on the prediction of the infinite time value of the process output.
+ 1 - k) + a_.
1
(24)
Its advantage over a parametric model is the capability to model unusual dynamic behavior which cannot be well represented by a reduced-order model. Moreover, many chemical processes are inherently infiniteorder and therefore a simple first- or second-order transfer function may introduce structural mismatch that is unfavorable to predictive control (Shook et al., 1992). However, on-line identification of such a large number of parameters is a difficult task. The process model in eq. (24) can be arranged into different predictive representations such as recursive form or state-space formulation (Garcia and Morari, 1982; Lim, 1988; Li et al., 1989). However, all of them assume a linear superposition of previous inputs to the total output such that (for unity time delay, i.e. k= I) jJ(t + jlt) = 5
giAu(t tj
- i)
i=l T +
,=F+, g&G
+j
- i) + ass,
(25)
The output comprises three terms: future control effects, past control effects, and a reference plus residual term. Equation (25) shows that even with steady-state weighting a complete step response sequence of up to
Long-range
T intervals is required regardless of the choice of output prediction, j. However, an approximation alternative, borrowed from Auslander et al. (1978), can be used so that a finite number of initial step response coefficients plus the steady-state gain are sufficient for long-range predictions. This approximation recognizes that most process responses are dominated by a single time constant at least in the low-frequency systems, the lowregion; i.e., for overdamped frequency dynamics can be well approximated by a first-order model. Consider that a unit step response from a stable second- or higher-order process is often characterized by an S-shaped curve, the low-frequency portion can be approximated by an exponential function from the nth point to the steady state. In terms of impulse response, the approximation is written as G,=
hlq-’
+ h2q-’
+ /&q-n + h,+Iq-n+’ =h,q-‘+h2q-2+
3. CONTROL
ALGORITHM
3.1. Derivation The design of a long-range, model-predictive controller is based on the predicted behavior of the process over the prediction horizon. The control objective is to minimize the errors between the trajectory of predicted values and setpoint over a horizon. The following quadratic function defines the performance index to be minimized to achieve the control objective:
r,(j) CP(t + j) - w(t +iW
.
+
...
to be A^= (1 - pq- ‘) such that the ARIMAX modeling conforms to an output error prediction type which is implicitly assumed in DMC type control algorithms that use the convolution model type.
JW, G n4 4 7) =
. . + hn_-lq-“+l
+
1291
predictive control with a terminal matching condition
+ F
hnq-”
+h,_lq-m+l+T
d(j)[Au(t
+j-
l)]’
ci -
1) - WI*
j=l
1 -pq(26) where p. being the rate of exponential decay, is a discrete characterization of the dominant time constant. If the steady-state gain of the process is known, for example, from steady-state design equations (Chesna and Ydstie, 1986), p is determined so that the model gain in eq. (26) is equal to the process gain:
h, LL - J$‘r;
hi’
(27)
The model in eq. (26) can be rearranged into a transfer function from G =b;q-LJ;b;q-2+ ... P 1 - pq-1
+b:q-” (28)
where b; = hI fori=2-n.
bi=hi-phi-1
A stable overdamped model is guaranteed as long as i= 1 It has been shown that the modeling error that stems from the exponential decay approximation can be made arbitrarily small by increasing the sampling period for a fixed n coefficients or vice versa (Takahashi et al., 1975). By transforming a “ limited” convolution model (a few step response coefficients and gain) into a transfer function form, the ARTMAX modeling procedure in Section 2.1 can be applied to make future predictions_ In doing so, the user-defined polynomial T is selected
5 r(j)CPM~ j=l
(29) where l
. l
G,(I) = 4s p=l-
+
l l l l
nl is the minimum output prediction horizon n2 is the maximum output prediction horizon nu is the control horizon such that Au@ + j) = 0, Vj >, nu y,(j) is an output weighting sequence d(j) is a control weighting sequence v(j) is a steady-state error weighting sequence s denotes a value at steady state.
The first two terms on the right-hand side form the standard generalized predictive control (GPC) objective. The last term corresponds to the additional terms penalizing the squares of steady-state predicted error. The summation from 1 to nu is required because, at each control interval j, where j = 1-nu, the steady-state prediction involves the sum of the first j consecutive control actions. It has been suggested that, for robust control, n2 should encompass the whole rise time of a process (Clarke et al., 1987), which often requires as many as 50 terms. Now with steadystate error weighting in place, n2 can be reduced to cover only the high-frequency dynamics. The process predicted values are required in eq. (29), from j = nl-n2 for j(t + j) and j = 1-nu for 9(slt + j - 1). To derive a control law from eq. (29), the predictive equations for the ARIMAX model in eqs (6) and (22) are used. The ARIMAX model can be expressed in the following compact vector/matrix equation (assuming unity delay): Y=GU+f
(30)
KUN-YU
1292
KWOK and SIRISH L. SHAH
GPC control law. Therefore, eq. (33) is referred to in short as “GPC with y-weighting”. For simplicity, the output weighting is assumed to be unity; thus, r, reduces to an identity matrix in the sequel.
where Y = [$(i + nl) P(t + nl + 1). -.
iL-1 9n1
..*
g, 0
0
*”
81
0 ..
Bo
-.
G:
9(t + n2)]’
.‘.
0
.-*
0
-.
90
-.-
CL2-
-
_..
clnl--I
9n2-2
II”
U = [Au(r) Au@ + 1) _. . Au@ + nu - l)lr f=[f(t+nl)f(t+nl+l) f(t +i) = H,Au,(t
...
f(t-tn2)]’
- 1) + Fj_Y,(t)
and Y, = GJJ + f,
(31)
where Y, = [J+Ic) j(slt
. Jvslt + nu - 1)lT
+ 1) .
rg. 0
3.2 Control law in linear form Although the control law in eq. (33) calculates the future control action from r = 1 to nu, only the current control action, i.e. Au(t), is implemented. At the next sample, the entire long-range minimization is repeated once again to calculate the next control action. This receding-horizon policy produces better control performance as the control law makes use of updated plant information. It has been shown that, in the absence of disturbance or model change, the control sequences in successive times are the same (Clarke and Mohtadi, 1989). Because of the receding-horizon policy, only the first element in vector U of eq. (33) is required. To enable closed-loop analysis of this GPC control law, the receding-horizon control law is expressed in a linear, difference equation form (Lambert, 1987; Mclntosh et al., 1991). A general linear structure of a controller with a setpoint and measurement variable is represented by
01
‘..
RAu(c) = VW(r) - Q(r) where R,V, and S are polynomials shift operator. Now define a vector
f&u x 1) = [l 1 . * . 13’CH.h, it - 1)+ F,Y, (dl. The equivalent of eq. (29) using vector/matrix notation is J = [Y - W]VJY
(32)
[GrG
W = [w(t
+ nl)
. - * y(m)]
w(r + nl + 1)
W,(nuxl)=
[1 1
..
. . .
w(r + n2)]r
+ A + Gfl-G,]
+ Gfr(W,
-
fJ1
Au(t) =
l]‘w(s).
- ‘[GTl-,(W
1~~~
w(t)
After substituting eqs. (30) and (31) into eq. (32), the minimization of the cost function yields the following control equation: U = [GTl-,G
(35)
.*a
+ (34)
+A2
+ A + GTl-G,] - ‘G:l-.
Then the current control action can be written as follows (assuming all future setpoints are the same as the present one),
A = diag [A(l) K(2) . . . A(m)]
r = diag [y(l) y(2)
a.21
...
where #Iicorresponds to the elements of the first row vector in the matrix
where r,. = diag Cry(nl) y,@l + 1) . . . y&2)3
in the backward
which denotes the first row of the matrix [GrG A + G:T&]-‘Gr, and a variable a,=Al+&+l
- W] + LJTAU
Wm-CY, - w.1
+ CY, -
Cm,, a,l+i
(34)
- f)
’ . . an2 a,]
- F.1 Y/N -F,i+,y#
44 x
an1 +1
- H., Au,+ - 1) -Kr+lAu,-(t-
1)
: - F,ay&)
w(t)
- H,,Au,(t
- 1)
- H,Au,(t i
- 1) I
- F,y, I (El
I w(t)
(33)
where GTT,G and GTTG. are both matrices of dimension nu x nu. The dynamic matrix, G, contains all the step response coefficients arranged in a lower triangular structure. The term G,TrG, is always of fun rank. It guarantees a nonsingular inverse even if the matrix G’T,,G is ill-conditioned due to large time delays or too short an output prediction horizon. When r is set to zero, eq. (33) reduces to the basic
(37) and, in compact Au(t) = w(r) 2
*=“I
- Aq(t
form, ai - yr(r)
- 1) 2
g F,a, I=“, ff{Kj
i=n, +
w(r)a. . , L-
v,(r)F.a. _, , _ - - Au&,_ - 1)H.a.. .__
(38)
Long-range predictivecontrol with a terminalmatchingcondition By comparing the terms in eq. (38) with eq. (34), the polynomials in eq. (34) can be shown to be as follows: R = Tfq-1
(
Hi& + H&.
z
i=nl
(39)
>
v=T(i$l~i+as)
WI
“2
S=
x
F;ai + FIas
(41)
i=nl
where the degrees are 6R = max(nb, nt), 6Y = nt, and SS = na. 4.
NONADAPTIVE
CLOSED-LOOP
ANALYSIS
OF
y-WEIGHTING
4.1. Closed-loop transfer function The equations for a general linear controller in eqs (39)-(41) can be combined with the ARIMAX model in eq. (1) to form a closed-loop system model. Assuming T to be 1, the designed closed-loop transfer function models for y(t) and u(t) are given by
Y(l) = u(t) =
ivq-‘w(t)
+ R&c(t)
RA^A •c q-‘lb _&w(t)
-
Sx(t)
RA^A + q- ‘6s
4.2. Elimination of the integrator A-’ with y-weighting Recall from Section 2.1 that an integrated process model is assumed for GPC strategy. The presence of the integrator, A-‘, in the noise model implicitly increases the overall closed-loop model order by one (McIntosh, 1988). In this section, y-weighting is shown to eliminate the additional model order caused by the integrator and retain offset-free performance at steady-state, so that the overall closed-loop model order using the ARIMAX model with y-weighting is the same as that for the open-loop model. In order to prove the model order reduction by using the steady-state error weighting, a minimum variance control tuning configuration is assumed for simplicity in this section, i.e. nl = 1 n2 = k
For any value of weighting on the steady-state error, co
the closed-loop characteristic polynomial can be derived from eqs (39) and (41) so that
(43)
where x(t) represents a nonstationary or integrated white noise sequence, and is defined as x(t) = &(t)/A. When there is a model-plant mismatch, particularly if the gain of the model is different from the actual gain of the plant due to inexact estimation of polynomials A and B, the presence of the steady-state error weighting term does not cause any offset in the closed-loop system. This property can be seen by analyzing the closed-loop transfer function at steady state. At steady state (i.e. q = l), the closed-loop transfer function in eq. (42) is reduced to
(4)
nu = 1.
o
1293
1 gh- 1Hk + g3.y
R=l+q-
s_
1
(gf-1
+ g:Y)
+gk-&+gsFsy
C&l + dY) The resulting characteristic polynomial ilrom eq. (42) becomes wMyss = RA^A + q-IL& = gk-IIA^A(gh-l
+ g-lHk)
+ q-‘B^F,]
+ gay[A^A(gs + q-‘H,)
+ q-‘BF,].
(45)
Recall from eq. (7) that
V(1)
Y(f) = -w(t) S(1)
and, forj = k and T=
where
1,
"2 V(l)=
C Cri+ a, i-n1
S(1) =
F alF1(I) + a,F,(l). i=“l
(46)
q’(&-gk-l)=H,+$-
Substituting eq. (46) into eq. (45) results in the follow-
ing characteristic polynomial:
Recall from eq. (2), for T = 1, that 1 = E,A^A + q-‘Fj. At steady state, it becomes 1 = FJl) = F,(l). Therefore, V(1) = S(1) =
nZ C oli+ a, i=ml
confirms that y(t) = w(t).
2 %?MY.S
=
8,
+
--&y‘x
(47)
The striking feature of eq. (47) is that the maximum order of the closed-loop model is na (assuming that nb is less than na as is usually the case). Unlike ordinary weighting on control increments in generalized minimum variance control, weighting on the output steady-state error does not incur an additional order in the closed-loop performance. The results in this section are summarized by a comparison of the following three cases:
KUN-YU KWOK and SIRISHL. SHAH
1294
Case 1: ARIMAX model ahead control increments:
with weighted one-step-
Model: Ay(t) = Bu(t - k) + 5B) A Control law: J, = [y(t + k) - w(t + k)]’ + [g.Au(t)]‘. The resulting closed-loop transfer function is Y(t) = Bk- 1Bq- l w(t) + (gk- l&B + s:NW Bk- iB + gf
AA
Although it provides ofieet-free servo control, an extra order in the denominator is incurred by the integrated mode, causing a larger phase shift in frequency domain which requires a lower gain and hence a more sluggish response. Case 2: ARMAX model ahead positional control:
4.3. Stability with y-weighting For any open-loop stable process, the controller with y-weighting tends to be mean-level control as y approaches infinity. This property makes y-weighting similar to a large n2 value in GPC. As n2 or y increases, the closed-loop poIes migrate towards the open-loop pole positions. Therefore, y-weighting is able to provide a stabilizing effect even if there is model-plant mismatch. The stabilizing ability of yweighting is formally stated in the following lemma. Lemma 2: For any open-loop stable process under long-range predictive control with a steady-state weighting term, y, ana’ a stable estimated model f&(q-‘). there exists a finite ualue, ym, such that y > yrn > 0 guarantees a stable closed-loop system. Proof: Although the following proof, for simplicity, is for the ease of nu = 1, similar proofs can be deduced for any size of the control horizon. Recall that the closed-loop characteristic equation from eq. (42) is Ws = RAA+
with weighted one-step-
q-l&.
Substituting eqs (39) and (41) into the closed-loop equation gives
Model:
Hjaj + H,tx,
Ay(t) = Eu(t - k) + e(t)
_
AA
Control law: 32 = [y(t + k) - w(t + k)]’ + [g,u(t)]‘.
\
j=nl
/
In the absence of the integrated mode, the modeling of the noise is assumed to be always stationary. The closed-loop transfer function is y(t)
=
(48)
gk-lBq-%r) + (gk- iJ-hB + g:)<(t) gk- 18 +
g,2A
The characteristic polynomial does not have an extra order. However, this approach is well known for having an oflet problem at steady state.
Then applying eqs (7) and (12) to eq. (48) turns eq. (48) into WS=~A+/iAq-‘[j~,($-Gj)qhj
Case 3: ARIMAX model with one-step-ahead control and steady-state error weighting:
+ (2)asli-ca]. As y increases to infinity, aj diminishes to zero and a, approaches I/g,. Then eq. (49) can be rearranged to the following format:
Model: Ay(t) = J3u(t - k) + 3 A
where
Control law: 33 = CY(~ + k) - w(t + WI2 + Y(Y.., - w.)’ Transfer function: y(t)
=
(gk- 1
+
YB.)4-1B”‘(t)+ (gk-lEkB + w,q-lHs + g:)<(t) a- 18 + w:A
The introduction of weighting on the steady-state error instead of on the control increments eliminates both the extra order in case 1 and the offset problem in case 2.
+ _ .)gj-1 = a,q-’
+&q-2
+ gsy + _. .
1
1
1295
Long-range predictivecontrol with a terminalmatchingcondition The coefficients Si diminish as y increases. According to Liu and Gertler (1987), there is a limiting value 6, > 0 so that the closed loop is guaranteed to be stable for 1Sil -c 6,. Therefore, a limiting value Y,,, 0 exists such that eq. (48) is stable (Kwok, 1992).
puted as follows:
Y 5. PERFORMANCE
TUNING
OF y-WEIGHTING
The important tuning parameters for GPC steady-state error weighting are listed below: l l l l
with
5.1. Strategy # 1: deadbeat approach Clarke and Mohtadi have shown that the equivalent of a stable deadbeat controller for an observable and controllable linear system is obtained by the following selections: n23
2N-
1, A=0
(50)
where N is the order of a state-space model (Mohtadi, 1987). When an input/output model is considered, the configuration of a deadbeat controller requires the following parameter settings: nu = na + 1,
nl = nb + 1,
y tr(G:G,)
= tr(mGrG)
g,2nu(nu+ 1) 2
= tr(mGTG)
2m tr(GTG) ’ =(nu
+ 1)nugf’
(52)
are roughly the sum of deadbeat output response
Since GPC without y-weighting can be reduced to many of the well-known control methods by different selections of the first three tuning parameters, two new tuning strategies are recommended by incorporating y as an “active” tuning parameter into two common control methods: deadbeat and mean-level control. The reason why it is called “active” is that, while keeping the other tuning parameters constant, only y is adjusted to obtain the desired response. Depending on the performance required, the first of the following approaches can be used for underdamped processes to achieve fast response, whereas the second can be used for overdamped systems to increase robustness.
nl = N,
= tr(mGrG)
Using the “formula” in eq. (52), the control increments
minimum output horizon nl maximum output horizon n2 control horizon nu steady-state error weighting y.
nu = N,
tr(GTrG,)
n2 = nu + nl -
1. (51)
The settings specified above are, therefore, fixed at the commissioning stage. Then y is adjusted to tailor the speed of response. In principle, the value of y can vary from zero to infinity. A very large value of y only makes the controller behave more like a mean-level controller. In order to find a y-weighting relative to the deadbeat control action, one can make the magnitude of GFIG. of the control law in eq. (33) “approximately equal” to mGTG, where m is the degree of detuning from deadbeat to mean level. Then the new control action is roughly l/(m + 1) times the original deadbeat action. If the magnitude of the matrix inverse is approximated by the trace of GrG as in McIntosh’s tuning strategy (McIntosh et al., 1991), y can be com-
weighted by l/(m + 1) plus the mean-level response weighted by m/(m + 1). The initial value of y can be chosen by setting m = 1 so that the control performance is about the average of the two responses for commissioning. Because the formulation of y-weighting already has g 3 as a factor in the inverse portion of eq. (33) the closed-loop performance is independent of any gain changes. This is an advantage of y over the use of L-weighting, which often has to be adjusted relative to gz by scaling. Another advantage is that the slowest control response resulting from a large value of y is limited to that of mean-level control; whereas increasing k can indefinitely decrease the response. 5.2. Strategy #2: mean-level approach For typical industrial plant models, a large value of n2 and a value of nu of 1 are usually sufficient to produce satisfactory control performance for openloop stable plants. The control equation in eq. (33) with nu = 1 and without any weighting can be written simply as
It is obvious from the equation that nonminimum phase model effect should be avoided by choosing a large n2 such that C;znl gi_i is strictly positive. Although a default output horizon window of n2 - nl + 1 = 10 has been found to produce robust response, a large n2 corresponding to a 90% rise time of the process is more desirable especially when commissioning a new controller. Depending on the sampling time, of n2 capturing the rise time can vary from 20 to 50, adding a lot of computational load to the controller. When y-weighting is used in conjunction with n2, it becomes the choice for performance tuning. Then the value of n2 can be reduced considerably to capture only the high-frequency dynamics which usually appear during the first “time constant” period. This is similar to the guideline suggested by Maurath et al. (1985) recommending n2 equivalent to the number of sampling periods required for the process open-loop step response to reach 50% of its final value. This is achieved by making n2 approximately equal to the number of sampling units for one overall process time constant. Therefore, the configuation for commission-
1296
KUN-YU
KWOK and SIRISH L. SHAH
ing a process is proposed to be nl = k,
n2 = 5 + k,
nu = 1,
/
=a
(54) where k is the process time delay T is the dominant or overall time constant process l T=is the sampling time 6 Elf., gi-1 + y > 0 ensures that the sign controller is compatible with the sign steady-state gain (to avoid nonminimum effect). l
.
of the
of the of the phase
For example, if the sampling period is chosen by taking one-fifth of a time constant, n2 will be 5 + k, where k is the units of time delay. Starting with a large value of y. the process response can be adjusted by gradually decreasing y until satisfactory results are observed. The control law for this approach is given below:
Fig. 2.
order over- or underdamped processes are discretized by the same factor of the individual time constants, i.e. zl/rl, = ~2/72,, the resulting step response models are exactly equal. Therefore, an approximate value of y and ns which will provide an effect similar to a larger n2 vlaue (with y = 0), can be computed from eq. (57) when nl and T/Z~are given. Figure 2 shows the trajectories of y corresponding to different n2 and ns. For this particular plot, nl is one and the sampling time is taken as one-tenth of the time constant, which is considered to be an appropriate discretization interval for most industrial processes.
Au = elf.,
si-
*Cw(t+ 4 -f( t + 41 + BJCM) -fMl
where ns denotes the reduced output prediction horizon in order to distinguish it from n2 in eq. (53). One can obtain a feel for the range of admissible y values by equating eqs (53) and (55) and finding the values of y corresponding to different n2s. Without any loss of generality, the setpoints in both equations are assumed zero. After some rearrangement, the resulting relationship between y and n2 is given as
Cz,,
!J-lfCt
Y=
i)C~lm,+l ST-1-C;1,+1
+
gi-lf(f
9j=
S*Sj
where the prime denotes a quantity corresponding to the case where the gain is unity. With this substitution, all gain terms in eq. (56) will be cancelled out. What remains is a relation between y and a combination of the model response coefficients with a standardized unity gain:
y, =
+
OCE,, CT?-1
kZC;2nlgi--lf(t+i)-kkS,C;2n,g:-1
It should be noted that the value of y in eq. (56) is dependent neither on the process gain nor on the process time constant (71, but only on the choice of nl, n2, ns, and the ratio r/~,. In order to show that y is independent of the process gain, the gain term is factored out from all step response coefficients, i.e.
CL
g:-lf(t
+
i,~;‘,+, C;f,,
Trajectories of y for different values of a2 for the special case of au = 1.
.
(56)
4. ADAPTIVE
@f-d2 - c12w+I g:- If’@ + Ok;“=., Cd- 3
gf-If’@
+ i) -_KClf.,
It should be noted again that the choice of y is dependent on the discretization factor, 7/7,, but not on 5. The reason is that when two equal-gain, equal-
GPC WITH y-WEIGHTING
One of the areas where closed-loop control has currently received much attention is in the control of biomedical systems (Linkens, 1984; Vozeh and Steimer, 1985; Linkens and Hacisalihzade, 1990). In particular, continuous infusion of a vasodilator for induced hypotension is often required for cardiac patients during and after surgery. The patient’s response to drug intervention is analogous to a complex and nonlinear chemical process subject to unmeasurable disturbance. The infusion of a vasodilator by a com-
(sI_ 1)’
.
(57)
Long-range
predictive control with a terminal matching condition
puter-controlled drug delivery pump (acting as the final control element) is the manipulative variable. Maintaining the mean arterial blood pressure (MAP) (which is the controlled variable) on the setpoint by the drug infusion is the control objective. Adaptive GPC has heen previously applied in this area by Kwok et al. (1991). The present investigation differs from the earlier work because it involves the use of GPC with steady-state error weighting and a control-relevant, long-range identification strategy for the control of MAP on a physiological system. This identification strategy was developed by Shook et al. (1991), designed to provide an optimal trajectory of output predictions instead of a single optimal prediction. Therefore, the identification provides the best match for a trajectory of predictions, which, in turn, is utilized by GfC with y-weighting. The initial tuning of the controller was done according to tuning strategy #2 in the previous section. Details of the theory, system description, and software package are available in Kwok et al. (1991), Shook et al. (1991), and Kwok (1992). 6.1. Experiment This study was approved by the Health Sciences Animal Care Committee of the University of Alberta. Experiments were performed, in an ethics-approved manner, on six healthy anesthetized mongrel dogs of either gender weighing between 16 and 18 kg. The infusion rate of sodium nitroprusside (SNP), which is a potent hypotensive agent, was used as the manipulative variable. Each run began with a controller initialization sequence which consisted of the infusion of SNP (7-15pgkg-‘min-‘) that produced a small hypotensive response. This also constituted the initial probing sequence for the on-line estimation scheme. Thereafter, the performance of the controller
was evaluated for setpoint changes by requesting a computer-controlled setpoint change of - 30 mmHg from the original baseline. Controller performance was also tested during a series of unpredictable disturbances produced by the administration of norepinephrine (NE) infusion to simulate a hypertensive disturbance and another SNP infusion to simulate a hypotensive disturbance. 6.2. Results For brevity, only one series of runs is presented. Detailed results are found in Kwok et al. (1992). For this particular run, the concentration of SNP was 200 pgml-‘. Figure 3 shows the initial probing effect and the subsequent servo response. Since the subject was highly sensitive to SNP infusion, the initial probing of 40 ml h- 1 already caused a large drop in MAP. During computer control, the MAP controller performed favorably because MAP was well within the acceptable boundary of f 5 mmHg around the setpoint. Two maintained disturbances were made consecutively, one in the absence and the other in the presence of closed-loop computer control. Deviations from baseline MAP under these two conditions were compared. In order to ensure that the disturbances in MAP were made from equivalent baseline pressures, MAP was lowered to similar levels by either a constant infusion of SNP (absence of closed-loop control) or by the controller (presence of closed-loop control). Figures 4 and 5 show the MAP trajectories in response to hypertensive and hypotensive disturbances, respectively. These results clearly quantify the performance of our closed-loop controller based on GPC with steady-state error weighting and the long-range identification strategy.
I
100’
a
011 0
1297
2
4
6
8
10
12
14
2
4
6
8
10
12
14
Tie. Fig. 3. MAP
minutes
response to the initial probing signal and setpoint change
I
KUN-YU KWOK and SIRISH L. SHA%I
1298
Without
closed-loop
control
+-
1-+ With closed-loop
control
601 0
10
20
30
40
50
60
I 70
405 0
10
20
30
‘IO
50
60
70
Time. minutes Fig. 4. MAP response
Without
to norepinephrine
closed-loop
control
+
causinghigh blood pressure.
(NE) infusion
1+ With closed-loop
control
90 85 _________________________________-________________________________ ____------_-
6.5 0
5
10
15
“0
5
IO
15
20
25
30
35
40
20
25
30
35
40
Time, minutes Fig. 5. MAP response
to another
nitroprusside
(SNP) infusion
7. CONCLUSIONS
A terminal matching condition is introduced into LRPC as a weighting on the square of the error between the predicted process steady-state value and the setpoint. The final control law only minimizes the square of the prediction error over a small future prediction horizon and the weighted terminal condition. A computational scheme is also derived for both parametric and convolution models so that, for the
causing
further blood pressure
drop.
ARIMAX model, the steady-state value is calculated without iterating the Diophantine identities and, for the convolution model, a complete set of step response coefficients is no longer required. The latter model type only relies on a few dynamic response coefficients plus the steady-state gain to make future predictions. The closed-loop analysis also shows that the weighted square of steady-state error in the control law not only provides a stabilizing effect as a large
Long-range predictivecontrol with a terminalmatchingcondition prediction horizon but also eliminates the extra closed-loop model order incurred as a result of the integral model. Two tuning strategies based on deadbeat and mean-level concept are provided for using the terminal matching condition with generalized predictive control, which has been successfully applied to the control of mean arterial pressure.
NOTATION
a,
Atq-‘1 bf
ml--‘) ctq-‘1 et Ej(q- '1 ftt+ i) f f,
Fj(q- ‘1 Si
Gj(q- ‘1 G,tq G
‘1
G hi h R;;q-
‘)
J k nl n2 na ns nb nt nu P --I 4 R(q- ‘)
:(4-l) T(q_‘) T u (0 U
approximation error in convolution model reference value in convolution model process model denominator polynomial ith coefficient of the numerator polynomial in G, process model numerator polynomial noise model polynomial ith coefficient of E,(q-I) quotient polynomial from jth-step-ahead Diophantine identity TJA ith open-loop process output prediction assuming no future control action vector off@ + i), where i is from nl to n2 vector off(s) remainder polynomial from jth stepahead Diophantine identity T/A ith coefficient of G,(q- ‘) quotient polynomial from jth-step-ahead Diophantine identity E,B/T transfer function of any process matrix of step response elements lower triangular matrix containing only &s ith impulse response coefficient ith coefficient of H,(q- ‘) remainder polynomial from jth-stepahead Diophantine identity EjBfA performance index time delay in the number of sampling intervals minimum output prediction horizon maximum output prediction horizon order of A(q- ‘) maximum output prediction horizon less than n2 order of B(q ‘) order of T(q- ‘) control horizon first-order coefficient in the denominator of G, backward shift operator polynomial in the linear form of GPC control law value at steady state polynomial in the linear form of GPC control law observer design polynomial model horizon for convolution model control input vector of incremental control inputs
v(q-‘) w(t) W x(r) v(t) Y
1299
polynomial in the linear form of GPC control law setpoint vector of setpoint trajectory noise input process output vector of output prediction values
Greek symbols ith element of the first row of ai [GTG + A + GfrG,] - ‘GT sum of pi from nl to n2 ith element of the first row of ; [GTG + A + G,TrG,] - ‘G:‘lcontrol steady-state weighting Y a critical steady-state weighting which Ya makes closed-loop control stable Idiagonal matrix containing control steady-state weightings 6 degree of the polynomial 4 A
1 A 6(t) Z r.
fth coefficient of L@ differencing operator
(= 1 - q-l) control weighting diagonal matrix containing control weightings uncorrelated random sequence of zero mean overall time constant sampling time
Superscripts
T
estimated value value corresponding to unity gain transposition
Subscripts /”
MV
S Others v a
d-step time delay signal filtered by l/T(q-‘) minimum variance steady state
characteristic equation a portion of %?s REFERENCES
K. J. and Wittenmark, B., 1973, On self tuning regulators.Automatica9, 185-199. Auslander,D. M., Takahashi,Y. and Tomizuka. M.. 1978. Direct diaital ntocesscontrol:nracticeand algorithmsfor micropro&& application.P&. IEEE 66, r99-208. Chesna, S. A. and Ydstie, B. E., 1986, Self tuning and adaptive control of chemical processes, in Adoptive and Learning Systems: Theory and Applications (Edited by Astrcm,
K. S. Narendra), pp. 119-130. Plenum Press, New York. Chesna. S. A. and Ydstie. B. E.. 1988. Controlline nonminimum chase svstems with nredictive controllers. in Workshop ‘on Ado&e Control’ Strategies for Industrial Use (Edited by S. L. Shah and G. A. Dumont), __ PD. 204-231. Lodge Kananaskis, Alberta. Clarke, D. W. and Gawthrop. P. J., 1975, Self-tuning controller. IEE Pmt. D 122.929-934. Clarke, D. W. and Mohtadi, C., 1989, Properties of generalized predictive control. Automatica 25, 859-875.
1300
KUN-YU KWOK and SIRISH L. SHAH
Clarke, D. W., Mohtadi, C. and Tuffs, P. S.. 1987. Gentralised predictive control-Part I. The basic algorithm. Automatica 23, 137-148. Cutler, C. R. and Ramaker, B. L., 1980, Dynamic matrix control-a computer control algorithm, in Proc. JACC, San Francisco. Francis. B. A. and Wonham. W. M.. 1976. The internal model principle of control thkory. Au&n&a 12457467. Garcia, C. E. and Morari, M., 1982, Internal model control. 1. A unifying review and some new results. Ind. Engng Chem Process Des. Den. 21, 308-323. Garcia, C. E., Prett, D. M. and Morari, M., 1989, Model medictive control: theorv and oracti-a survey. Auton&a 25, 335-348. Kwok, K. Y., 1992, Long-range adaptive predictive control: extensions and applications. Ph.D. thdsis, University of Alberta. Kwok, K. Y., Mutha, R K., Shah, S. L., Clanachan, A. S. and Finegan, B. A., 1991, Constrained long-range adaptive control of arterial blood pressure.. J. predictive Adaptive Conrrol Signal Processing 5, 363-374. Kwok; K. Y., Shah, s. L., Clanachan, A. S. and Finegan, B. A., 1992, Evaluation of a long-range adaptive predictive controller for computerized drug delivery systems, in Proceedings o/the IFAC International Symposium on Adaptive ControI and Signal Processing, Grenoble, France. Lambert, E. P., 1987, Process control applications of long range prediction. D.Phil. thesis, University of Oxford. Li, S., Lim, K. Y. and Fisher, D. G., 1989, A state space formulation for model predictive control. A.I.Ck.E. J. 35, 241-249. Lim, K. Y., 1988, Multivariable optimal constrained control aleorithm IMOCCAI. M.Sc. thesis. University of Alberta. Linkins, D. k., 1984, Computer codtrol in biimedicine, in Real-Time Computer ControL (Edited by S. Bennett and D. A. Linkens),Chap. 14. P. geregrin&London. Linkens, D. A. and Hacisalihzadc, S. S., 1990, Computer control systems and pharmalogical drug administration: a survey. .I. Med. Engng Technol. M(2), 41-54. Liu, K. and Gertler, J., 1987, On-line stabilization of adaptive controllers by detuning in a supervisory framework, in Proceedings of the American Control Confwence, Minneapolis, MN, pp. 194-200.
McIntosh, A. R., 1988, Performance and tuning of adaptive generalized predictive control. M.Sc. thesis, University of Alberta. McIntosh, A. R., Shah, S. L. and Fisher, D, G., 1991, Analysis and tuning of adaptive generalized predictive control. Can. J. them Engng 69,97-l 10. Maurath, P. R., Mellichamp, D. A. and Seborg, D. E., 1985, Predictive controller design for SISO systems, in Proceedings of the American Control Conference, Paper FP4, pp. 15461552. Mayne, D. Q. and Michalska. H., 1991, Model predictive control of nonlinear systems, in Proceedings ofthe American Control Conference, FA8, pp. 2343-2348. Mohtadi, C., 1987, Advanced self-tuning algorithms. D.Phil. thesis, University of Oxford. PaDadoulis. A. V.. Tsilixiannis. C. A. and Svoronos. S. A.. i987, A cautious self-&ning’controller for chemical pro: cesses. A.1.Ch.E. J. 33, 40149. Rawlings, J. B. and Muske, K. R., 1992, The stability of constrained receding horizon control. IEEE Trans. Automat. Control (submitted). Richalet, J., Rault, A., Testud, J. L. and Papon, J., 1978, Model predictive heuristic control: applications to industrial processes. Automatica 14, 413428. Shook, D. S., Mohtadi, C. and Shah, S. L., 1991, Identification for long range predictive control. IEE Proc. D 138,75-84. Shook, D. S., Mohtadi, C. and Shah, S. L., 1992, A controlrelevant identification strategy for GPC. ZEEE Trans. Automat. Control 37, 975-980. Sripada, N. R. and Fisher, D. G., 1985, Multivariable optimal constrained control algorithm (MOCCA): Part 1. Formulation and aDDlication. in Proceedinas of the Internatonal Conferench bn Industrial Process VMideling and Control, Vol. 1, Hangzhou, China. Takahashi, Y., Tomizuka, M. and Auslander, D. M., 1975, Simple discrete control of industrial processes (finite time setting control algorithm for single-loop digital controller). Trans. ASME, pp. 354-361. Vozeh, S. and Steimer, J. L., 1985, Feedback control methods for drug dosage optimisation. Clin. Pkarmacokin. 10, 457476.