Software of Identification and Modeling for Linear Multivariable Stochastic Systems

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

MODELING OF STOCHASTIC SYSTEMS

SOFTWARE OF IDENTIFICATION AND MODELING FOR LINEAR MUL TIV ARIABLE STOCHASTIC SYSTEMS Ruan Rong-Yao Department of Mathematics, East China Normal University, Shanghai, China

Abstract. The mathematical model of any multivariable linear time-invariant stochastic system with unknown order and parameters can be found off-line by use this software provided that the data are tested according to the technique included in this paper. The optimal feedback control for quadratic performance indices may be realized by the use of a microprocessor or a minicomputer on the basis of the mod el obtained. This paper describes the functions of the software and gives the explanation and the related agreement for use it. Computing formulae and block-diagram for real-time control are also given. Besides, the control efficiency by using this software, which has been applying successfully to two industrial production processes, is simply introduced. Keywords. Identification; mathematical models; multivariable linear time-invariant stochastio systems; parameter estimation; state-space method method; Kalman filters; stochastic optimal control. w

INTrtOnUCTION We consider only the discrete-time linear time-invariant stochastic system because a digital computer is used for testing and controlling. Assume that the system input-output relation can be represented aSI Y(k) - AlY(k-l)+ ••• +Apy(k-p)+BlU(k-l-d) + ••• +B u(k-p-d)+e(k) p

(1)

where u(k) and y(k) are respectively 1 and m dimensional vectors which denote the input and output signals of the system at k instant, A, and B. (i=l, ••• ,p) are mxm and ~

~

matrices recpectively. d is a parameter related to the pure del~ Td of the system

m~l

and sampling time interval T , that is i de [T/Ti where [x] denotes the integer

1,

e(k) - V(k)-Alv(k-l)- ••• -Apv(k-p)

(2)

where measuring noise v(k) and dynamic noice w(k) are white noisss which are independent. The input-output relation of such system m~ also be represented by Fig. 1, where G(z) - ( I -Alz-1 _ ••• -A Z-p)-l

m

G( z)

u

-+y

Fig.l It requires remarkable off-line computation efforts to identify the above system acco~ ding to the test data and establish a mathematical model for optimal feedback control, and also, usually, it might take us a great numb er of efforts for programming and computation on the computer. Hence the software package has been programmed for constructing mathematical models which will reduce greatly the unnecessary repeated work. And it would also be helpful to the popularization and application of the modern control theory. This software has been pro~rammed in both of ALGOI-60 and roRTRAN-IV languages. It mSiY be adopted to design mathematioal models by engineers and technicians or other workers. Even if the user is not acquainted with the programming, he can still obtain the desired mathematical model by this software when he has got the test data according to the explanation for use given below.

part of x. And +Blw(k-l-d)+ ••• +Bpw(k-p-d)

v

THE

EXPLANATION roR USE

In preparation for adopting this software to construot a mathemati.oal model of a certain system, one should not only be justified that the system is linear timeinvariant stoohastio but must also make some

p

( Blz-1+ ••• +Bpz-P) Z-d is a z-transfer function and I denotes the m-th order unit matrix. m 337

338

Ruan Rong-Yao

tes~s

in the ligh~ of oertain requireWhioh proTide the neoessary data for oonstruoting the mathematioal model. Tbe test oan be divided into two stepsl Tbe requirements of the first step arel 1. To find the admis.ible range of variation of the input-output variables; 2. To verify Whether the 8Ystem haS the properties of a linear time-invariant system with stochastio disturbanoes; 3. To estimate pure del~-time 'I'd and the setting time Ts of the ~stem b.1 testing. After doing that, we oan prooeed with the second step for experiment design so that the model m~ be constructed more acourately b,y the test data. It is simpler and oan meet the requirement of aocuracy for engineering to use step and impulse response tests (Both of them are carried out in open looP). For the two Cases the test m~ be well-done if we can determine the sampling time interval Ti and the amplitude a of experiment signal adequately. A period of test must be at least longer than Ts' .en~s

The m sequences m~ be taken as the experiment s~gnal which is superposed on oertain steady-state input value (The steady-state values of input and output of the system an comsidered as zero in calculations) and the impulse response matrix is found by the oorrelation method. The experimental signal is superposed on 1 inputs separately. One is superposed on the i-th input during the i-th test period (i-l, ••• ,l) while other 1-1 inputs are kept at the steady-state value. Tbe test should be oontinuously carried out for 1 periods. Such tests must be repeatedly made for several times under the same and the different conditions in order to choose a set of input-output data from them for identifying and modeling the system. The experiments signals determined acoording to certain experiment design oan also be superposed on all 1 inputs at the same time. If the conditions of identification experiment (Soderstom, Ljung and Gustarsson,1976, Hu and Lu, 1981) are satisfied the inputoutput data of the system operating in closed loop may as well be adopted for identification and modeling which are carried out b.1 the software. Assume that the number of sampling times is q during every test period, then we shall obtain lXq~(l+m) data of the input-output. The input data are given by experiment design, and the output data are obtained b.1 sampling and filter processing. In using this software, we agree to assign the input-output data in array x(lll+m, lllXq], where the first m rows are the output data, and the last 1 rows are the input data. Output must oorrespond to input aocording to the time relation in (1), that is, y{k+d) must correspond to u{k). If a m sequence is used as experiment signal, we need only be given the data of an impulse

response matrix obtained b.1 the oorrelation analysis method, and we agree to a.sign data in arr~ x[llm, lllXql. For step response test, also, we need only be given the output data x[llm, lllxqJ. In the last two oases, the input data will be a.signed automatioally b.1 the oomputer using this software. For a ~stem with arbitrary order p, dimensionalities 1 and m of input and output, this software will work without ~ change no matter which kind of signals is used as input, if only the related Values p, q, 1, m, t, u, v are punched at the beginning of the datum tape. Where q is the number of sampling times within a test period, v,a soale faotor, means that, before every output datum is to join operation, it must be divided b.1 v (Input data is not divided b.1 v). When the test data are step response ones take u-O, else u.l. When the test data are step or impulse response ones, take t-m, else t-l+m. At the beginning of using this software, the order of (I) p may be taken lower (e. g. p_l), beoause this software has suoh a functionl after estimating parameters of a p-th order difference equation (I), it oan automatically use the original data to estimate parameters of a {p+l)-th order one and then, {p+2)-th order and so on. So we Oan obtain several model. with different orders which offers a choice for using. THE FUNCTIONS OF THE SOFTWARE Estimate of Al, •• ,A ,Bl, •• ,B

p

P

and p in (I)

Acoording to the given informations p, q, 1, m, t, u, v we arrange the input-output data into matrioes X and Yl X _ (h(O),h(l), ••• ,h{s»T Y _ (y{O),y(l), ••• ,y{s»T where h{k). (y '1'( k-l), ••• ,y'I' (k-p),uT(k-l-d), 'I' T ••• ,u (k-p-d», Y{k) •

(Yl(~)"",ym(k»T

u(k) • {Ul(k), ••• ,u {k»'l', k - O,l, ••• ,s, s - lq-l, then the least square estimate of paraaeter eT _ (Al, ••• ,Ap,Bl, ••• ,B ) p

~LS - (XTX)-lxTy. If disturbanoe e in (l) is a Gaussian White noise, Ljung (1916) proTed that the least square esiimate @LS of e is oonsistent. r.

When e is a ooloured noise, 918

.~

not be

consistent. The instl'Ufental vUiable met.h od , Take 18

i . adopted to improve 9

Software of Identification and

'f Z - (z{O),z{I), ••• ,s{s» as the instrumental variable matrix, where z(k) _ (yT(k_I), ••• ,yT{k_P),uT{k_I), ••• ••• ,u'r (k-p» T , k-O,l, ••• ,s. and y(k-l), •.• ,y{k-p) are obtained b.1 the least square fitting. Then the instrumental variable estimate of e is A T -1 T SIV - (Z X) z Y. Using this software, we can obtain the least square estimate 9 and the instruLS A mental variable estimate SIV of a at the same time. In the rest of this paper, AI'" •• ,A p ,B , •• ,B p represent values est1mated l b.1 the instrumental variable method. When estimating the parameter e, we assume that the order p in (1) oan take on several different values. And we estimate the several sets of parameter values 9 (with the IV different orders) at the same time so that we oan get several models (I) with different orders. Then we oan determine the order p of model (I) according to Akaike (1974) Information Criterion, that is, choose p suoh that

Model~ng

339

the correlation function of e, that 1s s R(t) __1_ f(k)f'f(k_t), e-p kat

L

t.O,!I, ... ,!p when It I >p, R(t).O. Since the oovarianoe RaE[V{k)V'r(k») of measuring noise v shall be a diagonal matrix, from R(p)--A R, we p

obtain the estimator of R as follows {RJ i i - - [A;IR(P)] 11 (for i.j) or (for i~j). where £A]ij denotes the (i,j)-th element of matrix A.

°

From (2) and (4), it is derived that an approximate estimator of oovarianoe matrix Q.Elw(k)wT(k)] is

Q•

p

(iiTii)-liiTk~JR{k)-iliT.] B(ilTii)-1

where p

1- L

k-l

p

Ap -I m, B - I: Bp • k-l

Computation of A and B To adopt the state-spaoe method to make system analysis and synthesis, we take

AIC p - -210gP{e/eIV )+2pm(1+m) is to be minimum, where P(e/a) is a likelihood function of model (I). The order p estimated in this ~ m~ be a little higher. In practice, it is better to adopt control effeot as criterion to justify whether the estimate of order is aoourate. But it is not eas.y to establish a mathematioal model for dynamic system. Usually one oan not construct several models, and then compare one with another. However, it is eas.y to get several models with different orders b.1 this software. Therefore we oan try out one by one in praotioe and ohoose from the effeot the optimal or suboptimal model. Of oourse, the models to be tried out must be suitably ohosen beforehand in order to reduoe the number of tryings out.

T

X{k)-(x.{k), ••• ,xpm (k» as the state va.I. riable of the s.ystem, where T

(~(k)"",xm{k»

(Xi + l(k), ••. ,xi +m (k» ••• ,xi(k+l»

T

T

- (Xi + 1-m (k+l),

-HjU(k-d)-Bjw(k-d)

i • m,2m, ••• ,{p-I)m, HI • Bl ,

- Y(k)-v{k),

j . 1,2, ••• ,p-l,

H2 • B2+A 1HI ,···,

Hp - Bp+AIHp-l+ • . •+Ap-1B I then (I) and (2) are equivalent to the following state and measuring equationsl x(k+l) - Ax(k)+Bu(k-d)+Bw(k-d) Y(k) • Bx(k)+v(k)

Estimate of Q and R

(5) (6)

where Let Y{k) - AlY{k-I)+ •.• +Apy(k- p ) +BIU(k-I-d)+ ••• +BpU(k-P-d)

A.

(3)

Ap Al" Al , p-

f(k) - Y(k)-Y(k), R(t) ., Ele(k+t)eT(k)] then when SIV is equal to the true value of parameter e, subtracting (3) from (1), we obtain f(k)Ke(k), which is am-dimension weakly stationar,y stoohastic process and usually satisfies the ergodio condition. Henoe we use f(k) to estimate approximately

[~m.~ ..~m~ ...... ... ~ml ~

here O.and I . denote respeotively •• m null ~d

unit matrioes.

Computation of Kalman Filter Gain Matrix K

Ruan Rong-Yao

340

Let ~{k/j} denote the optimal linear estimate of x{k} based on y{O},y{l), ••• ,y{j}, then the optimal filtering of the system defined by {5} and {6} is ~{k} _ ~(k/k-l)+KItY{k), (7) where x{k/k-l) _ AX{k-l)+Bu(k-l-d) and 4y(k) - y(k)-H£(k/k-l) are the one-step optimal prediction and the innovations obtained from the k-th measuring Y(k). Kk in (7) is determined by the minimal variance criterion; that is, we should choose one such that the matrix

is to be minimum. This software adopts recursive formulae TAT Pk / k _ l - APk_IA +BQ'B (8) Kk - Pk/k_IH (HPk/k_lH +R)-l { Pk - Pk/k_l-~HP to find the optimal filter steady-state gain matrix K by cyclically iterative method where P k/k-l- E [x(k/k-l )x(k)] [x{k/k-l )-x{k)]

~ 5

The initial value Po may be taken as 10 1 n (n=pm). The number of iterating times is determined appropriately according to the converging speed in putations. If Kalman filter of the system is stable, the matrix sequence t~} should be convergent and the limit, lim KkcK, is the optimal filter steady-state

g~in

matrix to be found.

For set value control problem, generally, one needs only to consider the accuracy of steady-state filtering. Henoe K can replace Kk in (7) approximately, and the on-line recursive computation of Kalman filtering can be greatly reduoed (The computation of formulae (8) is omitted). Computation of L, L , Mi for Optimal d Feedback Control We adopt a quadratic performance index J{u). E ~ [xT{k)Q x{k)+uT{k-d-l)u{k-d-lY k.d+l 0 '] {9} As a criterion for choosing stoch3stio optimal control strategy, where

m

-1

-1 _1] m ••• -l

............

-1 -1 •••

m

Also, Q m~ be taken as a nonnegative deo finite matrix with other form (It is determined by the practical requirement of oontroll. In this case, the statement in the software which assigns the given values to Q should be accordingly revised. Weight q o is a scalar and the choioe of it depends on which is to be considered mainly. If we want to emphasize the control speed, q should be given a larger value. In order to make the power variation of the input small, q should be given a small value. To establish optimal stochastic control is to find the admissible oontrol strategy {u(O),u{l), ••• ,u(s-l-d)} for the system described by (5) such that J(u} defined by (9) is minimized. Acoording to the separation theorem (Witsenhawsen, 1971; Here it is generalized as the version of timedelay systems), the stochastic optimal oontrol strategy is given by u(k) .. -L(k)i(k+d/k), k=s-l-d, s-2-d, ••• , 0

{10}

where ~(k+d/k) is an optimal predictor of x(k+d} based on y(O),y(l}, ••• ,y(k) (When d-O,x(k+d/k) is the optimal filtering of x(k» and oan be represented as

x(k+d/k)_Ad~(k)+

cl

L

Ai-lBu(k+i_d_l} (11) i-I and L(k) is given by the recursive equations L(k) - (Il+BTQ(k+l)B)-lBTQ{k+l)A

(12)

R{k+l) .. Q(k+l)(A-BL(k»

(13)

Q(k) - Qo+ATR(k+l)

(14)

with the initial condition Q(e-d)-Q • When o

s is sufficiently large, L(k)-+L (a constant matrix) as k-+O, and L may reoursively be obtained from (12)-(14). The number of iterating times is given automatioally in using the software. Matrix L thus obtained is the optimal feedback coeffioi.nt matrix at steady-state. The value of L depends on the choice of weight q. In oomputing by the software, ten different values of q are being selected so as to compute ten corresponding L's for choosing and using. With the chosen matrix L, we calculate Mi • -LA~, i-O,l, ••• ,d-l Thus, the stochastic optimal control strategy at steady-state is Ld - _LAd,

d-l u(k) - Ldx(k)+.L -i u{k+i-d) 1.-0

(15)

Software of Identification and Mode1ing If the requirement on control accuracy i8 not so high before the system reaohing steady-state, formula (15) m~ also be used to compute control strategy.

341

Y(k+d)- Al(k-l+d)+ ••• +ApY(k-P+d) +BlU(k-l)+ ••• +Bpu(k-p)

(17)

Therefore the prediction error is Ay(k) • y(k+d)-Y(k+d)

Simulation of Feedback Control

.. Al AY(k-l)+ ••• +A p6 Y(k- P )

Take initial values of x(k) and u(k) as x(O)- (1,1, ••. ,1)

T

+B16~+ ••• +Bl'~_p+l·

and

(18)

This error by(k) can be replaced appreximately by innovation lIy(k) in real-time oontrol (see Fig. 2) so as to avoid the calculations in (17) and (18). We can also adopt other methods of approximate calculation for the corrections A~ of input

T

u(-l)_ ••• _ u(-d) - (0,0, ••• ,0) then compute according to the following formulae cyclically and recursevelyl d-l u(k) • Ldx(k)+ L M. u(k+i-d) i-O ~

setting value. Hence we obtain the correction for in put getting v alue as

x(k+l) - AX(k)+Bu(k-d), k-0, ••• ,12 (16)

A~ • ,,( C Ay(k)+C 4y(k-l)+ ••• +

O

Thus we obtain u(k) and x(k) which shows whether or how they approach zero as k becomes large. The simulation is made without considering any stochastic disturbances. However, its result can also reflect on the quality of the feedback control to a certain extent. For use in real-time control, th e better L is chosen based on the simUlation results.

1

Cp4Y(k-P)+D24~_1+···+Dpd~_P+l) ':here

CO_B~l(for I-m) or (BiBl)-lBT (for l~m) Cl--COAl' Ci",-COAi' Di--COBi' (i e 2, ••• ,p) and corre ction s

A~

include an attenuation

factor ~ , 0 <.A ~ 1, which m~ diminish the fluctuations of the outputs, but the speed of correcting devi a tions m~ be slowed down.

Computati on of Coefficient Matrices C i and Di Correcting the Input Setting Value We think the out put devition at steadystate as caused by the deviation ~~ of

COMPUTATION BLOCK-DIAGRAM FOR RBAL- 'l'UE CONTROL

the input setting value which has not been corrected. Hence we adopt the method of correcting the input setting value to make c ontrol system h ave s elf-adapt ability. y(k+d) can be approximately expressed as

A large numb er of computation results can be obtained by the software, most of them offer some reference informations for analysing performance of the system and choosing parameters for real-time control. We need OIJ ly to choose matrices A, B, CO'

y(k+d) "' AlY(k+d-l )+ ••• +Apy(k+d- P ) +Bl(U(k-l)+ ~)+ ••• +Bp(U(k-P)+ ~-p-l)

c l ,···,C p ,D 2 ,···,D p ' K, Ld ,

MO, ••• ,M d_ l each for on e in real-time control, and H and the predictor from a reference model is always taken as H - (1,0 , •.. ,0) • o f the 8,1stem iF m m mxpm After the established Output deviation model has been n e cesea Innovati ()n Computation L Unit del~ Y(k) = ~(k)-z(k) sarily che ched by di4Y(k)·Y(k)-H~(k/k-l) git simulation, the ~ ~eal-time control can I Compute Correction State prediction be carried out accor0 ding to the input si'''k- '[Co'Y(k)+C, ,y(k.-1)+. •• f-' [4 £(k+l/k)-A~(k) 0 ::sc+ gnal computed from +Bu(k-d)-BAU,_ +C PAY(k-p)+D 2 A\l,K- 1+ •• • +l\4~-T'- .], .., Fig. 2. 0 ..... In the block diagram, z(k) is the value ob0 eu . tained from the pra.state f11 tering Correct the setting Cl p. ctical output of the 0 X(k).~(k/k-l)+KAY(k) value ~=~_1-1I~ c+ m system by sampling and c+ Cl and filter processing 'l:I and z(k) is the ideal p. Feedback signal Cl output of the system ..... d-l Input signal at the instant k. For u(k)-Cd~(k)+ Di u(k+i-d) ~ ii( k) -'1c+ u( k) set value control, i.O z(k) should be taken as the desired idenFig. 2 tical value.

I

~

1

1

J-

J

1

T

I

l

1

L

-

Ruan Rong-Yao

342

EXPMPLES OF APPLICATION AND COMPUTATION TIME After working out the software, we have used it to identify the temperature control systems for creep-testing and covertype annealing furnaces and establish ed mathematical models which can be used to realize the optimal feedback control b,y a microproces sor or a minicomputer (Ruan, Xia and Wang, l~el; Ruan and others', I j el) Originally, the temperature control of these two types of furnaces was regulated by routine instruments manually. Each furnace of these two types was regulated in need of three routine instruments. And it aws n ot easy to work and the control accuracy was bad. Now the optimal feedback control of the temperature is carried out by making use of a minicomputer according to the established mathematical model for all furnaces (for tens of furnaces) in a workshop. Mot only the control acouraoy has been raised, but also the labour force and energy have been saved. It is convenient to apply the soft\-Iare to identification and modeling a practical process if it can be approximately considered as a multivariable linear timeinvariant etoohastic system. For example, the temperature control system for covertype furnace was considered a s a single input-single output one at first. But t here was no control to make the temperature inside the furnace well-distributed. So the furnace structure should be improved reasonably and then the optimal feedba ck control for furnace temperature could be realized based on it. However, the structure c ould not be determined once for all. Thus we were required to test several times with di f ferent contro l schemes. The system could be considered as a single input-single output one, it could also be altered as a two input-two output or thr e e input-thre e output even as a three input-five output one. These tests are very different but can be used to develo p several mathematical models, by the same software, which of f ers a choice for use. A variety of items can b ·, computed by this softwar e , \·:e may obtain many results from it. But the program is not long and the computing time is short. For example, we have used a com puter of floating points which may operate 10 5 times per second for modeling a 2nd order system with 3-input and 3-output by this softwa re. If the number of input-output testing data is not greater than 5000, the desired mathematical model ldll be computed completely within half an hour. Again, for a 2nd order model of a system with 2-input and 2-output, it only takes about 12 minutes for computing (When the input-output data is not greater than 1000 ). And, it only takes about an hour for computing all the models of 1st, 2nd, 3rd and 4th order.

CONCLUSIONS This software is of pra ctical use. Using the software, the mathematical models can be found off-line for any linear timeinvariant stochastic systems with unknown parameters and statistics of disturbances (The systems may also be approximately linear and time-delay), if the necessary input-output data are given. This software can be conveniently used in identifying and and modeling, the computing time is not long. According to the model established, the optimal feedback control for quadratic performance indices may be realized b,y the use of a microprocessor or a minicomputer. Furthermore, th e c ontrol system is of selfadaptability and can adapt to variations of of the system parameters and environment perturbances, its operating c an be kept optimality. Seve ral examples of application tions, which used the softwar e to identify and modeling, show that the control effects good. REFERENCES Akaike, H. (1 974). A new look at stochastical mode l identification, IEEE Trans. Autom. Control, ~, 716-723. Hu Yang-Zeng and Lu Wu-Sheng (1981). DesignA of experimental condition of identification of multivariable stochastic system in closed loop, The Sym.p osium Proceedings of the Sino-US Bilateral Meeting on Control Systems, Science Press of China. Ljun g , L. (1976). Consistency of the least squares identification method, IEEE Trans. Autom. Control, 21, 779-781. Ruan Ron~Yao, Xi a 'han-Chi and Wang XinWei (1981). The mathematical model of temperature control for cre ep-t esting furnaces, The Symposium Proceedings of the Sino-US Bilateral Meeting on Control Systems, Science Press of China. Ruan Rong-Yao and others' (198 1). Optimal feedb a ck control for cover-type annealin g furnaces, The Symposium Proceedings of the Sino-US Bilateral Meeting on Control Systems, Science Press of China . Soderstrom,T., Ljung, L. and GUst~88on, I. (1970). Identifiability conditions for line a r multivariable systems operating un de r fe edback, IEEE Trans. Autom. Control, 21, 837-840. Wit senhausen , H. S. (1971). Separation of estimation and control for discretetime sy stems, Proceedings I BEE, 22J 1557-1566.