Adaptive Predictive Control of Armax Plants with Unknown Deadtime

Copyright © IFAC Adaptive Systems in Control and Signal Processing. Grenoble. France. 1992

ADAPTIVE PREDICTIVE CONTROL OF ARMAX PLANTS WITH UNKNOWN DEADTIME L. Chisel and E. Mosca Dipartimento di Sistemi e lnformatica. Unillersitd di FiTenze. Via S. Marta 3. 50139 Firenze.ltaly

Abstract. Consideration is given to i/o versions of receding horizon controllers based on the minimization of multistep quadratic costs with the constraint that the terminal state goes to zero. Plant ARMAX/ ARIMAX models are taken into account. Particular attention is directed to develop a computationally efficient and numerically robust algorithm to be used for adaptive control of plants with unknown, possibly time-varying, deadtime. Keywords. Predictive control; stochastic control; adaptive control; computational methods; systems with delays.

1

2

Introduction

Background

Throughout the paper: prime denotes transpose; 11

Dynamic receding horizon controllers (RHC) with a guaranteed stabilizing property have been recently reported in [1, 2] relatively to SISO ARX/ ARIX plant models. The main result presented in [1, 2] ensures stability of the closed-loop system provided that the control horizon is large enough. This result is proved by extending the classical stabilizing properties of RHC [3] to non-minimal state-space representations.

v Ilb~ v'Qv; and

wt, ~

[w'(j) ... w'(k))'.

Given the discrete-time SISO plant

x(t + 1) yet)

=

=

~

x(t) H x(t)

+

G u(t)

(1)

with x(t) E ~n, consider the following problem. Zero Terminal State Regulation (ZTSR): Find, whenever it exists, an optimal open-loop input sequence 4-1 to the plant (1) initialized from x(O)

The present paper extends the above mentioned stabilizing ilo RHC (SIORHC) approach to ARk = 0,1, ... , T - 1 (2) u(k) = Fk x(O), MAX/ ARIMAX plants with deadtime. It is shown that closed-loop stability can be achieved under the minimizing the quadratic cost (Qy ~ 0 and Qu > 0) same conditions provided that the i/o variables filtered by the inverse of the noise-polynomial, instead J(UL1' x(O)) ~II y} + 11 U~_l (3) of the i/o variables, are constrained to zero. Special attention is paid to the fast and numerically reliable under the zero terminal state constraint computation of the SIORHC control law. (4) x(T) = On The paper is organized as follows. In Sect. 2, some background material on stabilizing receding horizon The RH regulation law is then given by regulation is put forward. In Sect. 3 the desired (5) u(t) = Fo x(t) o dynamic RH regulation problem is stated, its stabilization properties are analysed and a closed-form The following result concerns the stabilizing propsolution is found assuming, for the plant, an ARerties of the RH law based on ZTSR. MAX representation with known dead time. Then, the extension to the servo problem is considered. In Theorem 1 : Let the plant (1) be completely conSect. 4, an efficient algorithm to be used in prac- trollable and detectable and such that its minimal tice for the computation of the control law is given. state-space respresentations have non-singular state In Sect. 5, some simulation results are reported to transition matrix. Then, the RH law (5) relademonstrate the effectiveness of adaptive SIORHC tive to ZTSR exists unique and stabilizes (1) whento control plants with deadtime. Sect. 6 concludes ever the prediction horizon T satisfies the inequality T ~ n dim ell. 0 the paper.

IIby

=

221

lib.

Note that the above result, whose proof is not reported here due to space limitations, generalizes previous results reported in [3] and in [1] relatively to minimal and, respectively, completely reachable and reconstructible state-space representations.

ARMAX plants

3

Consider the stochastic SISO plant described by the ARMAX model

A(d) y(t)

=

B(d) u(t - r)

+ C(d) e(t)

(6)

where: y(t) is the output; u(t) is the input; {e(t)} is a zero-mean white noise sequence; r ~ 0 is the deadtime; d is the unit delay operator, viz. dy(t) ~ y(t - 1); and A, B, C are the following polynomials in d

A(d) = 1+ L~==l akdk, B(d) = L~~l bkdk , C(d) = 1 + L~~l ckdk with zeros in Idl > 1 (7)

Defining, for any scalar sequence w(t), w(t) ~ w(t+

r) and we(t) ~ C-1(d)w(t), the following ARMAX and, respectively, ARX representations for the plant (6) are obtained

- G [O~ Cn-l ... C2 Cl] [Cn-l ... Cl 1 O~_l]

e

He

(8) (9)

It is also convenient to introduce a non-minimal state-space representation of (9) as follows. Define the state vector

(14)

Fact 1: Provided that A and Bare coprime polynomials, ~e = {e,G,Hc} in (12),(14) is completely controllable and detectable, the unobservable modes being located either at the zeros of C or in the origin of the complex plane. Further, provided that 1'1 = max{ n a , ne} ~ nb and en - an # 0, the minimal state-space representations of ~e have a non-singular state transition matrix. [J Hence, by Theorem 1 and Fact 1, we can directly construct a stabilizing I/O RH regulator for the plant (6) by applying the stochastic counterpart 6. =

of ZTSR (cf. Sect. 2) to (13) - (14). Let N T - n + 1 and consider the following stochastic regulation problem.

Stabilizing Input/Output RH Regulation (SIORHR): Given the initial state se(O) of the plant (11), find an optimal open-loop input sequence u~_l

u(k)

= F\sc(O),

k

= 0,1, ... ,T-1

(15)

minimizing the conditional expectation E [J(UL1'Sc(0»

B(d) u(t) + C(d) e(t) B(d) uc(t) + e(t)

A(d) y(t) A(d) yc(t)

where

I se(O)] ,

J(u~_l,se(O»~" yj,~~ II~y Qy

~

0 and Qu

>0

+" u~_l ,,~u' (16)

subject to the constraint

E [se(T) I se(O)]

= 02n-l

(17)

Then, the SIORHR law is

u(t)

(10)

= Fa

sc(t)

[J

(18)

Then,

sc(t + 1) yc(t)

sc(t) + Guc(t) H se(t)

By virtue of Theorem 1, a sufficient condition for the above SIORHR law (18) to be stabilizing is given by the following theorem.

+ Ee(t + 1) (11)

where:

~=[

-an··· -

al

O(n-l)x(n-l) bn ... b2

0(n-2)xn

On-2 I n - 2

0'n

O~_l

G = [O~_l b1 0~_2 1]', H = IO~-l 1 O~_l]

E

1

= [O~_l

A2. A3. 1 o~-tl'

,

then the SIORHR law (18) stabilizes the plant (6),

It is immediate to get from (11) - (12) a nonminimal representation of (8) using the state vector defined in (10). In fact, it can be checked that

sc(t+1) y(t)

ese(t) + Gu(t) Hcsc(t)

A(d) and B(d) are copnme, 1'1 ~ max{n a , ne} ~ nb, Cn - an # 0,

AI.

(12)

{

Theorem 2 : Let the following assumptions be ful-

filled

+ Ee(t + 1) (13) 222

whenever N

~

6.

n = max{n a , ne

+ I}

Proof - The result is a direct consequence of Theorem 1 and Fact 1. [J A number of remarks on the above stated problem and its solution are in order. Remark 1. SIORHR involves costing the i/o variables y and u as well as constraining the filtered by C-1(d) - i/o variables Ye and Ue. Hence, to find

Problem. Given sc(O) = Ssc(O), find UC~_1 as a function of sc(O) minimizing

the solution we need to express either the cost (16) in terms of yc and Uc, or alternatively the constraints in terms of y and u. Choosing the former route, we observe that

Cy sc(O) C" sc(O)

+ +

C yc~ C UC~_1

(19) (20)

subject to the linear constraints

where the matrices C", Cy and

C=

o

[~~///l N-l

Using Lagrangian multipliers to solve the constrained optimization (26) - (27), one can derive a closed-form expression for the SIORHR solution. Only the final solution will be reported here, the derivation being omitted due to space limitations. Assumed N ~ nand L full row rank, the unique solution of (26) - (27) turns out to be

(21)

1

are defined in the Appendix.

0

Remark 2. The stochastic SIORHR problem can actually be reduced to a deterministic one. In fact, from (11), the following predictive model for yc~ given uC~_1 and sc(O) can be obtained

yc~ = W UCLl

+

f sc(t)

+ n e~

(28)

(22)

and the associated SIORHR law is therefore

where:

W f

= £ ([WI

= [' H'

n = £ ([WI

... WT

Uc(t)

]') = [WlL

... «!)T)' H']

OCN-I)Xcn-I)] W

t-

... WT]')

(23) Here: £(v) and U(v) are used to denote the lower and, respectively, upper triangular Toeplitz matrix having v as its first, respectively, last column. Further: Wk ~ Hk-IG and Wk ~ Hk-l E denote the samples of the impulse responses associated with the transfer functions dTB(d)/A(d) and, respectively, 1/A(d) . Since e~ is zero-mean and uncorrelated with sc(O), it can easily be seen that the stochastic SIORHR problem can be solved as if e(t) == 0, which amounts to consider the deterministic counterpart of SIORHR. 0 Remark 3. Note that the SIORHR law (18) is anticipative as sc(t) contains the variables yc;-T+l = YC~t; which are not available at time t. However, sc(t) can be expressed in terms of yt, u t - l and e~t~

-Mc

= Ssc(t) + Ze;t;

(24)

where

S

[T T-IG ... G G]

Z

[T-1E ... E E]

sc(t)

[

sC(tt-_TT) UCt_l

1

[YC~-n+l

t-n-T+l UCt_1

1(25)

From (24) the causal SIORHR law turns out to be

U(t)

= Fa

sc(t)

with

Fa

= FoS

6.

+

-

~(d)A(d)

6.

C(d) e(t)

-

(31) 6. =

= (1 - d)A(d) and c5u(t) ~(d)u(t) u(t) - u(t - 1). Further, assume that at time t the reference sequence {r( k)} is available up to time t + T + To Then, defining the tracking error E(t) ~ y(t) - r(t) and replacing y with E, and u with c5u in the SIORHR problem formulation, one gets the corresponding SIORH servo (SIORHS). As can be verified by repeating the same kind of calculations used to derive (28) - (30), the SIORHS solution is

=

o UCN-l

0

By the above observations, SIORHR amounts to solving the following quadratic optimization problem with linear equality constraints.

l

A(d) y(t) = B(d) c5u(t - T) with A(d) =

sc(t)

(29)

sc(t)

[QC f 2 + (IN - QcLMcl)Nc] S [1 0 ... 0]' Fb Mc [Cl C2 ]' Q" [Cl C2 ] + W{C~ QyC I W I Qc L' (LMcL,)-l Nc [Cl C2 ]' Q"C" + W{C~ Qy(Cy + Clft); (30) Cl, C2 , C", Cy are defined in the Appendix; and fl, f2' W l , W 2 , Land S are defined in (23) and (25). The formulae (28) - (30) coincide with the ones reported in [1] whenever C(d) = 1, Q" = if!yIr and Qy = if!"Ir· Remark 4. (Extension to the servo problem). The problem (15) - (18) can be easily extended to cover the servo case. In this connection, assume that the plant is modeled by the ARIMAX equation Fb ib

l

VIa

= ib

where:

2

= [~~]

(27)

= F b Sc (0) + F f

and the corresponding SIORHS law is

uc(t) = ib sc(t) 223

+ if

T+2-n rCT+ T

(32)

2-degrees-of-freedom rc:t;"t~-n

(33)

in which the feedback gains Fb and fb are the same as computed for SIORHR, i.e. by (30), while the feedforward gains are given by

Mc 1Qc]

F,

[Onx(T-l)

f,

W{C{QyCr [1 0 ... 0)' F,

+ (IN -

QcLMc1)

(34) where the matrix Cr is defined in the Appendix. Eq. (30) and (34) are closed-form expressions for computing SIORH control (SIORHC). Another algorithmic solution computationally and numerically more advantageous will be described in the next section.

in which: the matrix U21 is (2N -1) x (2N -1); 82 ~ [eT-l J-lT-1 ... e1 J-l1 eo]'; and 81 ~ OOT to denote that the entries of Vl1Xl + V12X2 are constrained to zero, i.e. are costed with an infinite weight.

A computation ally efficient method to solve ECLS has been proposed in [5] and therein used to solve quite a general class of optimal quadratic filtering and control problems. This technique is based on the triangularization

~12l

(40)

U22

in which U2l and U21 are (2N - 1) x (2N - 1) unit upper triangular matrices. The triangularization 4 SIORHC algorithm (40) must be obtained by Gauss and fast Givens elementary transformations on the rows of M suitIn (32), (30) and (34) a closed-form solu- ably chosen in such a way that the matrix in the tion to SIORHC is given for the general AR- RHS of (40) is an equivalent parametrization of the MAX/ ARIMAX case. The resulting calculation re- original ECLS problem (35). The details of such a quires O(N3) flops irrespectively of C(d). In [4], a triangularization procedure are deferred to [5]. Here fast algorithm requiring O(N) flops has been pre- it suffices to point out that o~ce (40) has been persented relatively to the ARX/ ARIX case. formed, the bottom row of U22 directly yields the Hereafter the algorithm in [4] will be extended coefficients of the SIORHC law, viz. to ARMAX/ ARIMAX plants with deadtime. SIORHS, to which we refer for greater generality (41) u(t) = -[0 0 ... (SIORHR is actually encompassed by SIORHS), can be recast into the following Equality Constrained Least Squares (ECLS) problem: Remark 5. Deadtime - For r > 0, the control law min 11 V 21 X l + V 22 X 2 II~ under V l1 X l + V12X2 = 0 (41) is anticipative due to the presence of yc(t + x, ' (35) 1), ... , yc(t + r). However, it can be made nonwhere Xl is the vector of the filtered i/o samples anticipative by exploiting costed in (26), viz. n

z(k)

~

[

Yc(t+k)

8uc(k) ] yc(k + r)

= 2: b; uc(Hk-i)-a; yc(Hk-i)

(42)

;=1

(36) The relation (42), subsequently applied for k

while X2 is the vector of the available reference and filtered i/o samples, viz.

Further the matrices Vij, i , j = 1,2, are defined in the Appendix. Finally, assumed that Qy diag{J-ll, J-l2, ... , J-lT } ~ 0 and Qu diag{eo, e1,· .. , eT-d ~ 0, the weight matrix ~2 is defined by

=

~2

=

= diag{eT_1, J-lT-l,""

e1, J-ll, eo}

(38)

Hence the problem (35) is fully parametrized by the matrix

M

=

r, ... , 1, allows to eliminate the dependence of (41) on yc(t + r), ... , yc(t + 1) and thus obtain a nonanticipative control law . 0 Remark 6. Direct computation of the control variable - So far we considered how to compute the control gain. This is actually more favourable than directly computing the control variable whenever the control law design has to be performed occasionally, i.e. not every sample time. If this is the case, the control variable can in fact be evaluated at each sample time by regressing the available control gain with the new control regressor without need of repeating the time-consuming design procedure. On the other hand, the direct computation of the control variable is computationally more parsimonious whenever the control law design has to be performed each sample time. The direct computation can be carried out as follows. The vector

(39) 224

is first computed. Whenever T > 0 the future samples yc(t + k), 1 ::; k ::; T, in X2(t) are evaluated by (42). Then the triangularization (40), in which U21 and U22 are now column vectors, is performed yielding the required value of the control variable as the bottom component of U22 . 0

plant dead time 0". Even if 0" is unknown and/or time-varying a good control performance is, however, expected provided that the identifier and controller regressors are large enough to include all the relevant i/o samples. More precisely, this occurs whenever na 2: n a , nb + f 2: nb + T and f::; T. Remark 7. Computational complexity - In the fol- The experiment in Fig. 1 concerns the case of exlowing considerations on computational complexity act mode ling with 0" = 0.95s. Fig. 2 shows the it will be assumed that N ~ n. The essential part behaviour in case one knows an upper bound for 0". of the proposed algorithm is the triangularization Here we chose f = 0 and nb = 6. As can be seen, the of [d Ut] in (39) . Should the matrix [d U1] be performance is only slightly degraded as compared full, this would take O(2N3) flops. Due to its high to Fig. 1. Larger values of nb yield a behaviour sparsity, however, only O(N 2 n) flops in general, comparable with that of Fig. 2. Finally, Figs. 3 O(2N n 2 ) in the ARX/ ARIX case, are required. To and 4 demonstrate SIORHC ability to follow sudcompute the control gain an extra load of O(N3), den dead time changes. The deadtime 0", initially set O(n 3 ), flops is required for the servo, respectively to 0.2 s, has been changed to 0.4 s at time tl = 60 s, 120 s. Here we chose nb 5 regulation, case. On the other hand, to directly then to 0.6 s at t2 compute the control variable only O(N n) additional and f = o. The transient output behaviour over flops are required. Note that, in any case, the the time interval 0 -:-180 s is shown in Fig . 3, while proposed algorithm is computation ally more con- the steady-state behaviour after the last change is venient than the direct computation of (30), (34) illustrated in Fig. 4. which always involves O(N3) flops . 0

=

Remark 8. Numerical issues - The use of the above algorithmic approach is also beneficial from a numerical point of view , as matrix inversions are avoided. By computer experiments it has , in fact, been demonstrated that the proposed algorithm yields the correct solution (within the computer accuracy) even when A(d) and B(d) have common stable factors , in which case the direct evaluation of (30) fails due to singularity of the matrix LMcL'. Note that the capability of handling common factors in A(d) and B(d) is a necessary requirement since in applications, where disturbances are usually coloured noises , the presence of stable common factors in A(d) and B(d) represents the generic case. 0

5

6

=

Conclusions

A stochastic predictive control approach to the design of SISO dynamic receding horizon controllers has been presented . A fast and numerically reliable UD-factorized SIORHC algorithm has been developed for possible adaptive control applications. Finally, simulation results have been reported to demonstrate the ability of adaptive SIORHC to control plants with time-varying deadtime.

References [1] Mosca E. , J .M. Lemos and J. Zhang (1990). Stabilizing I/O receding horizon control. Proc. 29-th CDC, Honolulu (U .S.A.), 2518-2523 .

Simulation results

SIORHC can be used as an adaptive controller if equipped with a recursive identifier for on-line estimation of the plant parameters (cf. [6]).

[2] Clarke D.W . and R. Scattolini (1991). Constrained receding horizon predictive control. lEE Proc. Pt. D.

The following continuous-time nonminimum-phase plant (ECC-91 benchmark example) has been considered :

[3] Kwon W.H. and A.E . Pearson (1975a) . On the stabilization of a discrete constant linear system. IEEE Trans. Aut. Control, AC-20 , 800801.

G(s) - e-'" 1- s (s + 1) (1 +

115 s )2

[4] Chisci L. , E. Mosca and S. Zengl (1991). A fast algorithm for a stabilizing receding horizon controller . Proc. 30-th CDC, Brighton (U .K.) , 517-518 .

(43)

with a sampling time T. = 0.25s , corresponding to the following discrete-time plant structure: na = 4,nb = 4 and T = LO"/0 .25J .

[5] Chisci 1. and G. Zappa (1991). Systolic architecture for predictive control. To appear in Proc. i-st IFAC Works. on Algorithms and A rchitectures for Real- Time Control, Bangor (U.K.) .

The SIORHC design knobs have been fixed as follows: N = 15, Qu = Qy = I . As for the assumed structural parameters, that are denoted by (n a , nb, ne, f) not to be confused with the true ones, we fixed na = 3 and ne = 0, while nb and f have been changed from experiment to experiment reflecting a different degree of prior knowledge on the

[6] Mosca E . and J . Zhang (1991) . Globally convergent predictive adaptive control. Proc. i-st ECC, Grenoble (France), 2169-2174. 225

'" 9

'" 9

.' A

8

! 1\

I\.

1\

/

8

ci

ci

Fig. 1

-'- -v

'"

"'0 00

-'----v

.l............,

3) . 0

15 . 0

Fig. 2

~5 . 0

~

'" '"

"'0 00

SO . O

9'"

l ..

11"

I

~

IV

3) . 0

15 . 0

~

~5 .

SO . O

h

h

h

0

8

.,

Fig. 3

.....,

..,

If-. ~ ~

If'!;

Fig. 4

If',

u.- -V

~

~

s,

s, ?'l

0 . 9)

• 0 . 00

l. 35

?'l1-4 . 20

1.80

~.

~.

35

SO

~.

65

Appendix - Definitions of Cl , C2 , Cu, Cy , Cr, Vll , Vi2, V21 and V22

C= .c ([1 Cl

...

O~-lJ') = [~~ III]

Cn -l

N-l

C =

[O (n-l)Xn

u

U ([Cl

ONxn Cn Nx(n-l)

Vu

T{ [ ~ 2n-l

_

[

U

d') I '--r ([1

~

2N-l

2n-l

-

2T -l

III {[ ~ 2n-1

2N-l

... Cn -l]')

ONx(n-l)

O(n-l)X(n-l)] ONx(n-l)

0'T-n ]')] n

b

all.. ... .

Onx (T-l )

-b

an

] E IDTx(T+n -l ) ;n

O(N-I)X(T-l)

0(N-l ) xn

V2l

ON

-

([an-l ... all]')

~

l

U ([Cl

[On-l

y

Cl· .. Cn-l

vR)]_ r 1 -

~

(2) V12

C =

ONx(n-l)

Cr -- [ U ([Cl0 ...

I II

... Cn-I]') ] ,

1

v.(l) ] _ 22 2n-l

r1

:

6

0

'"0'

)

V

_ [ ir(l)

12 -

v12

vg)]

'H

~

o o

v.22(2) --

-1

o

o -1

-Cn-l

o -Cn-l

226

o

o

-1

-Cn-l

o

o

E ~Tx(T+3n-2)

+. 80