Optimal input design for an AR model with output constraints

Optimal input design for an AR model with output constraints

0005-1098/84 $3.00+0.00 Pergamon Press Ltd. © 1984 International Federation of Automatic Control Automatica, Vol. 20, No. 3, pp. 359-363, 1984 Printe...

387KB Sizes 0 Downloads 145 Views

0005-1098/84 $3.00+0.00 Pergamon Press Ltd. © 1984 International Federation of Automatic Control

Automatica, Vol. 20, No. 3, pp. 359-363, 1984 Printed in Great Britain

Brief Paper

Optimal Input Design for an AR Model with Output Constraints* T. S. NG,I" Z. H. QURESHII" and Y. C. CHEAH? Key Wards--Identification; optimal systems; parameter estimation; input design; autoregressive model; output power constraints.

However, in many situations, a more appropriate constraint is on the system output fluctuations. For example, in an industrial process where' production cannot be disrupted, the quality of product output has to be regulated within certain prescribed tolerance limits. Alternatively in a biological experiment or a clinical trial, the well being of the subject under treatment may be of prime importance, Input design with constrained output has been considered by a number of authors. S6derstrbm, Ljung and Gustavsson (1974), S6derstrbm, Gustavsson and Ljung (1975) and S6derstrbm (1975) have examined a first-order system in detail and raised a number of interesting conjectures. Zarrop (1979) has shown that with certain system order constraints a D-optimal design could be achieved with finite number of input frequencies without feedback. Ng and Qureshi (1981) have shown that for an autoregressive model, these input frequencies can be obtained by solving a set of nonlinear equations. Ng, Goodwin and Payne (1977) and Ng, Goodwin and S6derstr6m (1977) considered the joint input and feedback design and showed that the minimum variance control law achieved optimal design for a large class of systems including the autoregressive model. The dependence of the optimal experimental conditions on the unknown parameters is a common occurrence in optimal experiment design problems for dynamic and nonlinear systems. Two approaches have been used to overcome this difficulty. The first approach uses the Bayesian formulation by assuming that an a priori probability distribution for the parameter vector is given. The second approach assumes that a nominal value of the parameter is given. In practice, this is usually obtained by some preliminary experiments. The main experiment is often divided into several stages where a design is carried out at each stage using the estimated parameter vector obtained from the previous stage. We have adopted the latter approach in this brief note. In this paper a method to obtain a D-optimal input design for an autoregressive model with output power constraints is given. D-optimal design is defined as the design that maximizes the determinant of the information matrix. We show that a Doptimal experiment consisting of sinusoidal input test signals without feedback can be obtained by solving two sets of p/2linear simultaneous equations, where p is the number of parameters to be estimated, and a simple polynomial equation. Furthermore, it is shown that under mild physical constraints, the achievable accuracy of parameter estimation is the same as in the case of minimum variance control law together with input test signals (Ng, Goodwin and Payne, 1977; Ng, Goodwin and S6derstr/Sm, 1977).

Abstract--The output power constraint problem of optimal input design for parameter estimation for an autoregressive model is considered. A simple method to obtain an optimal design by solving two sets of p/2-1inear simultaneous equations and a polynomial equation of p/2th order is proposed and two nontrivial examples are given to illustrate this methodology. 1. I n t r o d u c t i o n

WITH the current pace of development of science and technology, it is unavoidable that there will be increasing complexity in experimental investigation and, moreover, in the problem of obtaining maximal information from measured data. Experiment design for dynamic system identification has attracted considerable attention because most industrial, biological, economic and sociological systems are dynamic in nature. Different aspects of identification and experiment design problems, e.g. input signal synthesis, sampling instants, feedback settings, etc., have been discussed; see for example, survey papers Astrbm and Eykhoff (1971), Mehra (1974a) and Gustavsson, Ljung and St3derstriSm (1977). In the study of optimal input design problems, most authors imposed constraints on the input sequences (Mehra, 1974a; Gustavsson, Ljung and S6derstr~Sm, 1977). One of the earliest papers on optimal input design was that of Levin (1960). Levin showed that for a moving average model with uncorrelated disturbances and with constrained input power, the control sequence is uncorrelated in time. Litman and Huggins (1963) considered the use of growing exponential control signals. A number of authors also formulated the input design problem in the framework of optimal control theory, e.g. Nahi and Wallis (1969), Aoki and Staley (1970), Mehra (1974b) and LopezToledo and Athans (1975). Even though the optimal control approach leads to useful designs in specific cases, unfortunately it is frequently unsuitable for parameter estimation purposes (Zarrop and Goodwin, 1975). For large experimental times, it is possible to use asymptotic expressions for the information matrix (closely related to the covariance matrix of the parameter estimator) and many results in linear regression problems could be applied m u t a t i s m u t a n d i s to the dynamic system with suitable interpretation of the variables (Viort, 1972). Using this approach, the optimal experiment design problem has been studied extensively by Goodwin and Payne (1977) and Zarrop (1979). Constrained input design has a wide range of applications, and in many cases, leads to analytic solutions of the design problem.

2. P r o b l e m f o r m u l a t i o n

* Received 5 August 1982; revised 8 November 1982; revised 14 September 1982; revised 18 October 1983. The original version of this paper was presented at the 8th IFAC Congress on Control Science and Technology for the Progress of Society which was held in Kyoto, Japan during August 1981. The published proceedings of this IFAC meeting may be ordered from Pergamon Press Ltd, Headington Hill Hall, Oxford OX3 0BW, U.K. This paper was recommended for publication in revised form by associate editor M. Gevers under the direction of editor B. Anderson. t Department of Electrical and Computer Engineering, University of Wollongong, P.O. Box 1144, Wollongong, 2500, Australia.

Consider the autoregressive model Yt = a l Y t - a + ,.. + a . y ~ _ . + b u t - i + ~,

(1)

where {~,} is a white noise sequence having a Gaussian distribution with variance a 2. We assume the system is stable, i.e. the polynomial A ( z - 1 ) = 1 - al z - L - a2z - 2 - . . . - a . z - " has all its roots outside the unit circle. We consider the case where the variance of the system output {y,} is constrained to be less than or equal to some prescribed value (say W), i.e. 1

W>t E{y~}

tTM "

= ~ | d_

359

dFy(og)

(2)

360

Brief Paper

where Fy(o) is the spectral distribution function of the output sequence {y,}. Following Ng, Goodwin and Payne (1977) the average asymptotic information matrix for the system parameters (al . . . . . a,,b) r is given by

t I

"

= EF/ok~2t~= 1

i -

Result 2.1. (Ng, Goodwin and Payne, 1977.) For the system described by (1) subject to P0 ~< W, a necessary and sufficient condition to maximize log det M' is for Pt = p2 = . . . = p , = 0 .

(3)

where 0 = (al . . . . . a., b)T = (~r, b)T, Er/o denotes the conditional expectation over the data given the parameters and ~h is the residual sequence given by ?h = Yt - - a l y t -

(2) such that (15) is maximized. Results established for the problem are briefly summarized below.

... -- a.yt

. -

b u t - I.

(4)

The part of the information matrix for the parameter 0-2 is independent of any input sequence {ut} (Goodwin and Payne, 1977) and hence will not be considered here. Thus

Result 2.2. (Ng, Goodwin and Payne, 1977.) For the system described by (1) subject to P0-%< W, logdet ]~' (15) achieves its maximum by a minimum variance control law together with a white input signal. Result 2.3. (Zarrop, 1979.) For the system described by (1) subject to Po -%<14/,logdet M achieves its maximum by an input test signal without feedback comprising p/2 frequencies for p = even and (p + 1)/2 frequencies for p = odd, where p is the number of parameters to be estimated.

t%h - - =

-Y,-i,

c%f = db

-

Oar

--

i=

1 .....

n

(5).

u , _ 1.

(6)

Substitute (5) and (6) into (3) which yields

...Y,-lY,-,

-y2-1

f

1

= E--~ tr"

ILY,-lu, 1 I

(7)

I

[ Yt - nut - 1

LUt-lYt-l'ut

lYt-n

I at-i

LU,-~J

= ~

/8)

Result 2.4. (Ng and Qureshi, 1981.) For the system described by (1) subject to P0 -%
where the partitions correspond to the partitioning of 0 into

(~r, b)T,

det M*

-

P'~(Po - a2) b20-2~n + 11

We now define {Pl} as follows:

E{yt-jYt-k} = PU-kl.

(9)

Substitute (9) into (8) and using (1), we have Po F =

Pl

...

P.~-t]

Px

Pl Po

Pl

(10)

:

P.-1

..-

J

1

G =_(V-F~) o 1

H = 7~ (P0 - 2aXV+ ctTF~ - 0-2) b"

3. Main result Consider an input test signal given by (11)

(12)

where V is given by VT= [Pl, P2..... P.]"

which is also independent of the parameters a, i = 1, ...,n. However, there is a difference in regard to the minimum value of the constraint I4/. For the minimum variance control case, the minimum W is a 2 while in the open loop input design, the minimum W is (Co + max C~), i = 1. . . . . n which is dependent on the parameters al, i = 1. . . . . n. The definition of C~, i = 0, ..., n is given in Section 3. If W is chosen greater than or equal to 2Co, then the physical constraint (18) will always be satisfied.

(13)

U, = ~ i-I

m[coso)if

(16)

where m~ > 0, i = 1..... l; l = (n + 1)/2 for n = odd, or (n + 2)/2 for n = even, are chosen to satisfy P0 = W. Using (1), (7), (9) and Result 2.1, we arrive at the following equations for optimality

P0= Z gi+ C0= W i

1

From (8) we have t91 = Z g i c o s o ) i i=1

l o g d e t / ~ = logdet (F) + log(H - GTF-1G) + log(K) where K

= b2/0- 2In+

(14) P2 = ~, glcos2o~ + C2 = 0

11. Substitution of (11)-(13) into (14) yields

log det 10" = log det (F) + log (Po - 0-2 _ VTF- 1V - log(K).

+ C1 = 0

i=1

(15) p,= ~ g~cosnog~+C,=0

Our problem now is to find experimental conditions subject to

i=1

(17)

Brief Paper

361

where n=21-1

mi

i -- 1, ..., I

g~ = A(e~,,)A(e_j~,), mi = m[2b 2

(17)

0.2 f* e-Jko Ck = ~ j _ , , A ( e t ~ e _ j , ~ )

Co ~ [CjI, j = l ,

d co,

k = O, 1. . . . . n

2. . . . . n

1 . .1EIlli:o]

I

for n f odd and 2 1 - 2 for n f even

~s0)1" cos 2091"

cos 0)2*... cos 0)

cos(l- 1~1"

... c o s l l -

Yr = ayt- 1 + bu~_ 1 + e~

>I t A(z) A(z-1)J

which is equivalent to CodiCil,

t-I

I

(25) with g~, i = 1..... 1 known, mi or m[ can easily be calculated from (17'). Thus we have shown that a D-optimai input for an AR model can be obtained by solving two sets of/-linear equations together with a lth order polynomial equation.

Physical considerations (SSderstr6m, 1975) also requires

W-

l)to?.]

hi

4. Examples E x a m p l e 4.1. Consider the first-order example

A(z) = 1 - a l z - a2z z - ... - a,z ~.

~.~J

g2

j=

(18)

1..... n.

(26)

where - 1 < a < 1, E~,2 = #" and Ey~ = IV. Detailed analysis of this example has been given in St~derstrbm (1975). Note that optimal design requires at most one input frequency, Following Section 3 we arrive at the following equations 0.2

We now proceed to obtain a solution to (17) as follows: first rewrite (17) into the following form



(27)

gl = W - 7 - - - ~ = V o l--a"

a0. 2

g l c o s k % = hk,

k = 0 ..... n

gl cos co --

(19)

1 - a2

vl

(28)

i=1

where ho = W - Co, h~ = --Ck, k = 1, ...,n. For n ffi odd, we have 21 simultaneous equations and in the case of n = even, we introduce the additional equation ~ g ~ cos (n + 1)0)~ = h.+ 1 = - C,+ 1

(20)

i=1

so that 21 simultaneous equations are still available. Equations (19) and (20) can be transformed (see appendix) into ~ g ~ cos k tot = vk

where ml

gt ----

Let 0)* be the root of the polynomial do + cos 0) = O, then using (27) a n d (28) we obtain

which implies

(21)

131 COS0)* = - d o ~ - - :I::

V0

k = 0 ..... n for n = odd and k = 0 ..... n + 1 for n = even, where vt and hh are related by the recursive formula lnt (1/2) Int (1/2k-J)

~

~

j=l

r=0

det/~* ~ 2W(W-

(22)

and Int ( ' ) denotes the integer part of the n u m b e r enclosed by parentheses. Now let x ~ = c o s 0 ) , i = 1, ...,l be the ith roots of the polynomial P ( x ) = do + d l x + d2x 2 + ... + ~ = 0

0.a)

0.4b2

For comparison, we obtained the optimal frequency by the direct method as follows: first calculate the information matrix using (7) and (26), we have

1[w ( )

(23) = ~

where the coefficients of P(X), i.e. do..... di-1, satisfy

-- £10.2

W(I -- a 2) - 0.2

and

) k . ( - 1 ) ) (,r + j )(2,+2j)tk-2 /2 k - I

(29)

dovo= -vl

i=1

vh = hk--

1 + a 2 -- 2a cos 0)'

c o s co

- -

a

0"2

w -

b c o s 0) - a /

vl ,-~

v2

vl.

dx

v2,-~

ids_ 1

=

-.vl.x

.

a2

~

"1

(24)

l+a

/ - v ~: , - q

Equation (24) is derived from (21) and (23) as follows: multiply d~ to the ith equation of (21), i = 0,1 ..... l - 1, adding the first (l+ 1) equation and using (23), we obtain the first equation of (24). Repeat the above procedure I times with successive ( l + 1) equations of (21) yields (24). We note that the system {1,cos0), cos z to..... cos z0)} forms a Tchebycheff system and it follows that (23) has I distinct real roots in the interval [ - 1,1 ]. Solving for do ..... dr- ~ from (24), the roots of (23) can now be determined and the frequencies 0)~', i = 1..... I correspond to components of the optimal input. The weights of these components can then be obtained from the first l equations of (19), i.e. we solve the following set of linear equations

-2acos0)/

0.2

~l

/w- _--=TJJ 1

and -

1

a2

cos0) - a) 2 (W -

0-2

To obtain the optimal frequency, we differentiate (30) with respect to 0), thus

362

Brief Paper -- a0.2

d --

d09

det -M = 0 ~ cos co*

3.727 -0.5781F,01 = F-1.784 l -0.578 1.784JLd,] L+0.4393

W(1 -- 02) -- 0.2

back substitution into (30) gives detM*

or

04<

2 W ( W - 0.2) b20-,t.

0.096]"

We note that physical constraint (18) requires

Using (23), it follows that the optimal frequencies can be obtained from the solution of

0.2 W ~ > - 1 -lal

- 0 . 4 6 4 + 0.096cos 09 + cos2w = O.

and that the method presented reduces the computational burden substantially even for a first-order example.

which yields cos09' = 0.635, cos09~ = -0.731 or 09* = 0.883 and 09~ = 2.391. Using (25) it follows that

Example 4.2. Consider the following second-order example Yt = alyt-1 + a2yt-2 + bu~-I + t,

with values of al, a2 such that A (z- ~) have roots outside the unit circle, Ee~ = a 2, and Ey~ = W.. Following Section 3, we have gl + g 2 =

W-Co

gl COS3091 + g2 COS3092 =

-- C3

where

mi gi = ~ ,

i=1,2.

(32)

f(09) = (1 + a 2 + a 2) - 2ax(1 - a2)cos09 - 2a2 cos209

C1-

C2-

a2

0.2

1 +a2 1 - ( a 24-2a2) 4-a22 al

trz

l+a21-(at

2+2az)+a

a 2 + a2(1 -- a2)

a2

1 4- a2

1 - (a 2 + 2a2) + a22

(1 4- a2)

0.2 1 - (at2 + 2a2) + a22"

Physical constraints also require

W-Co>~lC~l,

j = 1,2,3.

Assuming al = 0.5, a2 = - 0 . 1 , b = 0.5 we obtain Co = 1.273, Ct = 0.578, C2 = 0.16, Ca = 0.021 and choosing W~> 2Co so that the physical constraints are satisfied, say W = 5, we have the following: gl + g2 = 3.727 gl COS(DI + g2COS092 = --0.578 g~ COS2091 + g2 COS209: = --0.160 gl COS3091 + g2 COS309Z = -- 0.021, Using (21) and (24), we obtain gl + g2 = 3.727 gl costol 4- g2cos092 = - 0 . 5 7 8 gl cos2091 + g2 cos2092 = 1.784 gl cos309~ + g2 cos3092 = - 0 . 4 3 9 and

l_rg,l_F

3.7271

Lg2/- L-0.578_1

or

g'l 1-1"571] 1_2.1563'

5. Conclusions Optimal input design problem for parameter estimation for an autoregressive model with constrained output variances has been considered. A method has been presented to obtain the optimal input by solving two sets of linear equations and a lth order polynomial equation. It does not require optimisation techniques nor the calculation of determinants. Two examples have also been given to illustrate this methodology.

2

al(a 2 + 2a2 - a22) C3 =

l

-0.731/-

Back substitution into (32) we obtain ml = 0.820 or m~ = 0.453; m2 = 4.480 or m[ = 1.058. We remark that this is a nontrivial problem and the method presented gives an optimal input design with very little computational burden. We also note that the calculation of C s, j = 0,..., n + 1 is nontrivial and a numerical method to evaluate Cj, j = 0, .... n + 1 can be found in Jury (1964).

gl cos 209~ + g2 cos 2092 = - C2

1 --

L0.635

g2A =

gt cos091 + g2cos092 = - C 1

Co =

F1

(31)

(33)

ReJerences Aoki, M. and R. M. Staley (1970). O n input signal synthesis in parameter identification. Automatica, 1,431. Astrfim, K. J. and P. Eykhoff (1971 ). A system identification - - a survey. Automatica, 7, 123. Goodwin, G. C. and R. L. Payne (1977). Dynamic System Identification: Experiment Design and Data Analysis. Academic Press, New York. Gustavsson, I., L. Ljung and T. S6derstr6m (1977). Identification of processes in closed loop - - identifiability and accuracy aspects. Automatica, 13, 59. Jury, E. I. (1964). Theory and Application of the Z-TransJbrm Method. John Wiley, New York. Levin, M. J. (1960). Optimal estimation of impulse response in the presence of noise. I R E Trans Circuit Theory, TC-7, 50. Litman, S. and W. H. Huggins (1963). Growing exponentials as a probing signal for system identification. Proc IEEE, 51,917. Lopez-Toledo, A. A. and M. Athans (1975). Optimal policies for identification of stochastic linear systems. IEEE Trans Aut. Control, AC-10, 754. Mehra, R. K. (1974a). Optimal input signals for parameter estimation in dynamic systems - - survey and new results. 1EEE Trans Aut. Control, AC-19, 753. Mehra, R. K. (1974b). Optimal inputs for linear system identification. IEEE Trans Aut. Control, AC-19, 192. Nahi, N. E. and D. E. Wallis Jr (1968). Optimal inputs for parameter estimation in dynamic systems with white observation noise. Paper IV-A5, Proc. JACC, Boulder, Colorado, pp. 506-512. Ng, T. S., G. C. Goodwin and R. U Payne (1977). O n maximal accuracy estimation with output power constraints. IEEE Trans Aut. Control, A C - 2 2 .

Brief Paper Ng, T. S., G. C. Goodwin and T. S6derstr6m (1977). Optimal experiment design for linear systems with input-output constraints. Automatica, 13, 571. Ng, T. S. and Z. H. Qureshi (1981). Optimal experiment design for autoregressive model with output power constraints. IEEE Trans Aut. Control, AC-26. Ortega, J. M. and W. C. Rheinboldt (1970). lterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York. S6derstrbm, T. (1975). Notes on the design of optimal identification experiments. Report UPTEC 7563R, Institute of Technology, Uppsala University, Sweden. S~Sderstrbm, T., I. Gustavsson and L. Ljung (1975). On the accuracy problem of identification. Proceedings of the 6th IFAC Congress, Boston, MA, paper 18-1. S6derstrbm, T., L. Ljung and I. Gustavsson (1974). On the accuracy of identification and the design of identification experiments. Report 7428, Lund Institute of Technology, Division of Automatic Control. Viort, B. (1972). D-optimal design for dynamic models: part I theory; part II applications. Department of Statistics, University of Wisconsin, Madison. TR 314, 316. Zarrop, M. B. (1979). Optimal Experiment Designfor Dynamic System Identification, Springer, Berlin. Zarrop, M. B. and G. C. Goodwin (1975). Comments on optimal inputs for system identification. IEEE Trans Aut. Control,AC20, 299.

363

(19) becomes

i=l

~

j=O

r=0

Vk= ~ gi cosk%

(A3)

Now let

i=0

and interchange the order of summation, we obtain lm (1/2k) Im (1~2k-j)

hA =

~

~

j=O

r=O

" ,+j )(2,+ k 2j)Vk. 2j. (-- 1)J(~

Int (1/2k) lnt

hkvk =

E J=1

(1/2k-J~

E ,~ o

(-ly(,

lnt (1/2k)

",+j (2r+2j)l;k-2j k . (A5)

(~,)

r=O

and finally noting

yields (22). (A 1)

(A4)

Rearrange terms in (A4) gives

~,

(2,+k2j) cosk- 2Jco

E r=O

(A2)

r=O

(-1)J(~ +~)

j=O

( - 1)J(;+J)(2,k+2j)COSk 2Jco,.

Int (l/2k)

Int (l/2k) Int (1/2k - j )

~

i=1

-

Appendix: Derivation of (21) and (22) Using the identity coskco=

Int (I/2k) Int (l/2k - j )

l

hk = 2 g, cosk~, = ~, g, E

(k,) = 2k=1