Constructing Continuous-Time Models From Discrete-Time Data

MODELING OF DETERMINISTIC SYSTEMS

Copyright © IFAC Identification and System Parameter Estimation 1982 , Washington D.C . , USA 1982

CONSTRUCTING CONTINUOUS-TIME MODELS FROM DISCRETE-TIME DATA T. B. Cline The Analytic Sciences Corporation, Reading, Massachusetts 01867, USA

Abstract . A continuous-time model must often be identified from sampled experimental data either because the underlying process is inherently continuous-time or because a model is desired for a sample-rate different than the data, as in design of multirate digital controllers. A method is given for obtaining a continuous-time, linear, time-invariant, state-space model from a previously identified discrete-time model. Based on an isomorphism between covariance functions of wide-sense stationary continuous-time and discrete-time processes, the method preserves causality, stability, and minimality. Impulse response, however, is not preserved. Used with the computationally efficient discrete-time canonical variates realization algorithm, the method solves the continuous-time canonical variates realization problem. Several related sampling issues are discussed, and it is shown that for exact reconstruction, the preferred sampling procedure is essentially trapezoidal integration with exponential windowing. Examples illustrate both use and misuse of the method. Keywords. Bilinear transform; canonical variates; identification; parameter estimation; reconstruction; sampling . INTRODUCTION

tion is properly chosen, the continuoustime model so obtained will have desirable properties .

In modeling a random process from experimental data, the goal is often to obtain either a continuous-time model or a discrete-time model corresponding to a different sample-rate than that of the collected data . This paper describes how a linear continuous-time state-space model can be obtained from a time-invariant linear discrete-time state-space model of a wide-sense stationary sequence.

Structure for the discrete-time model can be chosen several ways: using either the canonical variates method (Akaike, 1975), Akaike's information criterion (Akaike, 1974), or the recently developed structure estimator (Purviance and Cline, 1980) which estimates the structure (Kronecker invariants) directly. With the first and last of these, a structure is selected before any parameters are estimated. Using the information criterion requires estimating the parameters of all candidate structures before selecting the best.

One common approach attacks the problem directly, selecting a candidate continuous-time model structure and solving a parameter estimation problem to specify unknown parameters in the structure. When used, this approach raises two difficult new problems. First, good methods do not exist for choosing a model structure when little a priori information is available. Second, in order to evaluate the likelihood function, numerical integration is required which in most applications increases computational requirements dramatically over those for the corresponding discrete-time problem . Alternatively, available methods for selecting structure for a discrete-time model and estimating its parameters can be exploited . The resulting discrete-time model can then be transformed into a continuous-time model. As will be shown below, if the .transforma-

Whether a satisfactory discrete-time model can be directly transformed to a continuous-time model depends on the transformation used. For example, the dynamics matrix of a discrete-time state-space model should be the transition matrix (matrix exponential) of the continuous-time dynamics matrix, evaluated at the sample period . Applying a log transformation to recover the continuous-time dynamics matrix fails when the discrete-time dynamics matrix is singular. In such cases, the corresponding continuous-time process noise intensity matrix may turn out to be indefinite and the continuous-time model so obtained is not physically realizable.

521

522

T. B. eline

This occurs because the standard continuous-time state-space model is not appropriate in such cases; an implicit timedelay is present which must be mode led by a delay-differential equation rather than the standard ordinary differential equation model . The model then is not finite dimensional . Further, the log transformation is not suitable because it does not preserve stability. The basic difficulty with using a transformation which attempts to invert the sampling process is that sampling, as usually modeled, is not one-to-one. This is more apparent in the frequency domain. Let z denote the z-transform variable and s the Laplace transform variable . Suppose we try to invert sampling by replacing z by esT in the discrete-time model's transfer function, where T is the sample period. As is well-known, this results in frequency folding (Oppenheim and Schafer (1975), Kaiser (1963)), a direct manifestation that sampling is a many-to-one mapping . However, if sampling is modeled properly, a one-to-one mapping can be found. If an additional factor is included so that inner products are preserved, the bilinear transformation, used to design digital filters and digital controllers, gives the required isomorphism between continuoustime signals and discrete-time signals (Steiglitz, 1965). The modified transformation preserves stability, causality, and irreducibility. Hence, a minimal realization of a discrete-time process can be transformed into a minimal realization of an equivalent continuous-time process. Further, since the canonical variates method for selecting structure maximizes mutual information between the past and future of a process (Akaike, 1975), the continuous-time canonical variates problem can be solved by transforming the solution of the equivalent discrete-time problem. This obtains because mutual information can be expressed as a function of the process covariance function (Omatu, et al . 1976) which, being an inner product, is preserved by the modified bilinear transformation. The second section gives an isomorphism between covariance functions and transfer function representations of wide-sense stationary continuous-time processes and wide-sense stationary discrete-time processes. In the third section this transformation is used to determine a continuous-time state-space realization from a discrete-time state-space realization obtained by the canonical variates method . Examples show how improperly mode led sampling affects the reconstructed model and show how improper use of the modified bilinear transformation leads to reconstructed models which are not minimal . The final section states conclusions and describes topics for future research.

AN ISOMORPHISM BETWEEN CONTINUOUSTIME AND DISCRETE-TIME PROCESSES

Let {xc(t), t

£(-~,~)}

be a quadratic-mean

continuous wide-sense stationary process with covariance function Rc(l) and spectral density ~c(s).

Similarly, let {xd(k),

k = 0, ±1, ... } be a wide-sense stationary sequence with covariance function Rd(Q) and spectral density

~d(z).

A modified

form of the isomorphism between {xc} and {x } given by Steiglitz (1965)is used: d

~~~(a(z-~) (z+1)2 2

c

(1)

z+l a+s)

( a-s ~~ a -s

(2)

where a is a positive number. A similar form of the bilinear transformation was used by Hitz and Anderson (1972) to obtain an algorithm for solving the continuoustime matrix Riccati equation by solving an equivalent discrete-time matrix Riccati equation. A common choice of a used in digital control and digital filter design is a = 2/T. The isomorphism A will be used in the next section to reconstruct a continuous-time model from a discrete-time model such that the two are equivalent in the sense that second-order characteristics are preserved. ALGORITHM FOR TRANSFORMING TO CONTINUOUS-TIME Because the isometry described in the previous section preserves covariances, stability and casuality, the bilinear transformation can be used with spectral factorization to construct continuous-time models from discrete-time models . However, as will be seen, some subtleties are involved; some factorizations are better than others . Induced Transformation of Linear Filters Steiglitz (1965) showed that the isomorphism described in the previous section induces a bilinear transformation between transfer functions

wc (s)

= W (a+s) d a-s

Thus, a continuous-time model formally equivalent to the given discrete-time model can be obtained by evaluating (3) and constructing a minimal state-space realization of the resulting W (s). But c

(3)

523

Constructing Continuous-Time Models from Discrete-Time Data if care is not taken, this can give a model with some undesirable characteristics . Consider the following example.

As in the SISO case in Example 1, Wc is given by (3) . realization

Example 1. Suppose a single-input single-output (SI SO) discrete-time model is given which has a strictly proper transfer function m

n

K

i=1

(z+a ) i

(4)

Wd (z) = ----=--=-n----

(z-l).2

n

x t +1

Ax t + BUt

(8)

Yt = CX t + Ut where the innovations {u } has spectral density matrix 1. The tfansfer function matrix is

(z+b . )

Wd(z)

J

j=I

Now suppose Wd has the

= C(zI-A) -1 B +

I

(9)

Using (3), the equivalent continuous-time system's transfer function matrix can be written, after some algebra, as

Using (4), the transfer function of the equivalent continuous-time model is

Wc(s) = H(sI-F)

-1

K+ G

(10)

where Since the pole-zero excess of the discretetime system is .2 + n - m, the number of zeros introduced in Wc by transforming Wd equals the pole-zero excess . These new zeros lie in the right-half s-plane at a; hence, the continuous-time model will be nonminimum phase . Even though no feedforward term is present in Wd (it is assumed strictly proper), a feedforward term will always be present in W since W is only proper and its pOle-zefo excesscis zero . If ~d has a feedfo~ard term, as . d~ innovat10ns representat1ons, no nonm1n1mum phase zeros will be introduced in Wc ' Unfortunately, the transformed innovations representation is not an innovations representation for continuous-time. While the form of realization is essentially that of an innovations representation, the driving noise is not white; the model is not a shaping filter driven by white noise . This can be seen by examining the image of

F

a(A+I)-I(A-I)

(lIa)

G

I-C(A+I)-I B

(llb)

H = 2aC(A+I)-I

(llc)

K = (A+I)-1 B

(lId)

The required matrix inverse exists if -1 is not an eigenvalue of A. Since {Yt} is assumed w. s. stationary, this assumption is always satisfied. The model is invertible if a is not an eigenvalue of F . Factoring (10) the time-domain description of the reconstructed process is x = Fx + Kv (12)

Y

Hx + Gv

-1

the spectral density function under A . Suppose $d has been factored in the form

Unlike the innovations sequence {Ut}' the driving noise {v} is not white and has spectral density matr~

(6)

21

where Wd is the transfer function matrix of an innovations representation and 1 is the intensity matrix of the innovations. -1

Applying A ,

2

=~

a -s

$ (a+s) d a-s

= W (a+s) ~ wt(a-s) (7) d a-s 2 2 d a+s a -s

=~

(13)

a -s

The driving noise v can be thought of as the output of a multivariate "first-order Markov" model driven by white noise with spectral density matrix 21 . Thus, the transformed model is not an innovations representation. This~ not as surprising as first appears. In transforming the discrete-time model to continuous-time, a sample and hold operation must be used whose duration depends on a, the frequency warping factor . As a result, each discrete-time noise sample is essentially

T. B. Cline

524

integrated to give the equivalent driving noise for the continuum of times between sample instants, correlating the continuous-time noise. As expected, when a is chosen as 2/T, v is a white noise in the limit as T70. To model the sample-andhold effect with a shaping filter increases the order of the continuous-time model over its discrete-time equivalent by the dimension of the output .

(19)

where as before {u} is a white noise sequence with intensity ~ . The spectral density of the output of f,!d is

Consider the first-order Example 2 . discrete-time system

(20) -1

Using A ax

t

+ bu

t

,

(14)

2

2

a -s 2

bc + 1 z-a

$

(a+s) d a-s

(15)

Cl

/L a+s

z =

Evaluating (11), the equivalent continuoustime model is x

a(~) a+l

y

2ac a+l x +

b

x + a+l v (16)

(1 - ~) a+l

(21)

Thus, the transfer function matrix of the reconstructed continuous-time system is

v

W (s) = /L c[a+sI_A1-IB c a-s a-s 'j (22)

where {v} is the first order Markov process described by v

-ay

+

a,fi

w

(17)

with {w} a unit intensity white-noise. The transfer function for (16) is the same as would be obtained using (3), bc(a-s) (a+l)s-a(a-l) + 1

a+s a+s

-1

where F = a(A+I) (A-I) . A state-space realization for this transfer function is x = Fx + Gv

Mc :

(23) y = Hx

(18)

Note that the zero at a in $ (s) is cancelled by the pole of the sh~ping filter of the noise. Thus, the equivalent process noise is white, as required by an innovationS-representation. However, the feed through (observation) noise is still not white; the shaping filter po le is not cancelled in the spectrum . This suggests an alternate factorization of $ which overcomes the problems of nonminim~ phase and increased dimension and gives the desired minimal realization as a shaping filter driven by white noise.

where F = a(A+I)-l(A-I)

(24a)

G = (A+I) -l B

(24b)

H=,fiC

(24c)

and {v} is a white noise with intensity As before, the required inverse exists if Md is asymptotically stable. Further, it can be shown that if Md is a minimal realization then M is also. As long as Wd is strictly pro~er, Wc does not

matrix~.

Spectral Factorization

have the feedforward terms which gave problems in the previous section .

Suppose the discrete-time model has been obtained by the canonical va r iates method described in Akaike (1975). Then

Suppose the underlying conExample 3 . tinuous-time process is first-order Markov,

525

Constructing Continuous-Time Models from Discrete-Time Data x =

-~x

(25)

+ U

The steady-state variance of the reconstructed state is

where u is a unit intensity white noise. Sampling this process at a rate liT, the discrete-time description is (using the "usual" sampling method)

(34)

(26) OD

where w is a discrete white-noise sequence t with intensity

where P is the steady-state variance of the ori~inally sampled continuous-time process. OD

(27)

The transfer function of the sampled process is (28)

and its spectral density is

(z-e

-~T

)(z

-1

-e

-~T

A desirable characteristic of any reconstruction is that the correct underlying process be recovered in the limit as T~O. In general this is not possible for this example, but the steady-state variances can be made equal for some choices of sampling rate and frequency warp factor. For example, a common choice of Cl in digital control and digital filter design which makes the steady-state error constants invariant is Cl = 2/T. For this Cl, the sample period must satisfy

(29)

)

T( 1+e

Using the modified bilinear transformation to reconstruct a continuous-time model of the original process,

(l-e-~T)/~ '" c (s)

OD

In general, P~ 1 P because the reconstruction does notXinve~t the sampling used, but instead corresponds to a sample-andhold procedure.

(s+ClY) (-s+ClY)

t-'

)

= 2

(35)

Another choice, which allows T to take any positive value, is (36)

(30)

where

(31)

y

-RT

In both cases, the condition for keeping the steady-state variance invariant depends on ~, a parameter of the underlying continuous-time process generating the data. In general then, the sampling model used to get (36) will not give a very robust reconstruction.

A sampling procedure which will permit exact reconstruction, in the second-order sense, can be determined from the sample path version of the isomorphism A. Let Xd(z) denote the z-transform of xd(k) and

A spectral factorization of "'c gives

s+ClY

(32)

J.

For

finite time intervals, an isomorphism between the space of continuous-time signals and the space of discrete-time signals is

which has the state-space realization

X

Xc(s) the Laplace transform of xc(t).

-ClyX + v

~:

{x (t)} ~ X (s) ~ Xd(z) c

c

= CI~ X (CI(Z-l))~ {xd(k)}

where v is a white-noise process with intensity

z+1

c

z+1

(37) (33) lSPE-l- R"

T. B. Cline

526 ~

-1

{xd(k)} ~ Xd(z) ~ Xc(s)

=

~ X (a+s)~

a-s

d

a-s

{x (t)} c

(38) The second factor

Xd{~~:)in

(38) corre-

sponds to a trapezoidal-rule integration and the first factor corresponds to filtering by an filter with exponential impulse response. Hence, reconstruction should use the classical Tustin method with exponential smoothing. Conversely, sampling should discretize integrals using trapezoidal integration and follow with exponential windowing.

CONCLUSIONS Use of the bilinear transform to construct continuous-time models from discrete-time models has been described. As used in designing digital controllers and digital filters, properties of causality and stability are invariant under this transformation. Using a modification suggested by Steiglitz, the transformation also becomes norm-invariant. Norm-invariance enables connecting spectral factorization in the two time domains by the modified bilinear transformation. Using one spectral factorization, discrete-time innovations representations have been shown to be mapped by this transformation into innovationslike continuous-time representations which, however, are not driven by white noise. The transformed driving noise is a multivariate first-order Markov process with bandwidth equal to the frequency warp factor of the modified bilinear transformation. This smoothing of the discrete-time innovations accounts for the sample-andhold operation inherent in constructing a continuous-time signal from a discretetime one.

Using another spectral factorization of the transformed spectral density matrix, discrete-time canonical variates models have been shown to have an invertible and stable continuous-time minimal realization (Eqs. 23, 24). The implicit sampling model is trapezoidal integration with exponential windowing. REFERENCES Akaike, H. (1974). A new look at the statistical model identification . IEEE Trans. Auto. Contr., AC-19, no. 6, pp. 716-722 . -Akaike, H. (1975). Markovian representation of stochastic processes by canonical variables. SIAM J. Contri., 13, pp. 162-173. Hitz, K.L., and Anderson, B.O. (1972) . Iterative method of computing the limiting solution of the matrix Riccati differential equation. Proc. IEEE, ~, no. 9, pp. 1402-14~ Kaiser, J.F. (1963). Design methods for sampled data filters. Proc. First Allerton Conf. Circuit System Theory, pp. 221-236. Omatu, S., Tomita, Y., and Soeda, T. (1976). An alternative expression of th~ mutual information for Gaussian processes . IEEE Trans. Info. Theory, pp. 593-595. Oppenheim, A. V., and Schafer, R.W . (1975). Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, Inc. Purviance, J.E., and Cline, T. B. (1980). Direct structure identification of linear multivariable systems for operating records. Proc. IEEE Conf. Decision and Contr., Albuquerque, NM. Steiglitz, K. (1965). The equivalence of Digital and analog signal processing. Inform Contr., ~, pp. 455-467.