Optimal Input Design Using Linear Matrix Inequalities

Copyright © IFAC System Identification Santa Barbara, California, USA, 2000

OPTIMAL INPUT DESIGN USING LINEAR MATRIX INEQUALITIES Kristian Lindqvist Hakan Hjalmarsson

Department of Signals, Sensors and Systems, Royal Institute of Technology S-100 44 Stockholm, Sweden kristian. lindqvistQs3.kth. se, hakan.hjalmarssonQs3.kth.se

Abstract: The issue of experiment input spectrum design is studied in the framework of optimal experiment design for system identification. A main result is that, under the assumption that the identified parameters are sufficiently close to their true values, we show this problem to be convex in the input spectrum with linear matrix inequality constraints. Thus it is guaranteed to have a global solution (if feasible) and we can solve it with any demanded accuracy. We apply the results in the area of identification for control. The optimization problem that results is a weighted trace problem. We contribute a scheme for evaluating the weight matrix experimentally. This result indicate that the procedure for input spectrum design can be used for control design schemes previously unconsidered for optimal input spectrum design . Two numerical examples are given. These support the approximations taken in the theoretical approach and also show that great advantage may come from optimizing the input spectrum for the identification experiment. It is the authors' opinion that most gain comes from optimizing the input for systems with resonanse peaks, when there are constraints on both input- and output power spectrum. Copyright @2000 IFA C Keywords: Experiment Design, Identification for Control, System Identification

1. INTRODUCTION

ters). However, this problem can be overcome by two-step and iterative schemes.

The subject of optimal experiment design (OED) has been studied extensively, it involves solving the problem of finding the optimal conditions under which an experiment should be performed. Different aspects of the design may be considered, depending on what affects the outcome (quality) of the experiment most.

Recently, numerical solutions to some OEDproblems have been produced by the use of convex optimization, (Cooley et al. (1998)). The ~ED-problem has been shown to reduce to an LMI-constrained convex optimization problem, thus being solvable efficiently using Interior Point Methods developed in recent years , see Boyd et al. (1994) .

Classical results include many analytical solutions. The books by Fedorov (1972) and Goodwin and Payne (1977) give good insight in the type of results that can be expected. Typically, the results are, in practice, difficult to use since they involve the unknown quantity (e.g. the system parame-

The contributions of this paper are: Firstly, the general criteria of optimal experiment design is approximated into a form enabling the use of LMIconstrained convex optimization in the OEDproblem. Secondly, an example in the area of

1085

identification-for-control is given where we show how the results can be used in practice. Thirdly, we also propose a scheme for an experimental evaluation of the weights needed for the optimization problem, thus reducing the analytical complexity of the method considerably. Last, we present the results in the form of examples with simulations. This paper is composed as follows. Section 2 describes the notation and the basics of system identification. In Section 3, the experiment design problem is described. Section 4 contains the application in identification for control and in Section 5, numerical examples for this problem are presented. Conclusions are given in Section 6.

2. SYSTEM IDENTIFICATION In this Section a brief introduction to system identification is presented along with a presentation of the notation used in this paper. A vast literature can be found on the subject of system identification, e.g., the books Ljung (1987) and Soderstrom and Stoica (1989) .

2.4 The Prediction Error Method (PEM) The PEM procedure (Ljung (1987)) uses the prediction ofthe output to find the parameters O. For this we form y(tIO) = H- 1 (q), O)G(q, O)u(t) + (1H- 1 (q, O))y(t) (the one-step predictor of yet)) and also the prediction errors E(t, 0) = y(tIO) - yet). A natural way to find the parameters 0 is then to minimize the mean square prediction error N

VN(O) =

2

t=l

2.4.1. Properties of the estimates The quality of the estimates of the parameters is commonly measured by their asymptotic properties. That is, when the number of data N grows large, the estimates will belong to some distribution. The properties of the distribution will then determine the quality of the estimates. For the PEM procedure, the asymptotic distribution of the estimates will be Gaussian and is characterized by the mean and variance. It can be shown (under some conditions, see Ljung (1987) ) that (1)

lim E{ON} = 00 N-+oo

2.1 System

J~ooNE{(ON-OO)

The system is a real-life application. This may be a plant we want to control or similar.

(ON-Oo)T}=p

1

p- = AE{w(t,Oo)wT(t,Oo)} w(t,O)

2.2 Model

= :Oy(tIO)

(2) (3)

(4)

Notice that wet, 0) is affine in the input signal u(t) (since y(t,O) is affine in u(t) ), this gives that p-l is linear in the spectrum 4>u(w) of the input u(t).

We will use a mathematical model of the type

yet) = G(q,O)u(t)

LE (t,0)

+ H(q,O)e(t) 3. EXPERIMENT DESIGN

where yet), u(t) and e(t) are the output, input and noise in the system. 0 is a parameterization of the linear time invariant model (G (q, 0), H (q, 0)) and q is the forward shift operator commonly used in the literature. (Le., yet + 1) = qy(t) ). The noise e(t) is assumed to be white and Gaussian distributed with variance A.

For the identification, we need to consider the conditions under which the data will be collected. These considerations are collected under the term experiment conditions. These can be, for instance, Input signal waveform (deterministic input), input signal spectrum (stochastic input), possibility of closing the loop (control) and in that case what controller to use. We denote all the experiment conditions that are at the users disposal by X.

This model has been studied extensively, see e.g. Ljung (1987). We will assume that there exists a parameter 00 such that the model describes the system exactly.

2.3 Identification procedure

When the identified model is used in an application, the quality of the model will affect the performance for the application. For this we use the measure

The purpose of the identification is to find the best possible parameter 0 from a set of input- /output measurements.

where 00 are the parameters of the true system and ON are the (stochastic) parameters resulting

1086

from an experiment with condition X, thus ON = ON(X) (i.e. ON depends on X). See Section 4 for an example of J.

3.3 Formulating the OED as a LMI-optimization problem

The weighted trace optimal design problem is

Since ON is stochastic (due to the noise e( t) present in the data yet)) it is more natural to study the averaged (w.r.t. ON) criteria

min Tr{V pVT} u

s.t. a(w)

wE [-11",11")

where W = VTV is positive semidefinite i.e., V is the Cholesky factorization of W. a(w) and [J(w) are frequency by frequency bounds on the input signal spectrum for the experiment. ( Without bounds, the optimum would be reached by infinite input signal power).

3 .1 Approximation of the general criteria

In general, the criteria J(X, ( 0 ) is a very complicated function of x. However, typically ON is close to 0 0 which allows us to approximate the function J(X, ( 0 ) with a Taylor expension. Expanding the function around the true parameters 00 and involving only terms up to the second order, we get

J(O,Oo) ~J(Bo,Bo)

:s ~u(w) :S [J(w),

This may be reformulated as (see e.g. Cooley et al. (1998) ) min Tr{Z}

Z,u

s.t.

Z - V PV T

>0 a(w) :S 1>u(w) :S [J(w),

wE [-11",11")

Now, using the Schur-complement, we may rewrite this as min Tr{Z}

+ (B - Bo)T J'(Oo,B o)

+ ~ (B - Bo)T J"(Bo, Bo) (() - Bo)

Z,u

where J'(Bo,B o) = toJ((},(}o) 10=0 and J" is the second derivative in the same fashion. This gives

s.t .

0

a(w) :S ~u(w) :S [J(w),

J(x, ( 0 ) ~ J(Bo, ( 0 ) (5)

3.2 Solving the

= 0 i.e., a

~ED-problem

1 = 2Tr {P(X, (}o)W((}o)}

(6)

P(X,Bo) = E { (ON (X) - Bo) (ON(X) - (}o) T} (7)

W(Oo) = J"(B o, Bo)

wE [-11",11")

The above optimization involves infinite dimensional constraints. This is generally called semiinifinite optimization. A brief overview of the techniques for overcoming this problem is found in Wu et al. (1997). We will approximate the above problem, by sampling the (infinite-dimensional) constraints, giving us a standard optimization problem .

First, we rewrite the approximated criteria (5), omitting the constant term J(Bo, Bo)

J(X, (}o)

(9)

where R = p-l and since R is affine in 1>u(w), this is a convex optimization problem with LMI constraints. Thus this may be solved efficiently, to demanded accuracy, using recently developed interior point optimization methods, see Boyd et al. (1994). Also, a global solution is guaranteed. However, we still suffer from the fact that the weight matrix W depends on the true parameters, thus any approximation of W will lead to errors in the optimized input.

+ ~ E { (B - B0) T J" (B 0, B0) (() - ()0) } (with the assumption that E{ON - Bo} consistent estimate )

[;T~] ~O

(8) 4. IDENTIFICATION FOR CONTROL

Thus, we have the well known weighted trace design problem (equation 6) . This is known as a problem that is difficult to obtain analytical solutions for, see Goodwin and Payne (1977) and Fedorov (1972) . Only in very special cases can these be found. However in our setting, with p - 1 affine in the input spectrum (3 ), and only regarding the input signal spectrum as the experiment design variable(s) to optimize over, we can numerically find the solution easily as shown below.

In this part we will use the above framework to address one of the issues in identification for control. In identification for control, the aim is to unify the identification and control design procedure. We will only consider model-based control design schemes where first an identification experiment is performed and secondly a controller is designed

1087

from the model that results from the identification and is used to control the system. We will address the issue of choosing the conditions under which the identification experiment should be performed. We will solve this issue for a given controller design, assuming the use of a statistically efficient identification procedure (such as PEM).

is the fundamental problem of OED, the weight matrix here involves the true parameters. We may heuristically view this as a-priori knowledge and assume that the weights will depend slowly on variations in the parameters OD.

We wish to control an unknown plant [G(q, OD), H(q, OD)] using a controller designed from identified parameters Cid(q) = C(q, eN) .

4.1.1. Experimental evaluation the Weight Matrix There is an advantage (e.g., we do not need the derivatives of the controller) if we could experimentally evaluate the weight matrix. This can easily be done - with experiments in the form of simulations - by varying the system parameters round their true values (here, also a-priori knowledge is assumed). W

The optimal controller is the controller designed from the true parameters Co(q) = C(q, OD). Closing the system with this controller and running it under experiment conditions Xver (the operating conditions for the control application) we get the perfect output Yo(t) = y(t, OD' Xver).

Observing M realizations of the difference yi(t) = yi(t) - y}v (t) , i = 1 .. . M, where i denotes the i'th realization, and calculating the corresponding quality measures Ji = l:~l (yi(t))2 (with i denoting the i'th realization) . Then we may use some procedure to fit this data (w.r.t. W) to the model J i = O[WOi + ei, where ei denotes some disturbance of unknown distribution (e.g., we could us some ellipsoid data fitting method).

Using the identified controller Cid to close the system, we get a perturbed output Yid(t) = y(t,eN,Xver). We measure the quality of the controller Cid by comparing it to the optimal controller Co in the sense of the mean squared difference of the outputs (E{(Yo(t) - Yid(t))2}). Thus, we obtain the following criteria

Jidc(Xid)

k

=

With this approach we also have an interesting feature, a priori knowledge of the expected variance of the estimates should be used for the determination of the weight matrix. We will obtain different W's for different variances E{OiO[}. With proper use, this could be an advantage of this method.

E { (y(t,eN(Xid),Xver) - y(t,Oo,xver)f}

The problem we want to solve is, with above definitions,

minJidc(Xid) Xid

Figure 1 show a block-diagram for the closed loop system.

4.1.2. A procedure for evaluating W The three step procedure to experimentally evaluate the weight follows: • Generate random variables, Oi that are the perturbations of the estimates. These should be generated with care, to resemble the distribution they get from the identification (at least in order of magnitude). • Simulate the system with controllers designed from the perturbed parameters, Oi = (Jo + Oi, and compare to the optimal output in the sense of mean squared difference of the outputs, generating the measurements J i = E{Yo(t) - Yi(t)}. • Do data fitting w.r.t. W to the relation J i :::::: O[WO i . Preferably, divide the measurements in a fitting part and a validation part, such that a good measure of the performance of the approximation can be evaluated from the validation part. For instance, a data fitting method based on convex optimization with linear matrix inequalities can be used.

Fig. 1. The closed-loop system under experimental conditions Xver with controller C designed from parameters eN. We now fit this into the framework from Section 3, resulting in the criteria

Jidc(Xid) :::::: E

{~WeN } ~ ~ Tr {PW}

(10)

where the formula for W involves derivatives of the controller (dC :,9) evaluated at the true parameter values (0 = OD).

J

4.1 Calculating the Weight Matrix W A major difficulty is the calculation of the weight matrix W(Oo) = J"(Oo,Oo). This demands a second order derivative of the cost function J evaluated at the true values of the parameters OD . This

1088