- Email: [email protected]
Structural Identification and Software Package for Linear Multivariable Systems
STRUCTURAL IDENTIFICATION AND SOFTWARE PACKAGE FOR LINEAR MULTIVARIABLE SYSTEMS K. Furuta, S. Hatakeyama and H. Kominami Department of Control Engineering, Tokyo Institute of Technology, Oh-Okayama, Meguro-ku, Tokyo, Japan
Abstract. The aim of this paper is twofold. One of them is to present the interactive identification package developed by the authors and the other is concerned with the structural identification of linear multivariable systems. The identification package consists of supervisor. program modules and data of systems. The program modules are for the identification of the input-output relation of a linear multivariable system, for the structural identification based on realization, and for data management. The structural identification is done byE-minimal realization after the input-output relation is identified, and the procedure gives satisfactory results even in case of single variable plants comparing with conventional procedures. Keywords. Identification; computer aided design; computational methods; multivariable systems. fied a plant by using the minimal dimensional realization derived from Gopinath [6] and Budin [7] algorithm. This procedure does not give the acturate results in the presence of noise, because it is largely based on the noise free realization methods. Guidorzi[8] identified the observability indices from the measured input-output data by computing a certain criterion function for all possible combination of observability indices
INTRODUCTION This paper is concerned with the structural identification of a linear multivariable system from measured input and output data, and with the interactive identification program package to achieve the identification using the provided identification algorithms with already measured input-output data, to test the identified results and to transform the identified input-output relation in a desirable canonical form.
In this paper, the structural identification is considered to be done in the form of a minimal realization after the input-output relation is identified based on the generalized least squares identification which is generalized the algorithm of Clarks [9] to be applicable to the multivariable system [10]. where the structure of a model like observability indices is not necessary to be given precisely if they are chosen larger than the minimum values. The identification results obtained with larger values of observabili ty indices is found as accurate as that with the minimum values. For a linear multivariable system, two minimal realization approaches are compared from the view point of the structural identification. One of them is E-minimal realization which is considered to find a system which has a practically controllable and observable system as an equivalent system. The other is based on Gilbert's realization algor~thm. E-minimal realization could give better results for the examples in the paper. For a single variable system, structural identification procedures based on the E-minimal n~alization and AIC proposed by H.Akaike arc compared. For the verification of the results of the structural identification, statistical hypothesis test has been employed.
In the identification of a single variable plant, the number of parameters to be identified is determined by the order, i.e., a canonical form is a minimal representation. This problem is said the structural identification in the identification of a single variable system, and various methods have been proposed [1], [2] and studied their effectiveness [3]. Most of them are depending on the statistical test of residues resulted in the identification with several different order models. In such structural identification procedures, identification of a input-output relation should be repeated with above mentioned different order models, so lots of computation efforts have been required. For a multivariable system, however, the minimum number of parameters to be identified depends not only on the order of system but also on more precise structural information like observability indices. Furuta [4] proposed to practically minimally realize the identified results of the input-output relation which are not necessary to be minimally represented. Liu and Suen [5] identi415
K. Furuta, S. Hatakeyama and H. Kominami
416
In the practical identification including structural determination, most appropriate identification procedure is depending not only on the dynamic characteristics of the plant but also on the characteristics of measured input-output data [11], and can be found only by comparing the identification results obtained by various identification procedures. And until a plant is identified in the desirable form to be used, several kinds of computation must be done. To faciliate the task, the interactive identification package named PIPAC-F, (Parameter Identification Package Developed by Furuta Laboratory) has been developed. This software package has functions classified into following parts: 1) The command line interpreter and data management. 2) Identification of input-output relation from the input-output data contaminated by noise. 3) The minimal realization in a canonical form from the input-output relation. 4) The simulation of the behavior of the identified system for the given input. 5) The statistical test for the order determination and for the time invariance. This package makes it possible for researchers to choose appropriate identification algorithm effectively in the series of identification procedure. The brief description of the program system is explained and the numerical examples in this paper are computed by using PIPAC-F. In Chapter 2, the structual identification problem is discussed. Chapter 3 gives the details of PIPAC-F, and Chapter 4 gives the numerical examples which are calculated by PIPAC-F.
U(i,j)={
~
~~:~~}
Vj>vi (2) Vi-Vj+l Vj2Vi As one of typical identification algorithms, a modified generalized least squares approach is employed. The problem to determine vi for i-th output is same as the problem of order determination for a single variable system, and Guidorzi considered the determination of all {Vi} as the structural identification of a multivariable system. 2nd step) The identified input-output relation is expressed in a state space representa tion. {
x(t+l)= Ax(t)+Bu(t) y(t)= Cx(t)+Du(t)
(3a) (3b)
where state x is n-vector and p
n=i~l Vi 3rd step) The identified state space representation should be minimally realized, which gives the information concerning the observability indices and if the input-output relation should be identified with the properly chosen observability indices, its identification should be repeated with the newly identified observability indices from the first step. The state space representation of the plant is used for control and estimation.
Generalized Least Squares Estimation Algorithm The identification of the input-output relation of (1) can be achieved by using the generalized least squares identification as follows.
STRUCTURAL IDENTIFICATION OF LINEAR MULTI VARIABLE SYSTEM (4)
This section considers the problem of the structural identification in the identification of linear multivariable systems. For a linear multivariable plant with m-input u(t) and p-output y(t), the input and output are measured after contaminated by noise as v(t)=u(t)+m(t),
where 1fT u i (t)=[Yl(t)'···'Yi-l(t)'Yi+l(t),···, T yp(t),u (t)], ai(z
z(t)=y(t)+n(t)
The identification of the plant considered in this paper is achieved in the following steps. 1st step) From the measured data, the inputoutput relation of the plant is identified in the form Vi i-I Vi Yi(t)= k~laiikYi(t-k)-j~l k~U(i,j)aijkYj(t-k) P vi -. ~ k=L(i .)aijkYj(t-k) J=1+l U ,J m vi + j~l k~ObijkUj(t-k). (i=1,2,··· ,p)
-1 -1
ri(Z
vi
)=(l+k~laiikz
vi.
-k
)
)=[-k~U(i,l)ailkz
-k
, ... ,
vi
-k
vi
-k
-k~U(i,i-l)ai(i-l)kz -k~U(i,i+l)ai(i+l)kz
, ,
vi -k vi -k -k~U(i,p)aipkz 'k~Obilkz , ... , vi -k k~Obimkz ] (1)
where Yi(t) and u.(t) are i-th and j-th elements of outpul and input, {Vi}l=l are said to be observability indices and {U(i,j)} are defined as
The error of the identification can be written as (5)
Structural identifi c ation and software package
417
(12a)
blmk (12b)
T
~i=[-aiil,-aii2,··,-aiiVi,riO,ril,··,riVi]
ei(t) can be assumed to be generated by the white noise ~i(t) as Zi -k (l+k~lwikz )ei(t) = ~i(t). (6)
+
T
T-
~i= [ZiZi] Zizi(t)
are calculated. To the observability indices {Vi}' following relation exists.
T
[ailj , ···,aipj]=O , [bilj,···,bimj]=O
where Zi is chosen as 2Vi-l. The least squares estimate ~i is given by A
bpmk
(13) j~v i+l,
for
i=1,2,· ·· p.
(7)
using ~i' the least squares estimate of {wik} is obtained by
:fl[
3) The state space representation of (3) is expressed as
A-(8)
[All ... .
~lP1
.
Apl
w
iZi
and ei(t) is calculated from (5) using ~ i' Defining Zi A -1 A -k wi(Z )=(l+kfl wikz ) A
_
-1
B· [ : :
0 1
t 21 0
c· [ 0
A
(14a)
pp
1
where Ei=[ei(t-l),···,ei(t-Z i )]
and multiplying wi(z gives
T
d 1 o •..
0
t Pl 0
1 (14b)
~1
(14c)
where {t!j} is the (i,j) element of T~l and
o ....
) from the left to (5)
0
1
*
A
~ .(t)=Z~(t) - Zi(t)~1.' 1. 1. where
(9)
1
-* -1zi(t)=wi(z )zi(t), A
Z.* (t)=[z-*. (t-l) , •.. ,v-#* . (t-V i )], _1. A1. -1 _# 1. v~*(t)=Wi(Z )vi(t)
A . .=
1.J
The generalized least squares estimate of is given by
~ i= [Z~TZ~]+Z~Ti~(t)
~i
State Space Representation The algorithm to transform (1) to (3) is achieved as follows: 1) After the input-output relation is identified for i=1,2, •• ,p, it is expressed as -1
+ ... +Trz -1
-r
= [UO+Ul Z + .. ·+Urz
] Y(t) -r
]u(t)
(11)
where r is the observability index satisfying r = max{ v iL
-aijj.J(i,j)
[
bill
: •
b iml
o
(10)
By repeating the above procedure until ~ i(t) can be considered as the white noise . Th~ algorithm presented here is considered as the generalization of the algorithm by Cl arks to the identification of single output and multi inputs system .
[To+Tl z
bilVi " himv j
B.= : 1..
0
Guidorzi said that the determination of possibly minimum combination is the structural identification, and this identification was done by calculating the criterion function for all possible combination of the observability indices . This process requires a lot of computational time, and the sta te space representation obtained from the identified input-output relation with possibly minimal observability indices is not always minimal. And the identified state space representation should be minimally realized to be used for controi. In this paper, instead of identifying the input-output relation with possibly minimum observability indices, it is proposed to identify with larger values of observability indices and to minimally realize the identified input-output relation in the third step.
418
K. Furuta, S. Hatakeyama and H. Kominami £-Minimal Realization
x [sI.··· ,sn]
Even when the input-output relation is transformed into a complete controllable and complete observable state space representation, it is not always minimal from the practical point of view. Therefore, the £-minimal realization introducing £-controllability and £observability has been considered by Furuta [4]. The £-controllability (£-observability) is defined as follows. [Definition 1] A linear system (A,B,C,D) is said to be £-controllable (£-observabte) if and only if all eigenvalues of VV T (N N) are larger than £, where V (N) is the controllability (observability) matrix defined by n-l V B,AB, ... ,A B]. (N = [ CT, ATc T , ... , (AT)n-lCT ]T .) Furuta defined £-minimal realization as £-observable and £-controllable. But this definition gives that some equivalent systems to the £-minimal realization may not be £-minimal. This means that two systems with the same input-output relation have different properties in the minimal realization. Therefore, in this paper, £-minimal realization is defined as follows. [Definition 2] A linear system (A,B,C,D) is said to be £-minima!ll !e~lized if it has an equivalent system (A,B,C.D) for which all the eigenvalues of the matrix NTN are equal to 1 and all the eigenvalues of the matrix VVT are larger than £. where N and V are the observab~l~tl ~nd the controllability matrices for ( A,B,C,D), respectively. The state space representation from the identified input-output relation is not necessarily £-minimal realization which can be equivalently transformed into the £-observable and £-controllable system. since the measurement noise may construct the state model with £observable and £-uncontrollable state. According to the definition, the £-minimal realization of the identified state space representation (A.B.C,D) is achieved as follows: 1) Calculate T satisfying NTN = TTT
T
(18)
where £1~£2~£3 . . . . Let £no>£~£n +1. Then the £-minimal realization (AO~Bo,CO.DO) is calculated by _ T_ TA - SAS, BO- S B, O C = CS. DO= D. (19)
3)
o
where S
[SI' s2.···. s no]·
Above all, if the input-output relation is identified. it should be transformed into the state space form of (3). and then, it is £~ minimally realized by the algorithm. The £realization can be considered not only as a practical computation method of Mayne's realization algorithm [13]. but also as the one of of Ho-Kalman's algorithm [14] , since Hankel matrix
def[i~:d,b~AB
•..• ] 2 c~. CA B, ..•
Hn=
give"
H~Hn
p [ " . "n ••.
l
(20)
(21)
pT
where columns of P are normalized eigenvectors 'o f H~Hn. Therefore, by considering that {£.} (j~no+l)are taken negligible, it follows tftat rank H = no. and (AO.BO,CO,D O); of (19) is realiz~d from Rn. Practical Realization Based on Gilbert's Algorithm In the £-minimal realization procedure mentioned above, the input-output relation should be transformed into the state space in the first step. In order to avoid this process. the following realization method directly from the input-output relation is considered and compared with the £-minimal realization. In the following discussion, the system to be identified has distinct poles. Then the identified input-output relation can be written as G(z) d
(15)
where N is the observability matrix. Such T is given by T - [
2)
:Y'
o U~2
1 [ vl' .. ··vn
(22)
]-1
(16)
where U. and v. are eigenvalue and normalized ~igenve~tor of NTN. Let -1 ~ = TA!l' B = TB, (17) C = CT , D = D. And calculate the controllability matrix n-l V=T[B.AB . . . . . A BI. and determine the eigenvalu~~ and the normalized eigenvectors of VVT as VVI'= [sl, ... ,sn]
[ £1
£n
1
where i=l .... ,n . M c.b~ i 1. 1. As far as the value of /lM.11 2 = IIc.1I 2 11b.112 is not zero, the state corres~onding1.to tRe pole A. can be considered to be controllable and oBservable. d Let Pi= IIc.1I bb .1I • and the fraction expansion of (22Y 1.e1.done so that PI
~
P2
~
P3
~
..•
If Pn > P L Pno+l •.. for the small value of p, (2~) is practically considered as no G(z)
l:
i=l
I-A. z I 1.
+
D •
(22) ,
Structural identification and software package The practical minimal realization of (22)'
,. Obt:inl~'i1Y~,J.
B _
[~]
C = [cl, ... ,c no ], D = D The procedure presented is one of approaches like [15] by cancelling out the poles and the zeros in the transfer function. The results obtained by this method are not so good as those by E-minimal realization. Statistical Determination of the Stru c ture In the previous sections, it is assumed that E and p are to be given a priori. The problem of how to choose E appropriately is a key in order to succeed in the proposed identification procedure. One way to choose E,i.e. E is done by choosing no which satisfies no no n 'V y(no) = (i~l Ei ) / (i~l Ei ) = 1.(23) From the graph of y (no) versus no, the approapriate no is determined. Usually, the first no large epsilons, El,E2,' " ,Eno, are quite larger than the rest of epsilons. But, sometimes, such choice of epsilon is not possible. This paper proposes the following method using the statistical hypothesis test of Bartlett for determining the epsilon. In the proposed method, the total measured data are divided into L sections, similar to [17], and the input-output relation of the plant is identified L times using the data of each section. From the input-output relation identified using the data of the k-th section. the epsilons are calculated according to the algorithm of the E-minimal realization , which are denoted by El(k), E2(k), ... , E (k). Since the epsilons n {E i} are non zero eigenvalues of the matrix H~Hn' they do not depend on the data sections as far as the plant is identified accurately. However, because of the measurement noise, the epsilons {E. (k)} depend on the data sections, i.e. the id~ntified epsilons can be considered as the sum of the real epsilons and the data noise. The noise and the ratio of the noise to every non zero epsilon {E } are assumed to i be normally distributed. Then, Eno is determined by using the statistical hypothesis test as follows. The normalized value E.(k) of Ei(k) ~s defined by _ 1 Ei(k) = Ei(k)/E i , and its estimate of the variance is given by 2 1 L A 2 Si L k~l( E i(k)-l) (24) where Ei is theLestimate of E. given by 1 1 E. (25) L kh Ei(k), 1 and Ei ~k) for i=l, ... ,no can be written as Ei (k) = 1 + n (k), i=l, ..• no, i where the noise n (k) is supposed ~o be N(O, 2 2 2 i 2 0.) and 0 1=0 2= • • • =0 • Hence, no can be det1 { 2} ermined from s n no by finding no which i i=l
419
satisfies the hypothesis of 0~= .•• =02o' This statistical hypothesis test is done ~s follows. In order to test the hypothesis, H: O~=O~= •• = 2 0 , X2 which satisfies Pr(x 2~x 2)=a should be m 0 2 0 determined from X -distribution with the m-I degree of freedom, where a is the level of significance, and u
=
(-nnr)
/
( 1
+
(m+!))
3llll
should be calculated , where
r
(/L L-i
m i~l ~
)L-l / Si
r
is defined as
(~ Ls~ i=l
(L-l)/2
1
m(L-l)
If u~X 2 , the hypothesis Ho should be rejected, o and if u
K. Furuta, S. Hatakeyama and H. Kominami
420
[DATA MANAGEMENT] COMMANO 1) SYSIN/IO Input the input-output data to the database 2) SYSIN/MD Input the system data (A.B.C.D) of dis(MC) crete (continuous) system to the database 3) SYSLIS/1 Type out the System list 1 SYSLIS/2 Type out the System list 2 SYSLIS Type out the System list 1 and 2 4) LIST Type out the name list of the systems in the database with the number of input. output and order 5) DELETE Delete the system with its data 6) COpy Copy the system's data to the new system 7) RENAME Rename the system's name 8) TYPED Type out the data of the specified system 9) KILLD Erase the data of the specified system 10) fol)VED Copy a part of the data of the system to that of the other system 11) BYE Release the Package [IDENTIFICATION. REALIZATION AND SYSTEM ANALYSIS] 12) MKIO Store the input-output data of the system to the input of M-sequence 13) NOISE Add the noise to the input-output data 14) NORM Normalize the input-output data 15) GLS Generalized least squares estimate 16) ML Maximum likelihood estimate 17) GB Gopinath-budin's realization method 18) HANKEL Compose Hanke1 matrix from Markov parameters 19) MRLZ Minimal realization 20) MR Mode 1 reduct i 'In 21) RLZ Realize the system in the observable canonical form from the input-output relation 22) EPSRLZ Test the £-controllability and £-observability and £-minima1 realization 23) TRM Find the minimal order input-output relation from the identified input-output relation in practical sence [SIMULATION AND TEST] 24) DIGIT Discretize the system 25) SHLT Execute simulation 26) TEST Test for the residue to check the validity of the order Table 1. List of commands PIPAC-F
r
y(t)= [ 0.5 -0.5 is done.
e,
e, e,
~ $rsTEtf~
: I: 8 : tJ : S I 8 . D.
et , .
I: "
I : 8· I : I t
8 : I. 9 : I. I f
(i
0
1
iju(t)
The system has poles
0 0
0.24l0±0.0738 0.3927±0.1655 -0.2070±0.1089 -0.0049±0.0537 -0.0090±0.2535 -0.9748±0.1235 1. 5363±0. 0679
o o o 1 o -0 . 0885±0.0567 o o
D= -1.0115±0.0705 0.9930±0.0623
0.9935±0.106l] 0.9468±0.0943
8 ' "s · 8 1 9 : fll I' S: I' , -
• •YSTClf4 • I ' 8- ~ , . ' 8- 8 ' 0 ' 8 ' 0 - 9- I ' I ' . - 8- O' • • ••
Fig.l
1 0.6 J ( ) -1 -0.95 x t +
, 0.1477±0.1856 0.7084±0.140l 1.3794±0.2523 B= , 0.9740±0.1186 i-0.2736±0.1665 :-0.9488±0.1300 .-0 . 8803±0 .1042
o
, SY~TEIfZ : I ' S: 8: 13 ' SI S ' I' I : I . 8 -
[~ -~.61U(t)
and is a controllable and observable system with the observability indecies vl=2, v2=1. The M sequence signals (pseudo random binary si~nal) with the amplitude 1 and the period 21 -1 unit time are used as the input test signals to the system. The input and output are measured after contaminated by noises where ratios of the standard deviations of the noises and signals are 10% for the input and 30% for the output. The measured data are divided into 4 sections, and using the data of each section, the identification of the input-output relation is done using generalized least squares estimate by assuming observability indecies being vl=4, v2=3. The identification results obtained are represented in the observable canonical system of 7-th order as ro 0 0 -0.0976±0.0575 o 0 0.0295±0.053l 1 0 0 -0.16l6±0.1697 o 0 0.0089±0.0420 o 1 0 0.2305±0.0998 o 0 0.0105±0.1100 A= 0 0 1 0.1586±0.1265 o 0 -0.05486±0.073l o 0 0 0.0500±0.0493 o 0 -0.2069±0.08363 o 0 0 -0.0652±0.1757 1 0 0.0523l±0.16584 o 0 0 -0.1344±0.0525 o 1 0.3659±0.0749
C='-O
, $Y'ST£lfl ' 1 : J: t ; I t I f 9 , I ' 1 : I' I f
1 0 1.1 0.9 x(t)+ -0.45. 0
(-0.45, 0.55±j 0.444)
ISYSLlS C SYSTEIf L/ST-I OF PIPAC-FZ :1
0
l-~·5
x(t+l)=
Type out of system list
be given in the form COMMAND NAME SYSTEM NEWSYSTEM/O. Above command indicates that the data of SYSTEM are treated and the data generated by the command are stored in the database as Nm,SYSTEM. Thus a command name is corresponding to a program module. The list of commands is given in Table 1. The system list which indicates what systems and data are store is typed out by PIPACF in Fig.l. EXAMPLE I In this example the identification of a linear multivariable system represented by
where mij ±Omij of the (i,j) element of the above matrices denotes that the (i,j) element is identified with mean mij and the interval with 95 % confidence level is given by from mij- Omij and mij+omij. The poles of the mean of the identified system are (-0.4695 0.547l±j O.44l3 0.403 3± jO.4360 0.4534±jO.152l). The epsilons with means and variances calculated by using data of each section are tabulated in Table 2, the values in the parenthesis denote normalized values, and from hypothesis test, the order is determined as 3. The controllable canonical form of an £-minimal realization is obtained to the s ystem with parameters as mea n values. A=
0 -0.4693 , -0.02
I
1 0 1.079 0.855 0 -0.439
C=
0.5360 0.9763 ~ -0.498l -0.9629
D=
;' 1.0115 . 0.9930
0.9935 0.9468 .
0.5412' 0.9922
'0 B= 1 ·0
0 -0.563 1
Structural identification and software package
421 hypotllesis tests
1-Int.
2-Int.
3-Int.
4-Int
Var!lIlco
X.
5.04 '10-' (1.54'10")
'.44
3.14
accept
1.51
5.00 '10" (2.29 ·10 ")
0.7.
5.99
accept
1.42'10"
1.96'10-' (9.67'10")
20.67
7.14
nject
6.36'10"
1.01·10" (2.25 '10")
20.75
'.49
njec:t
2.27'10 ' (9.64-10")
2.44'10 ' (1.04)
2.35'10 '
3.74 ' 10" (6.74 '10 ")
5.80 (1.01)
5.75 (1.00)
5.73
3 1.66 Values in (l.05) parentheses are 4J 6.6110" normalized (4.65 ·10 .1) values by the 4.78·10" mean of epsilon 5) (7.52'10")
1.63 (1.03)
1.48 (9.38-10-')
(9.78·10~)
1.55
1.75'10" (1.23)
1.62'10" (1.14)
1.65 '10" (1.16)
6.53'10" (1.03)
6.54'10-1
test of Hf~l
u
)lean
Table 2 1) 2.34.10 ' 2.36 ' 10 ' (9.94 '10 ") (1.00) Identified epsilons and 2] 5.37 5.99 hypothesis test ! (9.37'10-') (1.05)
7.58 '10. 1
(1.03)
(1.19)
6) 1.23 '10" 2.75 '10" (5.04 '10 ") (1.13)
3.83'10- 1 (1.58)
1.93'10. 1 (7.94'10 ")
2.43'10"
9.43'10" (1.59 '10-')
26.17
11.07
nject
716.78·10" 1.57'10" (7.23'10-1) (1.68)
2.05 '10" (2.18·10")
1.30'10- ' (1.38)
'.38'10-'
2.84 '10" (3.23'10 ")
32.54
12.59
nject
where the poles of this system are (-0.4542, 0.5466±jO.4393). Comparing the results, the identified results are satisfactory. By applying the realization based on Giltert's method, the identified results contain large error in comparison with the E-minimal realization, where the p's calculated for the 7th order identified system are (3.5l2±1.422, 3.5l2±1.422, 0.969±0.550, 0.347±0.936, 0.174±0.459, 0.152±0.475), and it is found difficult the order accurately.
Normalized auto-correlation functions 1.0 -. ___
{F;l(t)} {F;l(t)}
o.S
0.0
EXAMPLE II For a sigle variable system, the identification of the system whose input-output relation is represented by (1+alz-l+a2z-2)y(t)=(bo+blZ-l+b2Z_2)u(t)
o Fig.2
+ ----.:.:."--- ~(t) l+CIZ- 1 +C2 Z-
2
where al=-l.5, a2=0.7, bo=l, bl=l.O, b2=0.5, cl=-l, c2=0.2, ,,=1.8. Gaussian white noises with zero mean and variance 1 are used for the input u(t) and the noise ~(t). Fig.3 shows AIC normalized by N and the loss function VQl/N E~=l~i(t) obtained by the maximum likelihood estimation. From Fig.2 the value of criterion seems to be flat after n=2. So we can roughly guess the order to be 2(true order). AIC takes minimal value at (n,n c )=(3,l), (4,2) ,(2,3), (2,4). Nest we applied the hypothesis test by Bartlett for the results obtained by E-minimal realization similar to the case of example I. The results of the test in Table 3 give the true order 2. From those results, we can say that the Bartlett's statistical hypothesis test method based on E-minimal realization is a powerful method for the structural identification. The identification results with the E-minimal realization from third order models and with the correct structural models are compared in Table 4. CONCLUSION This paper presents a method to identify a
(time)
-o.S
v
10
5
Normalized auto-correlation functions of {€l (t) land {€z(t)} with using identified parameters with using GLS.
L
v
v 4
oL-...............~--'...
o
2
4
AlC
AlC 4
l
4
o L-.~................---'..
o
Fig.3
-
-
Loss function V and AIC with n c =l,2,3.
422
K. Furuta, S. Hatakeyama and H. Kominami
linear multi variable system with structural identification and a software package to faciliate the offline identification. The proposed approach in structural identification is different from the already proposed statistical approaches in that the structural information of a linear multivariable system is used, and satisfactory results are obtained. REFERENCES
Table 3
1
aean: £1
1 2 3 4
1.330 10 1 8.820 9.609 10 -\ 2.422 10 -.
Data of hypothesis test (Example 2.) nonoalized variance : 51i 1. 839 10.1 5.239 10.1 0.7884 0.3915
u
x!
0.95 12.61 12.81
3.84 5 . 99 7.84
test of
Hl- 1
accept reject reject
Table 4 The variations of identified parameters by ML and CLS aethods usiDl first three data sections (Example 2: ~ -1.8)
[1] H.Akaike, "A New Look at the £ - ain. Statistical Model Identification" , n a 2 n- l ( 11. - 2 ) IEEE AC-19, pp.7l6-723 (1974). True GLS CLS ML HL ML CLS [2] T.Soderstrom, "On Model 1.0 1.055 1.058 1.030 1.039 1.039 1.030 h. Structure Testing in System Idento.086 to.096 10.162 ~0.112 tO~625 10.162 tification", Int.J .Control, vo1. -1.5 -1.575 -1.551 -1.518 -1.532 26, pp.1-18 (1977). al -2.073 -2.027 to.057 to . 076 to.151 to.101 to.625 to.134 [3] H.Unbehauen, B.Gohring, "Tests 1.0 1.030 1.046 1.030 1.016 0.150 0.502 h, for Determining Model Order in to.184 to . 170 to_320 to.119 10.255 to.127 Parameter Estimation", Automatica, 0.7 0.756 0.729 0.732 0.705 1.528 1.176 vol.lO, pp.223-244 (1974). aa to.073 to.128 to.164 to.125 d.051 10.292 [4] K.Furuta, "An Application of 0.5 0.495 0.517 0.489 0.495 Realization Theory to Identificaha -0.105 -0 . 027 to.037 to.017 to.090 to . 061 to.S57 to.I21 tion of Multivariable Processes", IFAC Symp.on Identification, pp. aa -1.371 -0.317 to.S11 to.155 939-942 (1973). [5] R.Liu, L.C.Suen, "Minimal h. -0.334 -0.280 to.162 to.167 Dimension Realization and Identifiability of Input/Output System", The value behind t denotes the interval with 9S percent Elec.Eng. Memorundum EE7504 Univ. confidence level of Notre Dame (1975). [6] B.Gopinath, "On the Control of Linear Multiple Input-Output System", Bell System Technical Journal, March, lation in Estimation Models", Automatica, (1971) . vol.ll, pp.537-54l (1975). [7] M.A.Budin, "Minimal Realization of Dis[16] I. Gustavsson, "Comparison of Defferent crete Linear Systems from Input-Output Methods for Identification of Industrial Observation", IEEE AC-16, No.5, pp.395-40l Processes", Automatica, vol . 8, pp.127-l42 (1971) . (1972) . [8] R.Guidorzi, "Canonical Structure in the [17] Y.Kaya, M.Ishidawa, "Test of Goodness of Identification of Multivariable Systems", Fit of a State Equation Model", IFAC Kyoto Automatica, vol.ll, pp.36l-374 (1975). Symp., pp.34l-346 (1970). [9] D.W . Clarks, "Generalized-Least-Square 'CI'$RLZ SYST£H2 SYST£/f3/0 Estimation of the Parameters of Dynamic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2 • • • • • • • • • • • • • • • • • • • • • Model", 1st IFAC Symp.on Identification 3.17 C OSYSTEIf le RSYSTEIf (1967) . o--·e CAHDIDRTES FOP. EPSILON ,·~.o [10] E.G. Gilbert, "Controllability and Obserle 112 . 38083307E I le 2J' $.83253409E 0 vability in Multivariable Control System", le 3J' 1.'9149007E 6 SIAM J.Control, vol.l, pp.128-l5l (1963). le 4J3.78221482E -3 le 1. '0274370E -4 [11] R.lsermann, et.al., "Comparison and Evaluation of Six On-Line Identification and ~IUF IfE THE HO . OF ACCEPTED EPSILOH' NP.3 Parameter Estimation Methods with Three SELECTED EPSILON. ..3782214817330047£ -2 Simulated Processes", Proc. 3rd IFAC Symp. ( 3. ;r) IfATRIX A on Identification, pp.l08l-ll02 (1973). Cl) (2) ( 3) [12] A.J.W.Van den Boom, A.W.M.Van den End en , ( ' ) '.800000000 • 1.800000000 o 8 . 000000800 8 ( 2) -4.713831'30 -I I. OS7Sg291 0 o 9.39Q297S30-1 "The Determination of the Orders of Process ( ;r) -2. 372639520 -2 8.888088880 • -~.$30S344S0 -I and Noise Dynamics", Automatica, vo1.l0, pp. ( 3. 2) IfATRIX . 245-256 (1974). ( I) ( Z) ( ' ) '.800009900 8 8 . 808000090 • [13] D.Q.Mayne, "Computational Procedure for (2) 1.809860D~0 8 -6 . 202103140 -I the Minimal Realization of Transfer Functions t ~> •. 880080900 • 1 ••00860800 8 Matrices", Proc. lEE, vo1.l55, pp.1363-l368 ( 2. ;r> IfATRIX C (1968). ( 1> ( -z> ( ;r) ( ' ) ~ .•2'8'86S0 -I 1.882080$70 8 6.188038730-1 [14] B.L.Ho, R.E . Kalman, "Effective Construc( 2) -4.113821150 -1 -J . BJB9~46S0 B S.999SsJsSD-1 tion of Linear State-Variable Models from - ( 2. 2) IfATRIX D Input/Output Functions", Regelungs technik, c I) ( Z) ( ' ) 9.982481820 -1 9 . 82181'420-1 vol.14, pp.545-548 (1966). ( 2 ) 1 .•84eJ678D • 1.819$7""0 • [15] T.Soderstrom, "Test of Pole-Zero CancelExample of computation by PIPAC-F
,J-
,.·····..··················· ..···0
Abstract. The aim of this paper is twofold. One of them is to present the interactive identification package developed by the authors and the other is concerned with the structural identification of linear multivariable systems. The identification package consists of supervisor. program modules and data of systems. The program modules are for the identification of the input-output relation of a linear multivariable system, for the structural identification based on realization, and for data management. The structural identification is done byE-minimal realization after the input-output relation is identified, and the procedure gives satisfactory results even in case of single variable plants comparing with conventional procedures. Keywords. Identification; computer aided design; computational methods; multivariable systems. fied a plant by using the minimal dimensional realization derived from Gopinath [6] and Budin [7] algorithm. This procedure does not give the acturate results in the presence of noise, because it is largely based on the noise free realization methods. Guidorzi[8] identified the observability indices from the measured input-output data by computing a certain criterion function for all possible combination of observability indices
INTRODUCTION This paper is concerned with the structural identification of a linear multivariable system from measured input and output data, and with the interactive identification program package to achieve the identification using the provided identification algorithms with already measured input-output data, to test the identified results and to transform the identified input-output relation in a desirable canonical form.
In this paper, the structural identification is considered to be done in the form of a minimal realization after the input-output relation is identified based on the generalized least squares identification which is generalized the algorithm of Clarks [9] to be applicable to the multivariable system [10]. where the structure of a model like observability indices is not necessary to be given precisely if they are chosen larger than the minimum values. The identification results obtained with larger values of observabili ty indices is found as accurate as that with the minimum values. For a linear multivariable system, two minimal realization approaches are compared from the view point of the structural identification. One of them is E-minimal realization which is considered to find a system which has a practically controllable and observable system as an equivalent system. The other is based on Gilbert's realization algor~thm. E-minimal realization could give better results for the examples in the paper. For a single variable system, structural identification procedures based on the E-minimal n~alization and AIC proposed by H.Akaike arc compared. For the verification of the results of the structural identification, statistical hypothesis test has been employed.
In the identification of a single variable plant, the number of parameters to be identified is determined by the order, i.e., a canonical form is a minimal representation. This problem is said the structural identification in the identification of a single variable system, and various methods have been proposed [1], [2] and studied their effectiveness [3]. Most of them are depending on the statistical test of residues resulted in the identification with several different order models. In such structural identification procedures, identification of a input-output relation should be repeated with above mentioned different order models, so lots of computation efforts have been required. For a multivariable system, however, the minimum number of parameters to be identified depends not only on the order of system but also on more precise structural information like observability indices. Furuta [4] proposed to practically minimally realize the identified results of the input-output relation which are not necessary to be minimally represented. Liu and Suen [5] identi415
K. Furuta, S. Hatakeyama and H. Kominami
416
In the practical identification including structural determination, most appropriate identification procedure is depending not only on the dynamic characteristics of the plant but also on the characteristics of measured input-output data [11], and can be found only by comparing the identification results obtained by various identification procedures. And until a plant is identified in the desirable form to be used, several kinds of computation must be done. To faciliate the task, the interactive identification package named PIPAC-F, (Parameter Identification Package Developed by Furuta Laboratory) has been developed. This software package has functions classified into following parts: 1) The command line interpreter and data management. 2) Identification of input-output relation from the input-output data contaminated by noise. 3) The minimal realization in a canonical form from the input-output relation. 4) The simulation of the behavior of the identified system for the given input. 5) The statistical test for the order determination and for the time invariance. This package makes it possible for researchers to choose appropriate identification algorithm effectively in the series of identification procedure. The brief description of the program system is explained and the numerical examples in this paper are computed by using PIPAC-F. In Chapter 2, the structual identification problem is discussed. Chapter 3 gives the details of PIPAC-F, and Chapter 4 gives the numerical examples which are calculated by PIPAC-F.
U(i,j)={
~
~~:~~}
Vj>vi (2) Vi-Vj+l Vj2Vi As one of typical identification algorithms, a modified generalized least squares approach is employed. The problem to determine vi for i-th output is same as the problem of order determination for a single variable system, and Guidorzi considered the determination of all {Vi} as the structural identification of a multivariable system. 2nd step) The identified input-output relation is expressed in a state space representa tion. {
x(t+l)= Ax(t)+Bu(t) y(t)= Cx(t)+Du(t)
(3a) (3b)
where state x is n-vector and p
n=i~l Vi 3rd step) The identified state space representation should be minimally realized, which gives the information concerning the observability indices and if the input-output relation should be identified with the properly chosen observability indices, its identification should be repeated with the newly identified observability indices from the first step. The state space representation of the plant is used for control and estimation.
Generalized Least Squares Estimation Algorithm The identification of the input-output relation of (1) can be achieved by using the generalized least squares identification as follows.
STRUCTURAL IDENTIFICATION OF LINEAR MULTI VARIABLE SYSTEM (4)
This section considers the problem of the structural identification in the identification of linear multivariable systems. For a linear multivariable plant with m-input u(t) and p-output y(t), the input and output are measured after contaminated by noise as v(t)=u(t)+m(t),
where 1fT u i (t)=[Yl(t)'···'Yi-l(t)'Yi+l(t),···, T yp(t),u (t)], ai(z
z(t)=y(t)+n(t)
The identification of the plant considered in this paper is achieved in the following steps. 1st step) From the measured data, the inputoutput relation of the plant is identified in the form Vi i-I Vi Yi(t)= k~laiikYi(t-k)-j~l k~U(i,j)aijkYj(t-k) P vi -. ~ k=L(i .)aijkYj(t-k) J=1+l U ,J m vi + j~l k~ObijkUj(t-k). (i=1,2,··· ,p)
-1 -1
ri(Z
vi
)=(l+k~laiikz
vi.
-k
)
)=[-k~U(i,l)ailkz
-k
, ... ,
vi
-k
vi
-k
-k~U(i,i-l)ai(i-l)kz -k~U(i,i+l)ai(i+l)kz
, ,
vi -k vi -k -k~U(i,p)aipkz 'k~Obilkz , ... , vi -k k~Obimkz ] (1)
where Yi(t) and u.(t) are i-th and j-th elements of outpul and input, {Vi}l=l are said to be observability indices and {U(i,j)} are defined as
The error of the identification can be written as (5)
Structural identifi c ation and software package
417
(12a)
blmk (12b)
T
~i=[-aiil,-aii2,··,-aiiVi,riO,ril,··,riVi]
ei(t) can be assumed to be generated by the white noise ~i(t) as Zi -k (l+k~lwikz )ei(t) = ~i(t). (6)
+
T
T-
~i= [ZiZi] Zizi(t)
are calculated. To the observability indices {Vi}' following relation exists.
T
[ailj , ···,aipj]=O , [bilj,···,bimj]=O
where Zi is chosen as 2Vi-l. The least squares estimate ~i is given by A
bpmk
(13) j~v i+l,
for
i=1,2,· ·· p.
(7)
using ~i' the least squares estimate of {wik} is obtained by
:fl[
3) The state space representation of (3) is expressed as
A-(8)
[All ... .
~lP1
.
Apl
w
iZi
and ei(t) is calculated from (5) using ~ i' Defining Zi A -1 A -k wi(Z )=(l+kfl wikz ) A
_
-1
B· [ : :
0 1
t 21 0
c· [ 0
A
(14a)
pp
1
where Ei=[ei(t-l),···,ei(t-Z i )]
and multiplying wi(z gives
T
d 1 o •..
0
t Pl 0
1 (14b)
~1
(14c)
where {t!j} is the (i,j) element of T~l and
o ....
) from the left to (5)
0
1
*
A
~ .(t)=Z~(t) - Zi(t)~1.' 1. 1. where
(9)
1
-* -1zi(t)=wi(z )zi(t), A
Z.* (t)=[z-*. (t-l) , •.. ,v-#* . (t-V i )], _1. A1. -1 _# 1. v~*(t)=Wi(Z )vi(t)
A . .=
1.J
The generalized least squares estimate of is given by
~ i= [Z~TZ~]+Z~Ti~(t)
~i
State Space Representation The algorithm to transform (1) to (3) is achieved as follows: 1) After the input-output relation is identified for i=1,2, •• ,p, it is expressed as -1
+ ... +Trz -1
-r
= [UO+Ul Z + .. ·+Urz
] Y(t) -r
]u(t)
(11)
where r is the observability index satisfying r = max{ v iL
-aijj.J(i,j)
[
bill
: •
b iml
o
(10)
By repeating the above procedure until ~ i(t) can be considered as the white noise . Th~ algorithm presented here is considered as the generalization of the algorithm by Cl arks to the identification of single output and multi inputs system .
[To+Tl z
bilVi " himv j
B.= : 1..
0
Guidorzi said that the determination of possibly minimum combination is the structural identification, and this identification was done by calculating the criterion function for all possible combination of the observability indices . This process requires a lot of computational time, and the sta te space representation obtained from the identified input-output relation with possibly minimal observability indices is not always minimal. And the identified state space representation should be minimally realized to be used for controi. In this paper, instead of identifying the input-output relation with possibly minimum observability indices, it is proposed to identify with larger values of observability indices and to minimally realize the identified input-output relation in the third step.
418
K. Furuta, S. Hatakeyama and H. Kominami £-Minimal Realization
x [sI.··· ,sn]
Even when the input-output relation is transformed into a complete controllable and complete observable state space representation, it is not always minimal from the practical point of view. Therefore, the £-minimal realization introducing £-controllability and £observability has been considered by Furuta [4]. The £-controllability (£-observability) is defined as follows. [Definition 1] A linear system (A,B,C,D) is said to be £-controllable (£-observabte) if and only if all eigenvalues of VV T (N N) are larger than £, where V (N) is the controllability (observability) matrix defined by n-l V B,AB, ... ,A B]. (N = [ CT, ATc T , ... , (AT)n-lCT ]T .) Furuta defined £-minimal realization as £-observable and £-controllable. But this definition gives that some equivalent systems to the £-minimal realization may not be £-minimal. This means that two systems with the same input-output relation have different properties in the minimal realization. Therefore, in this paper, £-minimal realization is defined as follows. [Definition 2] A linear system (A,B,C,D) is said to be £-minima!ll !e~lized if it has an equivalent system (A,B,C.D) for which all the eigenvalues of the matrix NTN are equal to 1 and all the eigenvalues of the matrix VVT are larger than £. where N and V are the observab~l~tl ~nd the controllability matrices for ( A,B,C,D), respectively. The state space representation from the identified input-output relation is not necessarily £-minimal realization which can be equivalently transformed into the £-observable and £-controllable system. since the measurement noise may construct the state model with £observable and £-uncontrollable state. According to the definition, the £-minimal realization of the identified state space representation (A.B.C,D) is achieved as follows: 1) Calculate T satisfying NTN = TTT
T
(18)
where £1~£2~£3 . . . . Let £no>£~£n +1. Then the £-minimal realization (AO~Bo,CO.DO) is calculated by _ T_ TA - SAS, BO- S B, O C = CS. DO= D. (19)
3)
o
where S
[SI' s2.···. s no]·
Above all, if the input-output relation is identified. it should be transformed into the state space form of (3). and then, it is £~ minimally realized by the algorithm. The £realization can be considered not only as a practical computation method of Mayne's realization algorithm [13]. but also as the one of of Ho-Kalman's algorithm [14] , since Hankel matrix
def[i~:d,b~AB
•..• ] 2 c~. CA B, ..•
Hn=
give"
H~Hn
p [ " . "n ••.
l
(20)
(21)
pT
where columns of P are normalized eigenvectors 'o f H~Hn. Therefore, by considering that {£.} (j~no+l)are taken negligible, it follows tftat rank H = no. and (AO.BO,CO,D O); of (19) is realiz~d from Rn. Practical Realization Based on Gilbert's Algorithm In the £-minimal realization procedure mentioned above, the input-output relation should be transformed into the state space in the first step. In order to avoid this process. the following realization method directly from the input-output relation is considered and compared with the £-minimal realization. In the following discussion, the system to be identified has distinct poles. Then the identified input-output relation can be written as G(z) d
(15)
where N is the observability matrix. Such T is given by T - [
2)
:Y'
o U~2
1 [ vl' .. ··vn
(22)
]-1
(16)
where U. and v. are eigenvalue and normalized ~igenve~tor of NTN. Let -1 ~ = TA!l' B = TB, (17) C = CT , D = D. And calculate the controllability matrix n-l V=T[B.AB . . . . . A BI. and determine the eigenvalu~~ and the normalized eigenvectors of VVT as VVI'= [sl, ... ,sn]
[ £1
£n
1
where i=l .... ,n . M c.b~ i 1. 1. As far as the value of /lM.11 2 = IIc.1I 2 11b.112 is not zero, the state corres~onding1.to tRe pole A. can be considered to be controllable and oBservable. d Let Pi= IIc.1I bb .1I • and the fraction expansion of (22Y 1.e1.done so that PI
~
P2
~
P3
~
..•
If Pn > P L Pno+l •.. for the small value of p, (2~) is practically considered as no G(z)
l:
i=l
I-A. z I 1.
+
D •
(22) ,
Structural identification and software package The practical minimal realization of (22)'
,. Obt:inl~'i1Y~,J.
B _
[~]
C = [cl, ... ,c no ], D = D The procedure presented is one of approaches like [15] by cancelling out the poles and the zeros in the transfer function. The results obtained by this method are not so good as those by E-minimal realization. Statistical Determination of the Stru c ture In the previous sections, it is assumed that E and p are to be given a priori. The problem of how to choose E appropriately is a key in order to succeed in the proposed identification procedure. One way to choose E,i.e. E is done by choosing no which satisfies no no n 'V y(no) = (i~l Ei ) / (i~l Ei ) = 1.(23) From the graph of y (no) versus no, the approapriate no is determined. Usually, the first no large epsilons, El,E2,' " ,Eno, are quite larger than the rest of epsilons. But, sometimes, such choice of epsilon is not possible. This paper proposes the following method using the statistical hypothesis test of Bartlett for determining the epsilon. In the proposed method, the total measured data are divided into L sections, similar to [17], and the input-output relation of the plant is identified L times using the data of each section. From the input-output relation identified using the data of the k-th section. the epsilons are calculated according to the algorithm of the E-minimal realization , which are denoted by El(k), E2(k), ... , E (k). Since the epsilons n {E i} are non zero eigenvalues of the matrix H~Hn' they do not depend on the data sections as far as the plant is identified accurately. However, because of the measurement noise, the epsilons {E. (k)} depend on the data sections, i.e. the id~ntified epsilons can be considered as the sum of the real epsilons and the data noise. The noise and the ratio of the noise to every non zero epsilon {E } are assumed to i be normally distributed. Then, Eno is determined by using the statistical hypothesis test as follows. The normalized value E.(k) of Ei(k) ~s defined by _ 1 Ei(k) = Ei(k)/E i , and its estimate of the variance is given by 2 1 L A 2 Si L k~l( E i(k)-l) (24) where Ei is theLestimate of E. given by 1 1 E. (25) L kh Ei(k), 1 and Ei ~k) for i=l, ... ,no can be written as Ei (k) = 1 + n (k), i=l, ..• no, i where the noise n (k) is supposed ~o be N(O, 2 2 2 i 2 0.) and 0 1=0 2= • • • =0 • Hence, no can be det1 { 2} ermined from s n no by finding no which i i=l
419
satisfies the hypothesis of 0~= .•• =02o' This statistical hypothesis test is done ~s follows. In order to test the hypothesis, H: O~=O~= •• = 2 0 , X2 which satisfies Pr(x 2~x 2)=a should be m 0 2 0 determined from X -distribution with the m-I degree of freedom, where a is the level of significance, and u
=
(-nnr)
/
( 1
+
(m+!))
3llll
should be calculated , where
r
(/L L-i
m i~l ~
)L-l / Si
r
is defined as
(~ Ls~ i=l
(L-l)/2
1
m(L-l)
If u~X 2 , the hypothesis Ho should be rejected, o and if u
K. Furuta, S. Hatakeyama and H. Kominami
420
[DATA MANAGEMENT] COMMANO 1) SYSIN/IO Input the input-output data to the database 2) SYSIN/MD Input the system data (A.B.C.D) of dis(MC) crete (continuous) system to the database 3) SYSLIS/1 Type out the System list 1 SYSLIS/2 Type out the System list 2 SYSLIS Type out the System list 1 and 2 4) LIST Type out the name list of the systems in the database with the number of input. output and order 5) DELETE Delete the system with its data 6) COpy Copy the system's data to the new system 7) RENAME Rename the system's name 8) TYPED Type out the data of the specified system 9) KILLD Erase the data of the specified system 10) fol)VED Copy a part of the data of the system to that of the other system 11) BYE Release the Package [IDENTIFICATION. REALIZATION AND SYSTEM ANALYSIS] 12) MKIO Store the input-output data of the system to the input of M-sequence 13) NOISE Add the noise to the input-output data 14) NORM Normalize the input-output data 15) GLS Generalized least squares estimate 16) ML Maximum likelihood estimate 17) GB Gopinath-budin's realization method 18) HANKEL Compose Hanke1 matrix from Markov parameters 19) MRLZ Minimal realization 20) MR Mode 1 reduct i 'In 21) RLZ Realize the system in the observable canonical form from the input-output relation 22) EPSRLZ Test the £-controllability and £-observability and £-minima1 realization 23) TRM Find the minimal order input-output relation from the identified input-output relation in practical sence [SIMULATION AND TEST] 24) DIGIT Discretize the system 25) SHLT Execute simulation 26) TEST Test for the residue to check the validity of the order Table 1. List of commands PIPAC-F
r
y(t)= [ 0.5 -0.5 is done.
e,
e, e,
~ $rsTEtf~
: I: 8 : tJ : S I 8 . D.
et , .
I: "
I : 8· I : I t
8 : I. 9 : I. I f
(i
0
1
iju(t)
The system has poles
0 0
0.24l0±0.0738 0.3927±0.1655 -0.2070±0.1089 -0.0049±0.0537 -0.0090±0.2535 -0.9748±0.1235 1. 5363±0. 0679
o o o 1 o -0 . 0885±0.0567 o o
D= -1.0115±0.0705 0.9930±0.0623
0.9935±0.106l] 0.9468±0.0943
8 ' "s · 8 1 9 : fll I' S: I' , -
• •YSTClf4 • I ' 8- ~ , . ' 8- 8 ' 0 ' 8 ' 0 - 9- I ' I ' . - 8- O' • • ••
Fig.l
1 0.6 J ( ) -1 -0.95 x t +
, 0.1477±0.1856 0.7084±0.140l 1.3794±0.2523 B= , 0.9740±0.1186 i-0.2736±0.1665 :-0.9488±0.1300 .-0 . 8803±0 .1042
o
, SY~TEIfZ : I ' S: 8: 13 ' SI S ' I' I : I . 8 -
[~ -~.61U(t)
and is a controllable and observable system with the observability indecies vl=2, v2=1. The M sequence signals (pseudo random binary si~nal) with the amplitude 1 and the period 21 -1 unit time are used as the input test signals to the system. The input and output are measured after contaminated by noises where ratios of the standard deviations of the noises and signals are 10% for the input and 30% for the output. The measured data are divided into 4 sections, and using the data of each section, the identification of the input-output relation is done using generalized least squares estimate by assuming observability indecies being vl=4, v2=3. The identification results obtained are represented in the observable canonical system of 7-th order as ro 0 0 -0.0976±0.0575 o 0 0.0295±0.053l 1 0 0 -0.16l6±0.1697 o 0 0.0089±0.0420 o 1 0 0.2305±0.0998 o 0 0.0105±0.1100 A= 0 0 1 0.1586±0.1265 o 0 -0.05486±0.073l o 0 0 0.0500±0.0493 o 0 -0.2069±0.08363 o 0 0 -0.0652±0.1757 1 0 0.0523l±0.16584 o 0 0 -0.1344±0.0525 o 1 0.3659±0.0749
C='-O
, $Y'ST£lfl ' 1 : J: t ; I t I f 9 , I ' 1 : I' I f
1 0 1.1 0.9 x(t)+ -0.45. 0
(-0.45, 0.55±j 0.444)
ISYSLlS C SYSTEIf L/ST-I OF PIPAC-FZ :1
0
l-~·5
x(t+l)=
Type out of system list
be given in the form COMMAND NAME SYSTEM NEWSYSTEM/O. Above command indicates that the data of SYSTEM are treated and the data generated by the command are stored in the database as Nm,SYSTEM. Thus a command name is corresponding to a program module. The list of commands is given in Table 1. The system list which indicates what systems and data are store is typed out by PIPACF in Fig.l. EXAMPLE I In this example the identification of a linear multivariable system represented by
where mij ±Omij of the (i,j) element of the above matrices denotes that the (i,j) element is identified with mean mij and the interval with 95 % confidence level is given by from mij- Omij and mij+omij. The poles of the mean of the identified system are (-0.4695 0.547l±j O.44l3 0.403 3± jO.4360 0.4534±jO.152l). The epsilons with means and variances calculated by using data of each section are tabulated in Table 2, the values in the parenthesis denote normalized values, and from hypothesis test, the order is determined as 3. The controllable canonical form of an £-minimal realization is obtained to the s ystem with parameters as mea n values. A=
0 -0.4693 , -0.02
I
1 0 1.079 0.855 0 -0.439
C=
0.5360 0.9763 ~ -0.498l -0.9629
D=
;' 1.0115 . 0.9930
0.9935 0.9468 .
0.5412' 0.9922
'0 B= 1 ·0
0 -0.563 1
Structural identification and software package
421 hypotllesis tests
1-Int.
2-Int.
3-Int.
4-Int
Var!lIlco
X.
5.04 '10-' (1.54'10")
'.44
3.14
accept
1.51
5.00 '10" (2.29 ·10 ")
0.7.
5.99
accept
1.42'10"
1.96'10-' (9.67'10")
20.67
7.14
nject
6.36'10"
1.01·10" (2.25 '10")
20.75
'.49
njec:t
2.27'10 ' (9.64-10")
2.44'10 ' (1.04)
2.35'10 '
3.74 ' 10" (6.74 '10 ")
5.80 (1.01)
5.75 (1.00)
5.73
3 1.66 Values in (l.05) parentheses are 4J 6.6110" normalized (4.65 ·10 .1) values by the 4.78·10" mean of epsilon 5) (7.52'10")
1.63 (1.03)
1.48 (9.38-10-')
(9.78·10~)
1.55
1.75'10" (1.23)
1.62'10" (1.14)
1.65 '10" (1.16)
6.53'10" (1.03)
6.54'10-1
test of Hf~l
u
)lean
Table 2 1) 2.34.10 ' 2.36 ' 10 ' (9.94 '10 ") (1.00) Identified epsilons and 2] 5.37 5.99 hypothesis test ! (9.37'10-') (1.05)
7.58 '10. 1
(1.03)
(1.19)
6) 1.23 '10" 2.75 '10" (5.04 '10 ") (1.13)
3.83'10- 1 (1.58)
1.93'10. 1 (7.94'10 ")
2.43'10"
9.43'10" (1.59 '10-')
26.17
11.07
nject
716.78·10" 1.57'10" (7.23'10-1) (1.68)
2.05 '10" (2.18·10")
1.30'10- ' (1.38)
'.38'10-'
2.84 '10" (3.23'10 ")
32.54
12.59
nject
where the poles of this system are (-0.4542, 0.5466±jO.4393). Comparing the results, the identified results are satisfactory. By applying the realization based on Giltert's method, the identified results contain large error in comparison with the E-minimal realization, where the p's calculated for the 7th order identified system are (3.5l2±1.422, 3.5l2±1.422, 0.969±0.550, 0.347±0.936, 0.174±0.459, 0.152±0.475), and it is found difficult the order accurately.
Normalized auto-correlation functions 1.0 -. ___
{F;l(t)} {F;l(t)}
o.S
0.0
EXAMPLE II For a sigle variable system, the identification of the system whose input-output relation is represented by (1+alz-l+a2z-2)y(t)=(bo+blZ-l+b2Z_2)u(t)
o Fig.2
+ ----.:.:."--- ~(t) l+CIZ- 1 +C2 Z-
2
where al=-l.5, a2=0.7, bo=l, bl=l.O, b2=0.5, cl=-l, c2=0.2, ,,=1.8. Gaussian white noises with zero mean and variance 1 are used for the input u(t) and the noise ~(t). Fig.3 shows AIC normalized by N and the loss function VQl/N E~=l~i(t) obtained by the maximum likelihood estimation. From Fig.2 the value of criterion seems to be flat after n=2. So we can roughly guess the order to be 2(true order). AIC takes minimal value at (n,n c )=(3,l), (4,2) ,(2,3), (2,4). Nest we applied the hypothesis test by Bartlett for the results obtained by E-minimal realization similar to the case of example I. The results of the test in Table 3 give the true order 2. From those results, we can say that the Bartlett's statistical hypothesis test method based on E-minimal realization is a powerful method for the structural identification. The identification results with the E-minimal realization from third order models and with the correct structural models are compared in Table 4. CONCLUSION This paper presents a method to identify a
(time)
-o.S
v
10
5
Normalized auto-correlation functions of {€l (t) land {€z(t)} with using identified parameters with using GLS.
L
v
v 4
oL-...............~--'...
o
2
4
AlC
AlC 4
l
4
o L-.~................---'..
o
Fig.3
-
-
Loss function V and AIC with n c =l,2,3.
422
K. Furuta, S. Hatakeyama and H. Kominami
linear multi variable system with structural identification and a software package to faciliate the offline identification. The proposed approach in structural identification is different from the already proposed statistical approaches in that the structural information of a linear multivariable system is used, and satisfactory results are obtained. REFERENCES
Table 3
1
aean: £1
1 2 3 4
1.330 10 1 8.820 9.609 10 -\ 2.422 10 -.
Data of hypothesis test (Example 2.) nonoalized variance : 51i 1. 839 10.1 5.239 10.1 0.7884 0.3915
u
x!
0.95 12.61 12.81
3.84 5 . 99 7.84
test of
Hl- 1
accept reject reject
Table 4 The variations of identified parameters by ML and CLS aethods usiDl first three data sections (Example 2: ~ -1.8)
[1] H.Akaike, "A New Look at the £ - ain. Statistical Model Identification" , n a 2 n- l ( 11. - 2 ) IEEE AC-19, pp.7l6-723 (1974). True GLS CLS ML HL ML CLS [2] T.Soderstrom, "On Model 1.0 1.055 1.058 1.030 1.039 1.039 1.030 h. Structure Testing in System Idento.086 to.096 10.162 ~0.112 tO~625 10.162 tification", Int.J .Control, vo1. -1.5 -1.575 -1.551 -1.518 -1.532 26, pp.1-18 (1977). al -2.073 -2.027 to.057 to . 076 to.151 to.101 to.625 to.134 [3] H.Unbehauen, B.Gohring, "Tests 1.0 1.030 1.046 1.030 1.016 0.150 0.502 h, for Determining Model Order in to.184 to . 170 to_320 to.119 10.255 to.127 Parameter Estimation", Automatica, 0.7 0.756 0.729 0.732 0.705 1.528 1.176 vol.lO, pp.223-244 (1974). aa to.073 to.128 to.164 to.125 d.051 10.292 [4] K.Furuta, "An Application of 0.5 0.495 0.517 0.489 0.495 Realization Theory to Identificaha -0.105 -0 . 027 to.037 to.017 to.090 to . 061 to.S57 to.I21 tion of Multivariable Processes", IFAC Symp.on Identification, pp. aa -1.371 -0.317 to.S11 to.155 939-942 (1973). [5] R.Liu, L.C.Suen, "Minimal h. -0.334 -0.280 to.162 to.167 Dimension Realization and Identifiability of Input/Output System", The value behind t denotes the interval with 9S percent Elec.Eng. Memorundum EE7504 Univ. confidence level of Notre Dame (1975). [6] B.Gopinath, "On the Control of Linear Multiple Input-Output System", Bell System Technical Journal, March, lation in Estimation Models", Automatica, (1971) . vol.ll, pp.537-54l (1975). [7] M.A.Budin, "Minimal Realization of Dis[16] I. Gustavsson, "Comparison of Defferent crete Linear Systems from Input-Output Methods for Identification of Industrial Observation", IEEE AC-16, No.5, pp.395-40l Processes", Automatica, vol . 8, pp.127-l42 (1971) . (1972) . [8] R.Guidorzi, "Canonical Structure in the [17] Y.Kaya, M.Ishidawa, "Test of Goodness of Identification of Multivariable Systems", Fit of a State Equation Model", IFAC Kyoto Automatica, vol.ll, pp.36l-374 (1975). Symp., pp.34l-346 (1970). [9] D.W . Clarks, "Generalized-Least-Square 'CI'$RLZ SYST£H2 SYST£/f3/0 Estimation of the Parameters of Dynamic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2 • • • • • • • • • • • • • • • • • • • • • Model", 1st IFAC Symp.on Identification 3.17 C OSYSTEIf le RSYSTEIf (1967) . o--·e CAHDIDRTES FOP. EPSILON ,·~.o [10] E.G. Gilbert, "Controllability and Obserle 112 . 38083307E I le 2J' $.83253409E 0 vability in Multivariable Control System", le 3J' 1.'9149007E 6 SIAM J.Control, vol.l, pp.128-l5l (1963). le 4J3.78221482E -3 le 1. '0274370E -4 [11] R.lsermann, et.al., "Comparison and Evaluation of Six On-Line Identification and ~IUF IfE THE HO . OF ACCEPTED EPSILOH' NP.3 Parameter Estimation Methods with Three SELECTED EPSILON. ..3782214817330047£ -2 Simulated Processes", Proc. 3rd IFAC Symp. ( 3. ;r) IfATRIX A on Identification, pp.l08l-ll02 (1973). Cl) (2) ( 3) [12] A.J.W.Van den Boom, A.W.M.Van den End en , ( ' ) '.800000000 • 1.800000000 o 8 . 000000800 8 ( 2) -4.713831'30 -I I. OS7Sg291 0 o 9.39Q297S30-1 "The Determination of the Orders of Process ( ;r) -2. 372639520 -2 8.888088880 • -~.$30S344S0 -I and Noise Dynamics", Automatica, vo1.l0, pp. ( 3. 2) IfATRIX . 245-256 (1974). ( I) ( Z) ( ' ) '.800009900 8 8 . 808000090 • [13] D.Q.Mayne, "Computational Procedure for (2) 1.809860D~0 8 -6 . 202103140 -I the Minimal Realization of Transfer Functions t ~> •. 880080900 • 1 ••00860800 8 Matrices", Proc. lEE, vo1.l55, pp.1363-l368 ( 2. ;r> IfATRIX C (1968). ( 1> ( -z> ( ;r) ( ' ) ~ .•2'8'86S0 -I 1.882080$70 8 6.188038730-1 [14] B.L.Ho, R.E . Kalman, "Effective Construc( 2) -4.113821150 -1 -J . BJB9~46S0 B S.999SsJsSD-1 tion of Linear State-Variable Models from - ( 2. 2) IfATRIX D Input/Output Functions", Regelungs technik, c I) ( Z) ( ' ) 9.982481820 -1 9 . 82181'420-1 vol.14, pp.545-548 (1966). ( 2 ) 1 .•84eJ678D • 1.819$7""0 • [15] T.Soderstrom, "Test of Pole-Zero CancelExample of computation by PIPAC-F
,J-
,.·····..··················· ..···0