- Email: [email protected]
Structural Identifiability of Nonlinear Systems: Algorithms Based on Differential Ideals
Copyright © IFAC System Identification, Copenhagen, Denmaric, 1994
STRUCTURAL IDENTIFIABILITY OF NONLINEAR SYSTEMS: ALGORITHMS BASED ON DIFFERENTIAL IDEALS.
L.d'Angio*, S.Audoly**, G.Bellu*, M.P. Saccomani*", C.Cobelli*** *Dip. di Matematica-Universita di Cagliari **Dip.di Ingegneria Strutturale-Universita di Cagliari *"Dip.di Elettronica e Infonnatica-Universita di Padova Abstract: The a priori structural identifiability problem for nonlinear dynamic systems is studied. Available methods for testing global identifiability via differential algebra are discussed. Two algorithms are proposed that utilise th.: infonnation provided by the characteristic set of the differential ideal associated to the system. Examples are presented. Key Words: Algebraic system theory, biomedical modelling. differential algebra, global identifiability, nonlinear systems, computer algebra.
been introduced originally by Ritt (1950) who also proposed, although not explicitly, an algorithm to get such a set. Ollivier (1990) and Glad and Ljung (1991), have recently shown how identifiability can be tested by analysing the characteristic set. They work in different differential rings, and the efficiency of their approaches is difficult to compare because Ollivier only describes a method to get a finite algebraic exhaustive summary, while Glad and Ljung suggest a complete procedure for testing identifiability, but do not provide algorithmic details for implementing it.
I. INTRODUCTION
A prion structural identifiability of dynamic nonlinear system models is an important problem in model identification since it is a prerequisite for well posedness of parameter estimation. While referring the reader to the literature (Cobelli et al. 1980; WaIter 1982; Godfrey et al. 1985; Vaida et al 1989) for the fundamentals on systems identifiability and a review of the methods to test it, the authors simply recall here that whatever method is chosen the problem boils down to find if a given polynomial map is injective, or, practically, if a set of algebraic simultaneous equations, the so-<:alled "exhaustive sununary " of the model, has one, more than one but finite, or an infinite number of solutions. For nonlinear systems, these equations are defined by the coefficients of the output series expansion and so they are not finite in number. Hence identifiability can only be tested in particular cases. The existence of a finite set of equations which contains all the information of the output series expansion, is guaranteed by a classical algebra theorem, (see e.g. D.Cox et al. 1992) however this theorem does not provide a method to derive such set of equations. Recently some novel approaches have been introduced in the identifiability area to overcome the difficulties of handling an infinite number of equations. They are based either on the "state isomorphism theorem" (Vaida and al. 1989), and on differential algebra (Glad and al 1994). The first has been applied to successfully testing a priori global identifiability of particular systems, but it is not easy to be automated in an algorithmic procedure. The methods based on differential algebra exploit the concept of characteristic set. The characteristic set of the differential ideal defined by the dynamic system equations, is a finite set of polynomials which sununarises all the information contained in the infinite differential ideal, so it is a finite exhaustive sununary of the model This concept has
The aim here is to develop algorithms, based on differential algebra, to test structural global identifiability of nonlinear dYl1amic models of relatively large dimension and reasonable complexity, typically those encountered in describing biological systems To do so algorithms have been constructed in house based on Ollivier's and Glad and Ljung's approaches. The procedures for reaching the characteristic set in the two different rings have been written in Reduce and state-of-the-art software has be used to implement Ritt's algoritlull. BOtll algoritluns seem heavy limited by computational burden: often they does not terminate when applied to nonlinear systems even of small dimensions (two or three state variables) and with relatively weak nonlinearities. In this paper a nunlber of strategies are proposed which, by taking into account the peculiar form of the dynamic system equations, allow a suitable extraction of the exhaustive sUOU11ary from the differential ideal defined by such equations. In tllis way the perfomlance of tlle algoritlmls improve particularly in the Ollivier's ring if applied to the models usually encowltered in biology and medicine.
The layout of the paper is tlle following: In section 2 tlle identifiability problem is briefly stated. 977
exhaustive summary of a model is a rational map such that:
In section 3 fWldarnentals on the differential ideals are recalled and Glad and LjWlg's and Ollivier's results are reported. In section 4 the two new algorithms are described. In section 5 examples together with a comparison between the two proposed algorithms are presented. Concluding remarks are in section 6.
R : C => Rk '<;f
The class of models of nonlinear dynamic systems is considered which can be described by:
R(p) = R(p·)
(h) (h)
Rand C will denote the set of real and complex numbers respectively. x E Rn is the state vector, u E R' and y E Rm the input and output vector respectively. p E P C RP the parameter vector, P the admissible parameter space, P its dimension A further equation:
3 DYNAMIC SYSTEMS AND DIFFERENTIAL IDEALS Consider the model defined by equations (1,,1 d 3), and assume that f,g,h are polynomial fWlctions, but the same also hold with slight changes in the rational case. For a formal description of basic objects utilised in differential algebra, the reader is referred to the specific literature (see e.g. Rill 1950, Kolkin 1973); here some definitions are only recalled. Let Z be the vector of the variables in Eqs. (1,,12.13). and K the set of their coefficients: the totality of polynomials in the variables z with coefficients on K is a Differential Polynomial Ring and will be denoted by K[Z]; K can be a differential ring itself.
F(p) = 0 can be added to Eqs. (I), describing known constraints on p; f,g,h,F are polynomial or rational fWlctions. In particular f,h and F describe the model structure and are thus dictated by a priori knowledge on the system, g is dictated by the input-output experiment designed to identify the model. Vectors x and p have usually a physical meaning, i.e. they represent e.g. masses and rate parameters respectively. Models belonging to the above class include the widely employed nonlinear compartmental models, in particular those exhibiting Michaelis-Menten and Languimur non linear kinetics. The system is supposed controllable and observable.
3.1. The characteristic set A ranking on the set of variables, must be first introduced, hence the derivatives are ranked according to a system that satisfies the following relations: u ~ 3 u u < v ==> 3 u < 3 v u and v are variables or derivatives of variables and 3 is an arbitrary derivation. The leader of a polynomial is the highest ranking derivative of that polynomial (it might be a derivative of order zero). The leader of a polynomial A well be denoted by llA. The initial of the polynomial A is the coefficient of the highest power of llA; the separant of A is 8A18llA: they will be denoted by I.. and SA respectively. If llA is lower than Us or llA = Us and deg(llA) < deg(Us), then the polynomial A is said to be of lower rank than B. The class of A is the greater p such that oz. is present in a term of A. If Ai is of class p > 0, Aj will be reduced with respect to Ai, if Aj is of lower rank than Ai in oz.. If a polynomial is not reduced with respect to another one, it can be reduced by. means of the pseudo division, that is: the polynomial of higher
2.1. Definitions Denote by B the map which associates u,y to each p. A model is called: structurally globally identifiable iff there is a set of measure zero V such that '<;f p E P-{V}, '<;fp' E CP
(2)
structurally locally identifiable iff there is a set of zero measure V such that '<;f pEP -{ V}, there is a neighbourhood U of p such that '<;f
(3,)
have one or more than one but finite, or an infinite number of solutions. The result will be true with probability one.
(1,)
B(p) = B(p') => P = p'.
(3)
If a model has a rational exhaustive sununary, the identifiability problem can be restated by utilising the map R in place of B. To check identifiability, one can choose p. in the parameter space P, and verify if equations
2 PROBLEM STATEMENT ANDDEFINmONS
x(t)=f1 x(t),u(t),p] x(O)=h(p) y(t)=g[ x(t),u(t),p]
p,p': R(p)=R(p') <=> B(p)=B(p').
p'E U, B(p) = B(p') ==> p--p'.
A system which is not locally neither globally identifiable is said nonidentifiable. A ratiol/al 978
rank is substituted by the rest of the division between the two polynomials, with respect to the leader of the first one. The system
A••A2•... Ar
4. TIIE NEW ALGORITI-lMS
The authors present here two new algoritJuns, denoted respectively CharsetJ and Charset2, which work in Ollivier's and Glad and LjWlg'S rings respectively. As is noted in the Introduction, experience shows that the success of an algoritJunic test of identifiability is heavily dependent on the ranking and reduction strategy. The new algoritJuns that will be described below in detail. have been constructed with particular care to this aspects. It should be noted that they are not general algoritJuns for getting the characteristic set of a generic differential ideal. but they are only aimed to reach a for testing suitable exhaustive sununary identifiability, of a dynamic system. Only a single or cumulated inputs are considered; in fact different inputs lead to different outputs and state vectors, so the related equations are completely independent one each other and the provided information can be finally assembled.
(4)
will be called a chain if either r = 1, and A. ~ 0 or r > 1, and for i > j Aj is of higher class than Ai and reduced with respect to Ai. The chain (4) will be said to be of higher rank of the chain BI.B2•.... B5 if. either: there is a j exceeding neither r nor s. such that Ai and Bi are of the same rank for i < j and that Aj is higher than Bj, or s > r and Ai and Bi are of the same rank for i ~ r. Let L be a finite or infinite set of differential polynomials. It is possible to form chains with differential polynorllials in I ; among all chains in L, some have lower rank. Each of such chain will be called a characteristic set. The peculiarity of the characteristic set is that of summarising all the information contained in the differential ideal into a finite number of differential polynomials.
4.1. AlgoritJun CharsetI The differential polynomials I, and I j are considered in the differential ring RI=R(p)[x,y,u,t) and a chain of minimum rank constructed by reducing thenl according to:
The search of the characteristic set is performed by a sequence of reductions between the polynomials of the given differential ideal until a particular set of minimum rank is reached. To do this the computational demand is heavily influenced by several factors such, the choice of the differential ring, the ranking of the indeterminates and the choice of the reduction sequence. The computational demand is heavily influenced by several factors such. the choice of the differential ring, the ranking of the indeterminates and the strategy in the reductions sequence.
u < y < x;
i < j <==> z{i) < z(j)
(5)
z(k) denote the k th derivative ofz. More in detail: the set of equations is ordered by increasing class, then reduction with respect to the x components is performed and a triangularised set on x , plus m equations free of x are obtained. Finally reduction with respect to y is made in these last equations, and a set of m equations is reached on u.y,p which is that part of the characteristic set that is informative on identifiability. In order to reduce the complexity, the process is divided into two steps with different strategies depending if the indeterminate x or y is involved in the reduction.
Ollivier (1990) suggests to analyse the characteristic set of Eqs. 11, and 13 in a ring where the only indeterminates are the output and input vectors components. Such ring will be denoted Ro = R(p)[y.u,t]. Once the characteristic set is known, the set of its coefficients on p is also known; these coefficients together with the initial conditions of the model, form an exhaustive summary The search is then performed with standard algebraic methods.
Reducing x. A matrix, Ix which describes the differential ideal under study is first associated to Eqs. 11 and h, with initial dimension (n+ 1)x(n+m). Its first n rows correspond to the n state components and the (n+ l)th to the derivative appearing in the related equation, while its columns correspond to m+n output and state equations. Because only one derivative appears in the original equations, in particular a first order derivative, only one row is sufficient to give information on it. The matrix is transformed step by step until a reduced form is reached in which the x components are eliminated. But during the reduction process, a number of derivatives of the original equations, that obviously belong to the same ideal, must be introduced, so the dimension of
Glad and Ljung (1991.1994) consider a differential ring in which all the u.y.x.P are indeterminates. Such ring will be denoted Ro = R[x.p.y.u,t]. They choose a ranking that produces a characteristic set whose simple observation allows to test identifiability. It has the form: C = A.(u.y)•........ Am(u,y). B.(u,y,P,), B2(U.Y.P•• J}2)•... Bp(u.y,pl, ... pp), C,(u.Y.P,XI), C2(U.Y.P.X2), .... Cn(u,y,p,Xn) B••B2•...Bp form a triangularised set on p from which the number of solutions can be derived. 979
the matrix could increase.
is a priori known, all the vector p components are different from zero, so the possible solutions pi = 0, have to be rejected. To do this a function res(o is created which performs the pseudodivision and automatical1y divides by the conunon factors leading to zero solutions; such division, not only eliminates a false solution, but significantly reduces the complexity.
The choice of the ranking and of the reduction sequence aims to avoid that, either derivatives of order higher than I appear, or, if i j, derivatives of xj are introduced into the ith state equation. In fact the observability of the system with the chosen ranking al10ws that derivatives of xi are eliminated by equations in which only derivatives of y appear. To this end the first order derivative of xi generated by derivation of the jth equation is directly eliminated by reducing it by the ith state equation. In this way during the reduction process, Ix summarises al1 the information on the ideal but its dimension does not increase.
*
An alternative strategy. If, despite of the strategies introduced for reducing complexity, tile above procedure fails, the problem is not computer time but the success of the test. In this case an alternative strategy is proposed, based on the fol1owing considerations:
Let Si the jth equation of S. If Si does not contain a monic monomial, the equation is suitably divided by a coefficient on p. so that a order zero term on p is created. Let denote gl,g2..gsi the resulting coefficients. The obtained equation, whose left hand side is a differential polynomial on y and its derivatives, is derived as many times as the nwnber of its coefficients on p minus one The si-I derivatives, together with Si foml a linear system of Si equations on the Si unknowns gk, k =I ,2 ..Si. that uniquely determine gk if the jacobian has ful1 rank; if the jacobian has rank si-di, di > 0, the system wil1 provide only si-di coefficients in terms of the remaining di ones; they will form that part of the exhaustive summary that can be extracted by the equation Si. So some coefficients of Si can be considered known and extracted without solving the system provided that the rank of the system matrix is calculated. Resorting to the initial point it is possible calculate a numerical value of y and its derivatives by evaluating the rank of the Jacobian in p. at time t·. In case of rank < Si the initial point will be changed in order to be sure p. does not belong to set of zero measure in which the jacobian has not full rank. The result will be true with probability one. If a similar strategy is applied to each equation of S a complete exhaustive summary is reached. The above procedure is longer than the previous one but it involves more simple equations, and so ends more easily.
The first stage of the algoritlun returns, in addition to a triangularised set on x, a set of m equations free of x. Let denote by S the set of these equations. The corresponding polynomials belong to the differential ring Ro, i.e. they are differential polynomials having y and u as indeterminates, and with coefficients in the algebraic ring R(p). If in this stage the number of different coefficients on p which appears in the equations set S is less than P, the system is nonidentifiable and the procedure stops. Reducing y. The set S is obtained by eliminating from the output and state equations the x components and their derivatives. In the above procedure at most n-m derivatives of each y component can be involved; hence a new matrix Iy is constructed which has initially m columns and m(n-m+2) rows at most. This matrix completely describing the y components and their derivatives that appear in the equations at this stage. Iy describes the reduction process until the characteristic set is reached. The exhaustive summary. The left hand side of the obtained m equations so obtained, can be considered as a characteristic set in the Ollivier's differential of the ring Ro. By dividing each element characteristic set by a suitable coefficient, a monic monomial is created; in this form the characteristic set uniquely provides its coefficients on p (see Ollivier )990). These coefficients can be extracted by a suitable procedure and form an exhaustive summary of the model. Identifiability can now be tested by (3) exploiting e.g. the Groebner basis. The choice of an initial point p. with numerical coordinates, and the use of (3,), further reduces the complexity in such a way the Groebner basis algoritlun works more easily. The result will be true with probability one. Trials with different initial points are necessary in order to make sure that p. does not belong to the set V.
4.2. Algorithm Charseh A second procedure has been implemented drawing inspiration from Glad's characteristic set, which once reached, tests identifiability directly, without the need of initial point and Groebner basis algorithm. Such procedure will be named Charse12. Consider the ring Ra with the fol1owing ranking:
u
j
< j <=> z{i) < z{j) "t i, j
The procedure Charseh coincides \vith Charset, until the set S of m equations free of x is reached.
On the pseudodivision. Because the model structure
980
equation in S together with its derivatives can be considered as a linear system on the gJt tmknown, that can be calculated by any technique for linear equations. Once the groups are calculated, the triangularisation can be more easily performed.
Charseu, according with the choice of the ring and ranking continues by reducing with respect to vector p, then y. Since p is a constant vector, the equations dp/dt = 0, have not been utilised as in the Glad and Ljung's procedure, having only the role of introducing the right number of equations for completing the test. If such equations are not reported, that is, one derives only y and u, the reduction process is faster, but because the ditTerential ideal has infinite elements, we are still faced with the problem of handling an infinite set of polynomials (see Beker, Weisspfening 1993). Thus a criterion is needed to decide how many derivatives of each equation of S it is necessary to get all the information available. To this aim the following consideration are utilised.
Some comments on the result. In addition to the improvement, because p, is not derived this method realises a number of advantages if compared to the simple reduction of equations with respect to p and y; in fact: the procedure stops with the right answer without performing such reduction either if the number of groups or the cardinality of T is less than P; if the system is identifiable, so the procedure does not stop, we avoid the further heavy reduction of T because only initials are reduced by the equations free of p. Moreover exploiting the linear equations on gk, the chance of final success further increases. The authors have not reached with such procedure, the complete characteristic set, but only that part of it which provides significant information for testing identifiability; no further equation on y and u is significant for this end; in fact any other derivative of S, would be of greater rank on y, than equations in T and so not able to reducing them. This algoritlun as the previous one has been constructed without taking into account the initial conditions and or any additional constraint; if they are knO\\ll some additional equations on the parameters must be added to the Charset, and Charseb issues to obtain the complete exhaustive
Reducing p. It can be seen from the Glad and Ljung's results that the polynomials that are informative on Identifiability, are the p polynomials Si i=1,2, ..p. Moreover other ID polynomials free of p, relating the vector y and u components, must be calculated, and the Si reduced with respect to them in order to be sure that no initial of Si vanishes. The vector p components are distributed in the equation belonging to S, into combinations multiplying polynomials on y and u and their derivatives. Thus if S is derived, being p a constant vector, such combinations does not change, so they are present in the same way in each derivative. It ensues that the number of derivatives of each equation which are significant for the procedure equates at most the number of such groups minus one, with the notations introduced in 4.1, Si-I. Hence the number of significant equations available for reducing p can not exceed S = L Si. If S < P, the system is nonidentifiable and the procedure stops; if S ~ P for each equation in S, Si - I derivatives are performed They will be sufficient, if the system is identifiable, to triangularise S with respect to p. Let T be the triangularised set on p components, and nT its cardinality, if nT < P, the system is nonidentifiable and the procedure stops.
summary.
5. EXAMPLES EXAMPLE I
Model from Godfrey et al.( 1991) The system is reduced to th.: integer foon : Am • DF(X1I ).nOXlI)· DF(X(\).T)"T3· XII)' "TI - XIIlOXI1'"T' - Xli )0U( I)m • XI I,"T' "T3 • XI' ,"T1 . X(1)"T3"T' -U(I)"T3m A(1)·· DF(XI1),n - XI I)"TI + X(1.'"T' + X(1'"T~ A()'- - Xli )"T6 • Y(I)
ALGORITHM Charset,
charsetl := _Y(I)' "TI"T~' YIIl' 0lJ( I)m"T6°IT' • n~ Y(I)' °DYlI,1) - (Y(I)' °OY",II)Om + T4. n, + YlI)' °00 1,I)m"T. • Y(I)' "T6°,-1"T1 "T3"T~ - n"T' - n"T~) + 1°Y(I)°lJ(I)"T3m"T6' °IT" n,-1°YII)oOYII.1)"T3"T6 - (1°Y(I)ODY(I,I)"T3"T6)'(T1 + T4. n) + 1°Y(I )°001,1)"T3""".' - (Y(I)"T3"T6' )'ITI "T3"T~ + n"T' + n"TSj + lJ(1)"T3' m-T6' 'IT' + T5) - OYlI.1)"T3' "T6' - (OY(I,I)"T3"T6' )'ITI"T3 + n + n"T' + n"T~) + OOI,I)"T3' m"T6'
On the initials. If nT = P, let be IT the vector of the initials of the equations in T; the elements of IT being polynomials on y, U and their derivatives. Some initial can vanish, if the related equation is reduced by the relations free of p, linking u,y, so it is necessary to verify if they are all ditTerent from zero. To do this the ID further derivatives of S of minimum rank on y, are added to the system, the vector p components eliminated by reducing them by S and the resulting equations utilised for reducing the initials with respect to y.
Coefficients are eXllacted and groebner Basis calulated. Let P'a (1(1,k2,k),k.,k,,Iu,k,.k,) Groebner Basis:
{T, - kl,Tz"k' - T'"k2, • T'"kl + T)"k7, T.·k., T,.Ju, T,"T, • k,"k,,) Time: 6700 ms; with a numerical initial point time 5980 ms. The system is not identifiable but if either T, or T, are known, it is globally identifiable. The algorilrn Char.;.:t2 has been also applied to this model, but the result (in symbolic foon). is practically unreadable because it involves a huge number oftenns.
The grouPS strategy. The reduction of S with respect to p, involves very hard symbolic calculations so it can be performed in two steps In fact, it ensues from the reasoning in 4.1 on the gJt, that each 981
dimension and reasonable complexity, like those usually encountered for describing biological systems and presented in the case studies, have been successfulIy tested. The algoritluns have been implemented in Reduce 3.4.1 on P.C. 486 33 MHz. Charset I appears more efficient than Charset2 which for state dimension three or greater often does not terminate. Charset I handles in acceptable computer time many models of the literature. At this stage it is difficult to define the domain of applicability of the new algoritluns. The authors can only claim that the largest nonlinear model successfully tested has n=4 and P=13, but, in their opinion these bounds can be relaxed by further refining the strategies for calculating a finite exhaustive sununary, and by the evolution of symbolic calculus.
EXAMPLE 2
Model from Chapcll et al. 1990 DF(X.~VMOXl(l(M+ X)-KOloX
Y-C°X
ALGORITHM Charset2
charset2:( _20DF(Y,T,4)ODF(Y,T,2)ODF(Y,1)' °Y' + 2"DF(Y,T,4)"DF(Y,1)" 0y. 3°DF(Y,T.3)' 0DF(Y,1)' oy' + 2°DF(Y,T.3)"DF(Y,T,2)' °DF(Y,:!I°y' - I4"DF(Y,T.3)ODF(Y,T,2)ODF(Y,1) 0y + 60DF(Y,T.3)ODF(Y,1)' - 300F(Y,T,2)" 0y' + 12°0F(Y,T,2)' °OF(Y,1)' 0y - 6°0F(Y,T,2)' °OF(Y,1)" , OF(Y,T.3)' °OFrY,1)' °C°y' °VM - 2°0F(Y,T.3)OOF(Y,T,2)' °OF(Y,1)°C°y' °VM _4°0F(Y,T.3)OOF(Y,T,2\OOF(Y,1)' °coY' °VM • 4°0F(Y,T.3)OOF(Y,1)f °C°y' °VM + OF(Y,T,2)" °C°y' °VM - 400F(Y,T,2)J 00F(Y,1)' °Y' + 400F(Y,T,2)' °DF(Y,1)' °coy J °VM + 12"OF(Y,T,2)'OOF(Y,1)' 0y' - 1200F(Y,T,2)OOF(Y,1)' 0y · 8°0F(Y,T,2rDF(Y,1)' oCOYOVM • 4°0F(Y,1)' • 4°0F(Y,1) °coVM, - DF(Y,T.3)ODF(Y,1)' 0y + OF(Y,T.3)OOF(Y,1)°y' °KOI + 3°0F(Y,T,2)' °OF(Y,1)"Y . OF(Y,T,2)' 0y' °KO) .2°0F(Y,T,2}"OF(Y,1)' - 2°0FrY,T,2)OOF(Y,1)' oYOKOI • 2°0F(Y,1) °KOI, DF(Y,T.3)' °OF(Y,1)' °Y' °VM - 2°0F(Y,T.3\OOF(Y,T,2)' °OHY,1)°Y' °VM · 200F(Y,T.3)OOF(Y,T,2\OOF\Y,Tj' 0y' 0VM • 2°0F(Y,T.3\OOF(Y,1)f °Y' °VM • OHY,T,2)" °Y' °VM J + 4°0F(Y,T,2)' °OF(Y,1)' °V °KM + 2°0F(Y,T,2)' °OF(Y,1)' 0y' °VM - 12°0F(Y,T,2\' °OF(Y,1)' °v' °KM - 2°0F(Y,T,2/o0F(Y,1)" °v' °VM + 12°0F(Y,T,2)OOF(Y,1)' "Y"KM· 4°0FrY,1)' °KM)
7. REFERENCES Becker, T, and V.Weispfelming (1993). Groebner Basis, Springer-Verlag, New York. Chappel, M.J., and K.Goffrey (1991). Structural Identifiability of non Linear Systems:Application 9th IFAC-IFORS to a Batch Reactor. SymPOSiunl- Cs.Banyasz and L.Keviczky editors Budapest - 8 July 1991. Cobelli, C., and J.1. DiStefano ill (1980). Parameter and structural identifiability concepts an anlbiguities; A critical review and analysis, AM.1.Physiol.VoI239:R7 - R24. Cox, D., 1. Little, and D.O·Shea (1992) Ideals varieties and algoritluns Springer-Verlag. Glad, T and L.Liung (1991). Testing Global Identifiability for Arbitrary Model Parametrisation 9th IFAC-IFORS SymposiwnCs.Banyasz and L.Keviczky editors. Glad, T, and L.Lijung (1994). On Global Identifiability for arbitrary model paranletrisation. Automatica Vo1.30 n.2 pp 265276. Elsevier Science. DiStefano (1985). Godfrey, K., and 1.1. Identifiability of model parameters E.Walter editor Pergamon Press. Kolchin, E.R., (1973). Differential Algebra and Algebraic Groups, Acadenlic Press, New York. Ollivier, F. (1990). Le problem de l'identifiabilite'... , These de doctor at en Sciences 1990, Laboratoir des signaux et systems. Palaiseau - Paris. Ritt, H., (1950). Differential Algebra. American Mathematical Society. Providence, RI. Vaida, S., K.R. Godfrey, and H.Rabitz (1989). Similarity transformation approach to structural identifiability of nonlinear models. Math. Biosci. 93:217-248. WaIter, E. (1982). Identifiability of State Space Models Springer-Verlag. New York.
Time: 2690 ms. The system is not identifiable but parameter KO J is globally identifiable. The algoritJun Charsetl has been applied to this example; it requires 440 ms; with a numerical initial point 330 ms. In the following two examples the characteristic sets are too large to be reponed; hence only the results are summarised. EXAMPLE 3
Nonlinear comp. model from Waiter 1982 pg. 167 ex.7.4. n = 4, P= 13 (unkno"n volume v,), input in comp. I, outputs in compartments I and 2. ALGORITIIM Charset I The system is globally identifiable. With a numerical initial point, time: 42735 ms. EXAMPLE 4
Linear comp. model from Ollivier 1990 pg 110 (synthesis of anunonia. mod of Horiuti). n= 5, P = 18 input in comp. 5 and outputs in comp. J and 5. ALGORITIIM Charset 1 Time for Characteristic set 123961 ms. Before calculating the Groebner basis, constraints on p due to the system steady state (like in Ollivier) have been introduced. while only the directly measurable parameters have been considered known. The remaining 9 parameters are globally identifiable. With a numerical initial point, Groebner basis time: 90885 ms.
6. CONCLUSIONS Two new algoritluns, based on differential algebra, for testing structural global identifiabi:ity of nonlinear dynamic models have been presented. The algoritluns Charset I and Charset2, which work in OlIivier's and Glad and Ljung's rings respectively, employ some novel strategies in the ranking and reduction sequence Models of relatively large
982
STRUCTURAL IDENTIFIABILITY OF NONLINEAR SYSTEMS: ALGORITHMS BASED ON DIFFERENTIAL IDEALS.
L.d'Angio*, S.Audoly**, G.Bellu*, M.P. Saccomani*", C.Cobelli*** *Dip. di Matematica-Universita di Cagliari **Dip.di Ingegneria Strutturale-Universita di Cagliari *"Dip.di Elettronica e Infonnatica-Universita di Padova Abstract: The a priori structural identifiability problem for nonlinear dynamic systems is studied. Available methods for testing global identifiability via differential algebra are discussed. Two algorithms are proposed that utilise th.: infonnation provided by the characteristic set of the differential ideal associated to the system. Examples are presented. Key Words: Algebraic system theory, biomedical modelling. differential algebra, global identifiability, nonlinear systems, computer algebra.
been introduced originally by Ritt (1950) who also proposed, although not explicitly, an algorithm to get such a set. Ollivier (1990) and Glad and Ljung (1991), have recently shown how identifiability can be tested by analysing the characteristic set. They work in different differential rings, and the efficiency of their approaches is difficult to compare because Ollivier only describes a method to get a finite algebraic exhaustive summary, while Glad and Ljung suggest a complete procedure for testing identifiability, but do not provide algorithmic details for implementing it.
I. INTRODUCTION
A prion structural identifiability of dynamic nonlinear system models is an important problem in model identification since it is a prerequisite for well posedness of parameter estimation. While referring the reader to the literature (Cobelli et al. 1980; WaIter 1982; Godfrey et al. 1985; Vaida et al 1989) for the fundamentals on systems identifiability and a review of the methods to test it, the authors simply recall here that whatever method is chosen the problem boils down to find if a given polynomial map is injective, or, practically, if a set of algebraic simultaneous equations, the so-<:alled "exhaustive sununary " of the model, has one, more than one but finite, or an infinite number of solutions. For nonlinear systems, these equations are defined by the coefficients of the output series expansion and so they are not finite in number. Hence identifiability can only be tested in particular cases. The existence of a finite set of equations which contains all the information of the output series expansion, is guaranteed by a classical algebra theorem, (see e.g. D.Cox et al. 1992) however this theorem does not provide a method to derive such set of equations. Recently some novel approaches have been introduced in the identifiability area to overcome the difficulties of handling an infinite number of equations. They are based either on the "state isomorphism theorem" (Vaida and al. 1989), and on differential algebra (Glad and al 1994). The first has been applied to successfully testing a priori global identifiability of particular systems, but it is not easy to be automated in an algorithmic procedure. The methods based on differential algebra exploit the concept of characteristic set. The characteristic set of the differential ideal defined by the dynamic system equations, is a finite set of polynomials which sununarises all the information contained in the infinite differential ideal, so it is a finite exhaustive sununary of the model This concept has
The aim here is to develop algorithms, based on differential algebra, to test structural global identifiability of nonlinear dYl1amic models of relatively large dimension and reasonable complexity, typically those encountered in describing biological systems To do so algorithms have been constructed in house based on Ollivier's and Glad and Ljung's approaches. The procedures for reaching the characteristic set in the two different rings have been written in Reduce and state-of-the-art software has be used to implement Ritt's algoritlull. BOtll algoritluns seem heavy limited by computational burden: often they does not terminate when applied to nonlinear systems even of small dimensions (two or three state variables) and with relatively weak nonlinearities. In this paper a nunlber of strategies are proposed which, by taking into account the peculiar form of the dynamic system equations, allow a suitable extraction of the exhaustive sUOU11ary from the differential ideal defined by such equations. In tllis way the perfomlance of tlle algoritlmls improve particularly in the Ollivier's ring if applied to the models usually encowltered in biology and medicine.
The layout of the paper is tlle following: In section 2 tlle identifiability problem is briefly stated. 977
exhaustive summary of a model is a rational map such that:
In section 3 fWldarnentals on the differential ideals are recalled and Glad and LjWlg's and Ollivier's results are reported. In section 4 the two new algorithms are described. In section 5 examples together with a comparison between the two proposed algorithms are presented. Concluding remarks are in section 6.
R : C => Rk '<;f
The class of models of nonlinear dynamic systems is considered which can be described by:
R(p) = R(p·)
(h) (h)
Rand C will denote the set of real and complex numbers respectively. x E Rn is the state vector, u E R' and y E Rm the input and output vector respectively. p E P C RP the parameter vector, P the admissible parameter space, P its dimension A further equation:
3 DYNAMIC SYSTEMS AND DIFFERENTIAL IDEALS Consider the model defined by equations (1,,1 d 3), and assume that f,g,h are polynomial fWlctions, but the same also hold with slight changes in the rational case. For a formal description of basic objects utilised in differential algebra, the reader is referred to the specific literature (see e.g. Rill 1950, Kolkin 1973); here some definitions are only recalled. Let Z be the vector of the variables in Eqs. (1,,12.13). and K the set of their coefficients: the totality of polynomials in the variables z with coefficients on K is a Differential Polynomial Ring and will be denoted by K[Z]; K can be a differential ring itself.
F(p) = 0 can be added to Eqs. (I), describing known constraints on p; f,g,h,F are polynomial or rational fWlctions. In particular f,h and F describe the model structure and are thus dictated by a priori knowledge on the system, g is dictated by the input-output experiment designed to identify the model. Vectors x and p have usually a physical meaning, i.e. they represent e.g. masses and rate parameters respectively. Models belonging to the above class include the widely employed nonlinear compartmental models, in particular those exhibiting Michaelis-Menten and Languimur non linear kinetics. The system is supposed controllable and observable.
3.1. The characteristic set A ranking on the set of variables, must be first introduced, hence the derivatives are ranked according to a system that satisfies the following relations: u ~ 3 u u < v ==> 3 u < 3 v u and v are variables or derivatives of variables and 3 is an arbitrary derivation. The leader of a polynomial is the highest ranking derivative of that polynomial (it might be a derivative of order zero). The leader of a polynomial A well be denoted by llA. The initial of the polynomial A is the coefficient of the highest power of llA; the separant of A is 8A18llA: they will be denoted by I.. and SA respectively. If llA is lower than Us or llA = Us and deg(llA) < deg(Us), then the polynomial A is said to be of lower rank than B. The class of A is the greater p such that oz. is present in a term of A. If Ai is of class p > 0, Aj will be reduced with respect to Ai, if Aj is of lower rank than Ai in oz.. If a polynomial is not reduced with respect to another one, it can be reduced by. means of the pseudo division, that is: the polynomial of higher
2.1. Definitions Denote by B the map which associates u,y to each p. A model is called: structurally globally identifiable iff there is a set of measure zero V such that '<;f p E P-{V}, '<;fp' E CP
(2)
structurally locally identifiable iff there is a set of zero measure V such that '<;f pEP -{ V}, there is a neighbourhood U of p such that '<;f
(3,)
have one or more than one but finite, or an infinite number of solutions. The result will be true with probability one.
(1,)
B(p) = B(p') => P = p'.
(3)
If a model has a rational exhaustive sununary, the identifiability problem can be restated by utilising the map R in place of B. To check identifiability, one can choose p. in the parameter space P, and verify if equations
2 PROBLEM STATEMENT ANDDEFINmONS
x(t)=f1 x(t),u(t),p] x(O)=h(p) y(t)=g[ x(t),u(t),p]
p,p': R(p)=R(p') <=> B(p)=B(p').
p'E U, B(p) = B(p') ==> p--p'.
A system which is not locally neither globally identifiable is said nonidentifiable. A ratiol/al 978
rank is substituted by the rest of the division between the two polynomials, with respect to the leader of the first one. The system
A••A2•... Ar
4. TIIE NEW ALGORITI-lMS
The authors present here two new algoritJuns, denoted respectively CharsetJ and Charset2, which work in Ollivier's and Glad and LjWlg'S rings respectively. As is noted in the Introduction, experience shows that the success of an algoritJunic test of identifiability is heavily dependent on the ranking and reduction strategy. The new algoritJuns that will be described below in detail. have been constructed with particular care to this aspects. It should be noted that they are not general algoritJuns for getting the characteristic set of a generic differential ideal. but they are only aimed to reach a for testing suitable exhaustive sununary identifiability, of a dynamic system. Only a single or cumulated inputs are considered; in fact different inputs lead to different outputs and state vectors, so the related equations are completely independent one each other and the provided information can be finally assembled.
(4)
will be called a chain if either r = 1, and A. ~ 0 or r > 1, and for i > j Aj is of higher class than Ai and reduced with respect to Ai. The chain (4) will be said to be of higher rank of the chain BI.B2•.... B5 if. either: there is a j exceeding neither r nor s. such that Ai and Bi are of the same rank for i < j and that Aj is higher than Bj, or s > r and Ai and Bi are of the same rank for i ~ r. Let L be a finite or infinite set of differential polynomials. It is possible to form chains with differential polynorllials in I ; among all chains in L, some have lower rank. Each of such chain will be called a characteristic set. The peculiarity of the characteristic set is that of summarising all the information contained in the differential ideal into a finite number of differential polynomials.
4.1. AlgoritJun CharsetI The differential polynomials I, and I j are considered in the differential ring RI=R(p)[x,y,u,t) and a chain of minimum rank constructed by reducing thenl according to:
The search of the characteristic set is performed by a sequence of reductions between the polynomials of the given differential ideal until a particular set of minimum rank is reached. To do this the computational demand is heavily influenced by several factors such, the choice of the differential ring, the ranking of the indeterminates and the choice of the reduction sequence. The computational demand is heavily influenced by several factors such. the choice of the differential ring, the ranking of the indeterminates and the strategy in the reductions sequence.
u < y < x;
i < j <==> z{i) < z(j)
(5)
z(k) denote the k th derivative ofz. More in detail: the set of equations is ordered by increasing class, then reduction with respect to the x components is performed and a triangularised set on x , plus m equations free of x are obtained. Finally reduction with respect to y is made in these last equations, and a set of m equations is reached on u.y,p which is that part of the characteristic set that is informative on identifiability. In order to reduce the complexity, the process is divided into two steps with different strategies depending if the indeterminate x or y is involved in the reduction.
Ollivier (1990) suggests to analyse the characteristic set of Eqs. 11, and 13 in a ring where the only indeterminates are the output and input vectors components. Such ring will be denoted Ro = R(p)[y.u,t]. Once the characteristic set is known, the set of its coefficients on p is also known; these coefficients together with the initial conditions of the model, form an exhaustive summary The search is then performed with standard algebraic methods.
Reducing x. A matrix, Ix which describes the differential ideal under study is first associated to Eqs. 11 and h, with initial dimension (n+ 1)x(n+m). Its first n rows correspond to the n state components and the (n+ l)th to the derivative appearing in the related equation, while its columns correspond to m+n output and state equations. Because only one derivative appears in the original equations, in particular a first order derivative, only one row is sufficient to give information on it. The matrix is transformed step by step until a reduced form is reached in which the x components are eliminated. But during the reduction process, a number of derivatives of the original equations, that obviously belong to the same ideal, must be introduced, so the dimension of
Glad and Ljung (1991.1994) consider a differential ring in which all the u.y.x.P are indeterminates. Such ring will be denoted Ro = R[x.p.y.u,t]. They choose a ranking that produces a characteristic set whose simple observation allows to test identifiability. It has the form: C = A.(u.y)•........ Am(u,y). B.(u,y,P,), B2(U.Y.P•• J}2)•... Bp(u.y,pl, ... pp), C,(u.Y.P,XI), C2(U.Y.P.X2), .... Cn(u,y,p,Xn) B••B2•...Bp form a triangularised set on p from which the number of solutions can be derived. 979
the matrix could increase.
is a priori known, all the vector p components are different from zero, so the possible solutions pi = 0, have to be rejected. To do this a function res(o is created which performs the pseudodivision and automatical1y divides by the conunon factors leading to zero solutions; such division, not only eliminates a false solution, but significantly reduces the complexity.
The choice of the ranking and of the reduction sequence aims to avoid that, either derivatives of order higher than I appear, or, if i j, derivatives of xj are introduced into the ith state equation. In fact the observability of the system with the chosen ranking al10ws that derivatives of xi are eliminated by equations in which only derivatives of y appear. To this end the first order derivative of xi generated by derivation of the jth equation is directly eliminated by reducing it by the ith state equation. In this way during the reduction process, Ix summarises al1 the information on the ideal but its dimension does not increase.
*
An alternative strategy. If, despite of the strategies introduced for reducing complexity, tile above procedure fails, the problem is not computer time but the success of the test. In this case an alternative strategy is proposed, based on the fol1owing considerations:
Let Si the jth equation of S. If Si does not contain a monic monomial, the equation is suitably divided by a coefficient on p. so that a order zero term on p is created. Let denote gl,g2..gsi the resulting coefficients. The obtained equation, whose left hand side is a differential polynomial on y and its derivatives, is derived as many times as the nwnber of its coefficients on p minus one The si-I derivatives, together with Si foml a linear system of Si equations on the Si unknowns gk, k =I ,2 ..Si. that uniquely determine gk if the jacobian has ful1 rank; if the jacobian has rank si-di, di > 0, the system wil1 provide only si-di coefficients in terms of the remaining di ones; they will form that part of the exhaustive summary that can be extracted by the equation Si. So some coefficients of Si can be considered known and extracted without solving the system provided that the rank of the system matrix is calculated. Resorting to the initial point it is possible calculate a numerical value of y and its derivatives by evaluating the rank of the Jacobian in p. at time t·. In case of rank < Si the initial point will be changed in order to be sure p. does not belong to set of zero measure in which the jacobian has not full rank. The result will be true with probability one. If a similar strategy is applied to each equation of S a complete exhaustive summary is reached. The above procedure is longer than the previous one but it involves more simple equations, and so ends more easily.
The first stage of the algoritlun returns, in addition to a triangularised set on x, a set of m equations free of x. Let denote by S the set of these equations. The corresponding polynomials belong to the differential ring Ro, i.e. they are differential polynomials having y and u as indeterminates, and with coefficients in the algebraic ring R(p). If in this stage the number of different coefficients on p which appears in the equations set S is less than P, the system is nonidentifiable and the procedure stops. Reducing y. The set S is obtained by eliminating from the output and state equations the x components and their derivatives. In the above procedure at most n-m derivatives of each y component can be involved; hence a new matrix Iy is constructed which has initially m columns and m(n-m+2) rows at most. This matrix completely describing the y components and their derivatives that appear in the equations at this stage. Iy describes the reduction process until the characteristic set is reached. The exhaustive summary. The left hand side of the obtained m equations so obtained, can be considered as a characteristic set in the Ollivier's differential of the ring Ro. By dividing each element characteristic set by a suitable coefficient, a monic monomial is created; in this form the characteristic set uniquely provides its coefficients on p (see Ollivier )990). These coefficients can be extracted by a suitable procedure and form an exhaustive summary of the model. Identifiability can now be tested by (3) exploiting e.g. the Groebner basis. The choice of an initial point p. with numerical coordinates, and the use of (3,), further reduces the complexity in such a way the Groebner basis algoritlun works more easily. The result will be true with probability one. Trials with different initial points are necessary in order to make sure that p. does not belong to the set V.
4.2. Algorithm Charseh A second procedure has been implemented drawing inspiration from Glad's characteristic set, which once reached, tests identifiability directly, without the need of initial point and Groebner basis algorithm. Such procedure will be named Charse12. Consider the ring Ra with the fol1owing ranking:
u
j
< j <=> z{i) < z{j) "t i, j
The procedure Charseh coincides \vith Charset, until the set S of m equations free of x is reached.
On the pseudodivision. Because the model structure
980
equation in S together with its derivatives can be considered as a linear system on the gJt tmknown, that can be calculated by any technique for linear equations. Once the groups are calculated, the triangularisation can be more easily performed.
Charseu, according with the choice of the ring and ranking continues by reducing with respect to vector p, then y. Since p is a constant vector, the equations dp/dt = 0, have not been utilised as in the Glad and Ljung's procedure, having only the role of introducing the right number of equations for completing the test. If such equations are not reported, that is, one derives only y and u, the reduction process is faster, but because the ditTerential ideal has infinite elements, we are still faced with the problem of handling an infinite set of polynomials (see Beker, Weisspfening 1993). Thus a criterion is needed to decide how many derivatives of each equation of S it is necessary to get all the information available. To this aim the following consideration are utilised.
Some comments on the result. In addition to the improvement, because p, is not derived this method realises a number of advantages if compared to the simple reduction of equations with respect to p and y; in fact: the procedure stops with the right answer without performing such reduction either if the number of groups or the cardinality of T is less than P; if the system is identifiable, so the procedure does not stop, we avoid the further heavy reduction of T because only initials are reduced by the equations free of p. Moreover exploiting the linear equations on gk, the chance of final success further increases. The authors have not reached with such procedure, the complete characteristic set, but only that part of it which provides significant information for testing identifiability; no further equation on y and u is significant for this end; in fact any other derivative of S, would be of greater rank on y, than equations in T and so not able to reducing them. This algoritlun as the previous one has been constructed without taking into account the initial conditions and or any additional constraint; if they are knO\\ll some additional equations on the parameters must be added to the Charset, and Charseb issues to obtain the complete exhaustive
Reducing p. It can be seen from the Glad and Ljung's results that the polynomials that are informative on Identifiability, are the p polynomials Si i=1,2, ..p. Moreover other ID polynomials free of p, relating the vector y and u components, must be calculated, and the Si reduced with respect to them in order to be sure that no initial of Si vanishes. The vector p components are distributed in the equation belonging to S, into combinations multiplying polynomials on y and u and their derivatives. Thus if S is derived, being p a constant vector, such combinations does not change, so they are present in the same way in each derivative. It ensues that the number of derivatives of each equation which are significant for the procedure equates at most the number of such groups minus one, with the notations introduced in 4.1, Si-I. Hence the number of significant equations available for reducing p can not exceed S = L Si. If S < P, the system is nonidentifiable and the procedure stops; if S ~ P for each equation in S, Si - I derivatives are performed They will be sufficient, if the system is identifiable, to triangularise S with respect to p. Let T be the triangularised set on p components, and nT its cardinality, if nT < P, the system is nonidentifiable and the procedure stops.
summary.
5. EXAMPLES EXAMPLE I
Model from Godfrey et al.( 1991) The system is reduced to th.: integer foon : Am • DF(X1I ).nOXlI)· DF(X(\).T)"T3· XII)' "TI - XIIlOXI1'"T' - Xli )0U( I)m • XI I,"T' "T3 • XI' ,"T1 . X(1)"T3"T' -U(I)"T3m A(1)·· DF(XI1),n - XI I)"TI + X(1.'"T' + X(1'"T~ A()'- - Xli )"T6 • Y(I)
ALGORITHM Charset,
charsetl := _Y(I)' "TI"T~' YIIl' 0lJ( I)m"T6°IT' • n~ Y(I)' °DYlI,1) - (Y(I)' °OY",II)Om + T4. n, + YlI)' °00 1,I)m"T. • Y(I)' "T6°,-1"T1 "T3"T~ - n"T' - n"T~) + 1°Y(I)°lJ(I)"T3m"T6' °IT" n,-1°YII)oOYII.1)"T3"T6 - (1°Y(I)ODY(I,I)"T3"T6)'(T1 + T4. n) + 1°Y(I )°001,1)"T3""".' - (Y(I)"T3"T6' )'ITI "T3"T~ + n"T' + n"TSj + lJ(1)"T3' m-T6' 'IT' + T5) - OYlI.1)"T3' "T6' - (OY(I,I)"T3"T6' )'ITI"T3 + n + n"T' + n"T~) + OOI,I)"T3' m"T6'
On the initials. If nT = P, let be IT the vector of the initials of the equations in T; the elements of IT being polynomials on y, U and their derivatives. Some initial can vanish, if the related equation is reduced by the relations free of p, linking u,y, so it is necessary to verify if they are all ditTerent from zero. To do this the ID further derivatives of S of minimum rank on y, are added to the system, the vector p components eliminated by reducing them by S and the resulting equations utilised for reducing the initials with respect to y.
Coefficients are eXllacted and groebner Basis calulated. Let P'a (1(1,k2,k),k.,k,,Iu,k,.k,) Groebner Basis:
{T, - kl,Tz"k' - T'"k2, • T'"kl + T)"k7, T.·k., T,.Ju, T,"T, • k,"k,,) Time: 6700 ms; with a numerical initial point time 5980 ms. The system is not identifiable but if either T, or T, are known, it is globally identifiable. The algorilrn Char.;.:t2 has been also applied to this model, but the result (in symbolic foon). is practically unreadable because it involves a huge number oftenns.
The grouPS strategy. The reduction of S with respect to p, involves very hard symbolic calculations so it can be performed in two steps In fact, it ensues from the reasoning in 4.1 on the gJt, that each 981
dimension and reasonable complexity, like those usually encountered for describing biological systems and presented in the case studies, have been successfulIy tested. The algoritluns have been implemented in Reduce 3.4.1 on P.C. 486 33 MHz. Charset I appears more efficient than Charset2 which for state dimension three or greater often does not terminate. Charset I handles in acceptable computer time many models of the literature. At this stage it is difficult to define the domain of applicability of the new algoritluns. The authors can only claim that the largest nonlinear model successfully tested has n=4 and P=13, but, in their opinion these bounds can be relaxed by further refining the strategies for calculating a finite exhaustive sununary, and by the evolution of symbolic calculus.
EXAMPLE 2
Model from Chapcll et al. 1990 DF(X.~VMOXl(l(M+ X)-KOloX
Y-C°X
ALGORITHM Charset2
charset2:( _20DF(Y,T,4)ODF(Y,T,2)ODF(Y,1)' °Y' + 2"DF(Y,T,4)"DF(Y,1)" 0y. 3°DF(Y,T.3)' 0DF(Y,1)' oy' + 2°DF(Y,T.3)"DF(Y,T,2)' °DF(Y,:!I°y' - I4"DF(Y,T.3)ODF(Y,T,2)ODF(Y,1) 0y + 60DF(Y,T.3)ODF(Y,1)' - 300F(Y,T,2)" 0y' + 12°0F(Y,T,2)' °OF(Y,1)' 0y - 6°0F(Y,T,2)' °OF(Y,1)" , OF(Y,T.3)' °OFrY,1)' °C°y' °VM - 2°0F(Y,T.3)OOF(Y,T,2)' °OF(Y,1)°C°y' °VM _4°0F(Y,T.3)OOF(Y,T,2\OOF(Y,1)' °coY' °VM • 4°0F(Y,T.3)OOF(Y,1)f °C°y' °VM + OF(Y,T,2)" °C°y' °VM - 400F(Y,T,2)J 00F(Y,1)' °Y' + 400F(Y,T,2)' °DF(Y,1)' °coy J °VM + 12"OF(Y,T,2)'OOF(Y,1)' 0y' - 1200F(Y,T,2)OOF(Y,1)' 0y · 8°0F(Y,T,2rDF(Y,1)' oCOYOVM • 4°0F(Y,1)' • 4°0F(Y,1) °coVM, - DF(Y,T.3)ODF(Y,1)' 0y + OF(Y,T.3)OOF(Y,1)°y' °KOI + 3°0F(Y,T,2)' °OF(Y,1)"Y . OF(Y,T,2)' 0y' °KO) .2°0F(Y,T,2}"OF(Y,1)' - 2°0FrY,T,2)OOF(Y,1)' oYOKOI • 2°0F(Y,1) °KOI, DF(Y,T.3)' °OF(Y,1)' °Y' °VM - 2°0F(Y,T.3\OOF(Y,T,2)' °OHY,1)°Y' °VM · 200F(Y,T.3)OOF(Y,T,2\OOF\Y,Tj' 0y' 0VM • 2°0F(Y,T.3\OOF(Y,1)f °Y' °VM • OHY,T,2)" °Y' °VM J + 4°0F(Y,T,2)' °OF(Y,1)' °V °KM + 2°0F(Y,T,2)' °OF(Y,1)' 0y' °VM - 12°0F(Y,T,2\' °OF(Y,1)' °v' °KM - 2°0F(Y,T,2/o0F(Y,1)" °v' °VM + 12°0F(Y,T,2)OOF(Y,1)' "Y"KM· 4°0FrY,1)' °KM)
7. REFERENCES Becker, T, and V.Weispfelming (1993). Groebner Basis, Springer-Verlag, New York. Chappel, M.J., and K.Goffrey (1991). Structural Identifiability of non Linear Systems:Application 9th IFAC-IFORS to a Batch Reactor. SymPOSiunl- Cs.Banyasz and L.Keviczky editors Budapest - 8 July 1991. Cobelli, C., and J.1. DiStefano ill (1980). Parameter and structural identifiability concepts an anlbiguities; A critical review and analysis, AM.1.Physiol.VoI239:R7 - R24. Cox, D., 1. Little, and D.O·Shea (1992) Ideals varieties and algoritluns Springer-Verlag. Glad, T and L.Liung (1991). Testing Global Identifiability for Arbitrary Model Parametrisation 9th IFAC-IFORS SymposiwnCs.Banyasz and L.Keviczky editors. Glad, T, and L.Lijung (1994). On Global Identifiability for arbitrary model paranletrisation. Automatica Vo1.30 n.2 pp 265276. Elsevier Science. DiStefano (1985). Godfrey, K., and 1.1. Identifiability of model parameters E.Walter editor Pergamon Press. Kolchin, E.R., (1973). Differential Algebra and Algebraic Groups, Acadenlic Press, New York. Ollivier, F. (1990). Le problem de l'identifiabilite'... , These de doctor at en Sciences 1990, Laboratoir des signaux et systems. Palaiseau - Paris. Ritt, H., (1950). Differential Algebra. American Mathematical Society. Providence, RI. Vaida, S., K.R. Godfrey, and H.Rabitz (1989). Similarity transformation approach to structural identifiability of nonlinear models. Math. Biosci. 93:217-248. WaIter, E. (1982). Identifiability of State Space Models Springer-Verlag. New York.
Time: 2690 ms. The system is not identifiable but parameter KO J is globally identifiable. The algoritJun Charsetl has been applied to this example; it requires 440 ms; with a numerical initial point 330 ms. In the following two examples the characteristic sets are too large to be reponed; hence only the results are summarised. EXAMPLE 3
Nonlinear comp. model from Waiter 1982 pg. 167 ex.7.4. n = 4, P= 13 (unkno"n volume v,), input in comp. I, outputs in compartments I and 2. ALGORITIIM Charset I The system is globally identifiable. With a numerical initial point, time: 42735 ms. EXAMPLE 4
Linear comp. model from Ollivier 1990 pg 110 (synthesis of anunonia. mod of Horiuti). n= 5, P = 18 input in comp. 5 and outputs in comp. J and 5. ALGORITIIM Charset 1 Time for Characteristic set 123961 ms. Before calculating the Groebner basis, constraints on p due to the system steady state (like in Ollivier) have been introduced. while only the directly measurable parameters have been considered known. The remaining 9 parameters are globally identifiable. With a numerical initial point, Groebner basis time: 90885 ms.
6. CONCLUSIONS Two new algoritluns, based on differential algebra, for testing structural global identifiabi:ity of nonlinear dynamic models have been presented. The algoritluns Charset I and Charset2, which work in OlIivier's and Glad and Ljung's rings respectively, employ some novel strategies in the ranking and reduction sequence Models of relatively large
982