Optimal Input Design for Dynamic Model Prediction Accuracy

15th IFAC Workshop on Control Applications of Optimization The International Federation of Automatic Control September 13-16, 2012. Rimini, Italy

Optimal Input Design for Dynamic Model Prediction Accuracy Ke Fang ∗ Tom Shenton ∗ ∗

Centre for Engineering Dynamics The University of Liverpool, Liverpool, UK. e-mail: [email protected]

Abstract: A new optimality criterion for the optimal input design for discrete dynamic nonlinear system identification is presented. The criterion weights the parameter covariances by the magnitude of regressors in order to reduce the prediction error. An iterative optimization procedure to the optimal signal design is proposed. The effectiveness of the optimal test signal design is demonstrated by the system identification and validation of a nonlinear multiple-input singleoutput (MISO) automotive engine fuelling model. The output-error simulation accuracy of the resulting model is compared with those of models identified by means of commonly used nonoptimal test signals and optimal signals designed by alternative optimality criteria. The proposed test signal design method is shown to provide superior outcomes in output prediction fit and also to allow the application of input and output constraints as required for experimental industrial applications. Keywords: Constraints, Dynamic system, Design-of-Experiments, Nonlinear, System-Identification. 1. INTRODUCTION Many industrial applications of nonlinear system identification, such as in aircraft systems and automotive engine calibration, require high efficiency of data capture, high model prediction accuracy and protection from the exceedence operational limits. Dynamic design of experiment (DoE) methodologies are accordingly sought to address these requirements for nonlinear dynamic experimental testing (Ljung, 1999; Klein and Morelli, 2006; Nelles, 2001). Optimality-criteria for DoE based on the variance of parameters have been investigated by many authors. It has been found that the accuracy of parameter estimation has a significant influence on the model prediction error, especially when the model structure is known or can be physically determined. Aoki and Staley (1970) selected inputs using A-optimality, which minimises the trace of the inverse information matrix and briefly indicated the convergence of estimation. Mehra (1974) found that Doptimality, which minimises the determinant of the information matrix, is not affected by the scaling of parameters and showed implementations in both time-domain and frequency domain. In industry, Morelli (1994) demonstrated the practical feasibility of optimal input design in aircraft testing by comparing the optimal signals designed based on A-optimality to non-optimal signals. Zaglauer and Delflorian (2011) proposed a Bayesian modification of the D-optimal design which can be used in dynamic engine calibration in order to prevent bias between the experimental constraints. In this work we argue against the practice of evaluating a test signal design method solely by assessing the level 978-3-902823-14-4/12/$20.00 © 2012 IFAC

118

of accuracy of a set of estimated parameters obtained from the simulation model. In an black box experimental situation, appropriate model structures are in fact, difficult to obtain and are not unique so that assessment based on the covariance of estimated parameters becomes less meaningful. Therefore, if the objective of system identification for experimental based work is to obtain models which make less error in output prediction, minimising the output prediction accuracy might be a more appropriate criterion. Another issue that has been emphasized is that any practical experiment design must also take into account the constraints on the allowable experimental conditions. Goodwin (1971) introduced typical limits that might be found in practice in system identification and applied input and state amplitude constraints to the optimal test signal design for non-linear system identification. Ng et al. (1977) discussed the achievable maximum identification accuracy with a constraint of output power. Ng et al. (1984) considered a problem of optimal input design for parameter estimation for an autoregressive model with constrained output variances. Forssell and Ljung (1998) studied an experiment design problem in identification for control with constraints on both the input and output power and provided an explicit solution. Techniques that can be applied to dynamical systems with practical constraints imposed on input and output amplitudes have been discussed by Morelli (2003). In the Sec. 2, this paper first proposes an iterative procedure of constrained optimal input design for system identification. A new performance index for optimal identification-test-signal design is then introduced in Sec. 3, which weights the parameter covariance by the square 10.3182/20120913-4-IT-4027.00034

CAO 2012 September 13-16, 2012. Rimini, Italy

disturbance, the output at the sample instant t is given by: y(t) = x(t)θ + ǫ(t)

(1)

= [u(t − 1), ..., u(t − a), y(t − 1), ..., y(t − b)]θ (2) +ǫ(t) = f {u, y} θ + ǫ(t)

(3)

where t is the time index from 1 to k, x is the vector of regressor and a + b = n. It can also be presented in a general form: Y = Xθ + E

(4)

T

Fig. 1. Flow chart of the iterative identification process of regressor norm in order to improve the accuracy of the output prediction, and its benefits are discussed. In Sec. 4 an experimentally derived automotive torquefuelling model is used as an illustrative example. The model is used to investigate the ability of the proposed approach to optimize input signals within specified input output limits and thereby obtain improved system identification. The pattern search algorithm is preferred as the optimisation algorithm because its non-stochastic and repeatable characteristics aid the comparisons in the study. The signal obtained with the new performance index is compared with the more typical test inputs used in experimental system identification practice, such as Pseudo Random Binary Sequence signal (PRBS), Pseudo Random Ternary Sequence (PRTS), Uniform Distributed Random Number (UDRN) signals as well as with optimal inputs obtained from conventional performance indices. The model prediction accuracy corresponding to each input type are validated in order to investigate the effectiveness of the new performance index and the proposed optimisation approach. 2. ITERATIVE OPTIMAL INPUT DESIGN WITH EXPERIMENTAL CONSTRAINTS An iterative procedure of constrained optimal input design is illustrated in Figure 1. In the first iteration, an initial model needs to be identified by non-optimal inputs and corresponding outputs. The initial condition for the test signals in the optimisation are taken as the non-optimal signals. The constraints in the current iteration can be conservative and will gradually approach the ultimate experiment constraints in next iterations. The resulting optimal inputs are then applied to the system and another model identified accordingly. The model is updated and the procedure is repeated until an acceptable model is determined using optimal signals obtained with nonconservative constraints. 3. DEVELOPMENT OF PERFORMANCE INDEX For a dynamic time-invariant system which is composed of nonlinear regressors and linear parameters, with a white 119

where Y = [y(1), y(2), ...y(k)] is the vector of output, X = [x1 , x2 , ...xn ] is the matrix of regressors and xi is a k × 1 vector, θ = [θ1 , θ2 , ...θn ]T is the vector of parameters and E = [ǫ(1), ǫ(2), ...ǫ(k)] is the vector of disturbance and ǫ(i) ∼ N (0, σ 2 ); i = 1, . . . , k. 3.1 Information Measurement and parameter covariance Our measure of data information is the Fisher information matrix T   ∂logp(Y |θ) ∂logp(Y |θ) M = EY |θ ∂θ ∂θ     N T ∂y(t) 1 X ∂y(t) = 2 σ t=1 ∂θ ∂θ 

(5)

where the partial derivatives ∂y/∂θ represent the output sensitivity and can be determined from ∂θ ∂f {u, y} ∂y ∂y(t) = f {u, y} + θ ∂θi ∂θi ∂y ∂θi

(6)

Clearly, in an optimized experiment, the system dynamic behaviour needs be excited sufficiently. In other words, the optimal input should maximise the data information. A measure of this information is conveniently taken as a scalar function of the information matrix. The output sensitivity associated with the information matrix indicates the influence of each parameter on the model output. A small change in a parameter will have a considerable effect on the model output when the corresponding output sensitivity is high and vice versa (Klein and Morelli, 2006). According to the Cramer-Rao law, ˆ = E[(θ − θ)(θ ˆ − θ) ˆ T ] ≥ M −1 cov(θ) ( T  )−1 N  ∂y(t) 1 X ∂y(t) = σ 2 t=1 ∂θ ∂θ

(7)

The information matrix can be related to the covariance matrix of the estimated parameter and the equality will be asymptotically achieved if the estimator is efficient. Therefore, optimising the parameter covariance corresponds to optimising the information matrix in some sense. Many scaler functions of the information matrix such as tr(M −1 ) (A-optimal) or log(det(M )) (D-optimal) have been invented for the purpose of optimization.

CAO 2012 September 13-16, 2012. Rimini, Italy

The proposed performance function Jo with respect to the entire data sequence k is accordingly

3.2 Output Covariance In contrast to static models, parameters and regressors in dynamic models often have an influence on each other. Consequently the elements in the parameter covariance matrix are correlated and so minimising a scalar function of the matrix, such as the trace, does not necessarily find a minimum value of each parameter covariance. Moreover, in most experimental identification problems, the true model structure is unknown and in black box modelling a number of different possible model structures may be effective. Generally, the purpose of identifying a black box model is to find a model to simulate the dynamics of the system accurately, but an optimal input design based only on minimising the estimated parameter covariance may not lead to an improved model prediction accuracy since the model structure might be inappropriate. Consider a system described by (4), the ordinary least square (OLS) parameter estimate and parameter variance and output prediction are then given respectively by θˆ = (X T X)−1 X T Y ˆ = σ 2 (X T X)−1 = M −1 cov(θ) Yˆ = X θˆ

(8) (9) (10)

Over the data length k, the output prediction error can be expressed as the covariance matrix of yˆ(i), obtained by a first-order expansion around θ in the form

Jo =

n k X X

xp (i)2 cov(θˆpp )

i=1 p=1

=

k n X X

xp (i)2 cov(θˆpp )

p=1 i=1

=

n X

2

kxp k cov(θˆpp )

(13)

p=1

It can be seen from equation (13) that the approximated output covariance is affected by both the norms of the regressors and the parameter covariance of the estimated model. The regressor norms appear as a weighting to the variances of the parameter estimates. The approach requires significantly less computing time because of the dimensional reduction, especially when the data sequences are extremely long. Initially, any non-optimal signals which are used as initial conditions in the optimization will determine the initial norms. The norms will then be updated iteratively for successive optimal inputs. In practice the test signals and associated norms are found to converge very quickly. The effectiveness of the proposed performance index and this iterative identification process in improving output prediction error is discussed in later sections. 3.4 Measurement of output prediction accuracy

cov(ˆ y (i)) = E[y(i) − yˆ(i))(y − yˆ(i))T ] ˆ ˆ T] = E[(x(i)θ − x(i)θ)(x(i)θ − x(i)θ)

A method of determining the model quality is to check how well the data collected from the system can be reproduced by the identified model. Many criteria can be employed such as mean square error and output fitness (Mathwork, 2009). In this paper, we choose a multiple correlation coefficient function:

ˆ − θ) ˆ T x(i)T ] = E[x(i)(θ − θ)(θ T ˆ = x(i)cov(θ)x(i) = [x1 (i), · · · xn (i)]   cov(θˆ11 ), . . . , cov(θˆ1n )  x1 (i)   cov(θˆ21 ), . . . , cov(θˆ2n )   x2 (i)       ..  .. .. . .   . .  . . xn (i) cov(θˆn1 ), . . . , cov(θˆnn ) n n XX xp (i)xq (i)cov(θˆpq ) =

||Yˆ − Y ||2 R2 (Y, Y¯ ) = 1 − ||Y − Y¯ ||2

(11)

p=1 q=1

where x(i) = [x1 (i), x2 (i), · · · xn (i)] is the vector of regressors at time i. 3.3 An output error based optimality criterion In equation (11), if the regressors are well chosen, the parameters covariances will be uncorrelated. Therefore a sensible approximation of the output covariance at data sample instance i will be given by cov(ˆ y (i)) =

n n X X

xp (i)xq (i)cov(θˆpq )

p=1 q=1

≈

n X

xp (i)2 cov(θˆpp )

where Y = [y(1), y(2)..., y(k)]T is the measured output, Yˆ = [ˆ y (1), yˆ(2)..., yˆ(k)]T is the simulated model output, ¯ and Y = [¯ y (1), y¯(2)..., y¯(k)]T is the mean value of Y . The feasible value of R2 is from 1 to −∞ and R2 → 1 indicates a perfect model. The maximum value of R2 that can be achieved is dependent upon the specific identification case. There is no universal value of R2 which is considered to be good in all applications so that R2 is relative and an acceptance criteria must be judged on a case by case basis. 4. ILLUSTRATIVE EXAMPLE 4.1 I/O equations and constraints An dynamic discrete nonlinear model of a 1.6L Zetec engine obtained experimentally in the University of Liverpool, Powertrain Control Laboratory is presented as an example. The model is determined with the following model structure and parameters

(12) y(k) = θ1 + θ2 u1 (k) + θ3 u1 (k)u2 (k − 10)

p=1

120

(14)

CAO 2012 September 13-16, 2012. Rimini, Italy

U1PRBS

+θ4 u1 (k)u3 (k) + θ5 u1 (k)u3 (k − 3) +θ6 u1 (k − 10)u2 (k − 10) +θ7 u1 (k − 10)u3 (k − 9) 2 U1UDRN

+θ10 u2 (k − 10)u3 (k − 9) + θ11 u3 (k)

U1PRTS

+θ8 u1 (k − 10)2 + θ9 u2 (k − 10)2 +θ12 u3 (k − 2)u3 (k − 4) + θ13 u2 (k − 10)u3 (k) +θ14 u1 (k − 10) + θ15 u3 (k)u3 (k − 5)

6000 4000 2000

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

0

10

20

30

40

50 Time(0.1s)

60

70

80

90

100

6000 4000 2000

6000 4000 2000

U1Jp

6000

+θ16 u3 (k − 1)u3 (k − 10)

4000 2000

+θ17 u2 (k − 10)u3 (k − 5) U1Jo

6000

+θ18 u3 (k − 2)u3 (k − 9)

4000 2000

+θ19 u3 (k − 2)u3 (k − 6) +θ20 u1 (k − 6)u2 (k − 7) U2PRBS

+θ21 u1 (k − 6)u1 (k − 10)

U2UDRN

where u1 (k), u2 (k) and u3 (k) present the fuel pulse width (µS), air-bleed valve and engine speed (RPM), y(k) is the corrupted output response of torque and ǫ(k) is a white noise signal with cov(ǫ) = σ 2 = 80. The sample time is taken as 0.1s. The parameters determined by engine experiments are

U2PRTS

(15)

0.5 0.4

0.6 0.5 0.4

0.6 0.5 0.4

0.6 U2Jp

+θ22 u1 (k − 6)u1 (k) + θ23 y(k − 10) + ǫ(k)

0.6

0.5 0.4

θ = [θ1 , θ2 , ...θ23 ] U2Jo

0.6

= [64.17, −0.025, 0.062, −4.06 × 10−7 , 2.73 × 10−6 ,

0.5 0.4

−0.0069, 8.5 × 10−7 , 1.54 × 10−6 , −327.7, 0.0016,

, 0.0021, 8.71 × 10

−6

, −1.36 × 10

, 0.15]

The experimental input constraints are 2000 ≤ u1 (k) ≤ 6000

(16)

0.4 ≤ u2 (k) ≤ 0.6

(17)

1000 ≤ u3 (k) ≤ 2000

(18)

The sensitivity equations of (15) are determined from (6) as ∂y(k) ∂θ1 ∂y(k) ∂θ2 ∂y(k) ∂θ3 .. . ∂y(k) ∂θ22 ∂y(k) ∂θ23

∂y(k − 10) ∂θ1 ∂y(k − 10) = u1 (k) + θ23 ∂θ2 ∂y(k − 10) = u1 (k)u2 (k − 10) + θ23 ∂θ3 = 1 + θ23

∂y(k − 10) ∂θ22 ∂y(k − 10) = y(k − 10) + θ23 ∂θ23 = u1 (k − 6)u1 (k) + θ23

2000

2000 1500 1000

2000 1500 1000

U3Jo

4.2 Sensitivity equations

U3PRBS

6.4 × 10

−7

U3PRTS

−7

U3Jp

−1.56 × 10−6 , −2.25 × 10−6 , 0.029, 5.67 × 10−7 ,

U3UDRN

−1.52 × 10−5 , −4.22 × 10−8 , 0.069, −0.014,

(19)

1500 1000

2000 1500 1000

2000 1500 1000

Fig. 2. Inputs of the MISO model 4.3 Parameter covariance vs proposed performance index

(20) (21)

(22) (23)

121

The A-optimality and new criterion are represented by the perfomance indices Jp = tr(M −1 ) n X 2 −1 kxp k Mpp Jo =

(24) (25)

p=1

−1 where xp is the pth regressor and Mpp is the pth diagonal element of inverse information matrix. Optimal inputs Up and Uo are generated by minimizing Jp and Jo with global

CAO 2012 September 13-16, 2012. Rimini, Italy

algorithm pattern search. Figure 2 gives an example of non-optimal inputs and optimal inputs. The evaluation of the results of all of the different types of inputs against Jp and Jo are shown in table 1, where it can be seen that the Jo -optimal inputs have the smallest values in each case. Table 1. Evaluation of performance indices on original model Inputs PRBS PRTS UDRN Up Uo

Jp 3.76×1012 σˆ2 79.13σˆ2 231.73σˆ2 25.2σˆ2 49.18σˆ2

Jo 1.18×1017 σˆ2 3.85×105 σˆ2 1.02×106 σˆ2 6.22×105 σˆ2 1.97×105 σˆ2

By applying inputs of the type in table 1 to the original model, five different models are identified and these are then validated by signals of each type, to give the results shown in table 2. It can be seen that the models identified by optimal input Uo produces a better R2 than all other models in all cases. The results also highlight the known general unsuitability of PRBS inputs for nonlinear identifications (Leontaritis and Billings, 1987). 4.5 Optimisation with additional output constraints In experiments where the outputs have operational magnitude limits, the amplitudes of the inputs often need to be adjusted in order to satisfy the output constraints. For non-optimal inputs, scaling the amplitude of inputs manually can be an effective approach in linear models. The input-output relationship in nonlinear systems is not that simple, making the process of tuning the input amplitude much more difficult. However, global constrained optimization can help to find solutions to the inputs quickly. In this example problem, the appropriate output constraint and associated input constraints are 1800 ≤ u1 ≤ 6200

(26)

0.35 ≤ u2 ≤ 0.65

(27)

800 ≤ u3 ≤ 2200

(28)

−10 ≤ y ≤ 40 (29) For the non-optimal inputs feasible values are found by trial and error by test on the original model. Models for simulation are identified accordingly and table 3 shows the assciated validation results. Optimal signals satisfying the constraints are found by the constrained optimisation algorithm. The optimal inputs designed by the proposed Jo -performance index lead to a better R2 than any other inputs. Table 2. R2 of models identified by input constraints Model

Model

2 RP RBS

2 RP RT S

2 RU DRN

2 RU p

2 RU o

MP RBS MP RT S MU DRN MU p MU o

49.37% 27.99% 36.73% 43.59% 51.48 %

-62.31% 53.61% 52.57% 51.37% 60.22%

-38.15% 53.84% 47.17% 52.6617% 56.95%

-9.4% 66.42% 55.54% 72.59% 76.88%

-98.75% 52.46% 51.7% 68.65% 79.33%

5. CONCLUSIONS

4.4 Model identification and validation

MP RBS MP RT S MU DRN MU p MU o

Table 3. R2 of models identified by additional output constraints

2 RP RBS

2 RP RT S

2 RU DRN

2 RU p

2 RU o

-284.24% 65.45% 47.38% 44.46% 74.16%

-654.77% 55.01% 46.82% 52.92% 61.59 %

-1048.3% 37.22% 33.06% 41.46% 45.08%

-442.4% 55.22% 36.2% 56.96% 60.42%

-278.99% 74.13% 59.47% 60.07% 75.7%

122

A proposed criterion and assocaited performance index is proposed for optimal inputs design based on the weighted parameter covariance. The criterion is proved to be effective in reducing model prediction error. Optimizations performed subject to input and output amplitude constraints have been studied. The proposed iterative procedure could be of use in the industrial application of constrained optimal test signal design for black box nonlinear dynamic system identification such as in experimental engine testing where high fidelity dynamic models must be acquired subject to safety based operational constraints. REFERENCES Aoki, M. and Staley, R. (1970). On input signal synthesis in parameter identification. Automatica, 6, 431–440. Forssell, U. and Ljung, L. (1998). Identification for control: Some results on optimal experiment design. In Proc. 37th IEEE Conf. on Decision and Control. Goodwin, G. (1971). Optimal input signal for nonlinearsystem identification. In Proc. IEEE, 7, 922–926. Klein, V. and Morelli, E. (2006). Aircraft System Identification: Theory and Practice, chapter Experiment Design, 289–329. A. I. A. A., 1801 Alexander Bell Drive, Reston, VA. Leontaritis, I. and Billings, S. (1987). Experimental design and identifiability for non-linear system. International Journal of Systems Science, 18(1), 189–202. Ljung, L. (1999). System Identification: Theory for the User, chapter Experiment Design, 408–454. Prentice Hall PRT, New Jersey. Mathwork (2009). System Identification Toolbox. MATLAB2009a. Mehra, R. (1974). Optimal input signals for parameter estimation in dynamic systems-survey and new results. IEEE Tran. Aut. Cont., 19(6), 753–768. Morelli, E. (1994). Flight test validation of optimal input design using pilot implementation. In 10th IFAC Symposium on System Identification. Morelli, E. (2003). Multiple input design for real-time parameter estimation in the frequency domain. In 13th IFAC Symposium on System Identification. Nelles, O. (2001). Nonlinear System Identification. Springer-Verlag, Berlin. Ng, T., Goodwin, G., and Payne, R. (1977). On maximal accuracy estimation with output power constraints. IEEE Tran. Aut. Cont, 22(1), 133–134. Ng, T., Qureshii, Z., and Cheah, Y. (1984). Optimal input design for an ar model with output constraints. Automatica, 20(3), 359–363. Zaglauer, S. and Delflorian, M. (2011). Bayesian d-optimal design. In 6th Conference Design of Experiments in Engine Development.