Ind. Eng. Chem. Res. 2004, 43, 1251-1259
1251
Identifiability of Linear Time-Invariant Differential-Algebraic Systems. 2. The Differential-Algebraic Approach Amos Ben-Zvi, P. James McLellan,* and Kim B. McAuley Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6
A mathematical model is identifiable if and only if there is a unique relationship between each parameter value and the input-output behavior of the model. If a model is not identifiable, there is no unique solution to the parameter estimation problem, regardless of the number and type of experiments that are performed. A method for testing the identifiability of linear timeinvariant (LTI) differential-algebraic equation (DAE) systems, based on differential algebra, is presented. In the proposed approach, the LTI DAE system is treated as a set of linear mappings in the input, output, and state variables. The proposed treatment allows the input-output representation of the system to be obtained by combining and differentiating elements of this set. The identifiability of the system is tested by checking whether the relationship between the model parameters and the coefficients in the input-output representation of the system is one-to-one. One benefit of the proposed method is that it readily produces a simplified realization of the system that is identifiable even when the original LTI DAE model is not identifiable. Necessary and sufficient conditions for local and global identifiability are presented, and the application of the proposed method is illustrated using a simplified gas-phase reaction model. 1. Introduction Differential-algebraic equation systems (DAEs), also called implicit or semistate systems, are common in the study of mechanical and electrical systems and chemical processes.1-5 Accurate estimation of model parameters is important for the design, scale-up, and optimization of these systems.3,6 DAE models often contain many unknown parameters that must be determined from experimental data. If a model is not properly formulated, the parameter estimation problem may be illposed. Attempts to estimate the parameters in a poorly parametrized model will either fail or lead to nonunique parameter estimates. It is important, therefore, to check a priori whether parameters in a system are uniquely identifiable. A parameter set is identifiable if identical input-output behavior for a model implies identical values of the individual parameters for any two candidate parameter vectors on a dense and open subset of the set of possible parameter values. A model is identifiable if its parameters are identifiable. If the parameters are not identifiable, then they cannot be accurately estimated using input-output data, even using welldesigned experiments. This is the second of a pair of papers addressing the identifiability of linear time-invariant (LTI) DAE systems. Although LTI DAE systems arise less frequently in the modeling literature, their study has important applications to nonlinear DAEs. Specifically, it has been shown that a linearized LTI DAE system retains many of the qualitative aspects of the original nonlinear DAE system7,8 and that if the model parameters in the linearized DAE are locally identifiable, then the parameters in the corresponding nonlinear DAE model are also locally identifiable.9 A test for the identifiability of LTI DAE systems based on a generalization of the Markov parameter test for * To whom correspondence should be addressed. Tel.: (613)533-6343. Fax: (613)533-6637. E-mail: mclellnj@ chee.queensu.ca.
ordinary differential equations (ODEs) has recently been proposed.10 This test, which is based on partitioning of the DAE system into algebraic and ordinary differential equation parts, while always theoretically possible, is often computationally intensive. In addition, if a system is found to be not identifiable, the generalized Markov parameter (GMP) method does not provide a clear indication of how to reparametrize the system in terms of a reduced set of parameters that are identifiable. In the current paper, a differential algebra approach for assessing identifiability is developed that is computationally simpler than the GMP method. In addition, the proposed approach provides valuable information about how to reparametrize locally unidentifiable systems to make them locally identifiable. The proposed method produces an ODE realization of the original LTI DAE, which is of equal or lower dimension than the original system and has the same inputoutput behavior, using a minimal number of parameters and states. The class of models (or systems) under consideration is given by the following set of equations:
E(p) x′(t) + M(p) x(t) + B(p) u(t) ) 0 y(t) - C(p) x(t) ) 0 x(0) ) x0
(D)
where x(t) is the semistate vector, which lies in Rn. We will assume that system D has been expressed in deviation variables about some initial steady state, so that x0 ) 0. The output y(t) is an m-dimensional vector of outputs. The index of the DAE system D is ϑ. The input, u, is a Cϑ curve mapping a time value t in an open interval, u, containing the origin to a point u(t) in Rγ. The set of such inputs, which, in addition, allow the initial condition x0 ) 0, is denoted U. A further discussion of the properties of the inputs under consideration can be found in the Part 1 of the series.10 Any experi-
10.1021/ie030534j CCC: $27.50 © 2004 American Chemical Society Published on Web 01/31/2004
1252 Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004
ment done on system D will start at t ) 0. As a result, the output will be observed from t ) 0 to some strictly positive final time tf. Let u + be the set of all of the strictly positive elements of u. The parameter vector p is assumed to be in a subset P of Rs. We also assume that P is open and dense in a subset of Rs and that the system D is solvable for all p ∈ P.10,11 The quantities E(p), M(p), B(p), and C(p) are matrices of appropriate dimensions whose entries are continuously differentiable (C1) functions defined on the parameter space P. In addition, we assume that the n × n matrix E(p) is singular, the matrix C(p) has full row rank, and the matrix B(p) has full column rank. A brief review of key ideas from the study of DAEs is provided in section 2. In section 3, the system is expressed as a set of linear mappings in the variables x, y, and u. These mappings are then used to describe the input-output behavior of the system, which is used in the proposed test for identifiability presented in section 4. The problem of constructing a minimal (in the sense of the number of parameters and states used) ODE system based on the input-output properties of system D is treated in section 5. In section 6, the method is applied to a simple example involving an idealized gas-phase reactor. A comparison of the proposed method with the GMP method previously proposed,10 as well as a discussion of the application and limitations of this method to large-scale systems, is provided in section 7. 2. Identifiability and Input-Output Behavior Our proposed method is based on capturing the input-output behavior of a system using a set of equations generated by system D. The definitions in this section have previously been presented and discussed10 and are summarized below for completeness. The inputoutput behavior is defined as the relationship between a specified input curve and the value of the output at a final time, for a specified parameter vector and a specified final time. Definition 1 (Input-Output Behavior). The inputoutput behavior at x0 ) 0 of system D is the vectorvalued function ioD(p) mapping from the set U of possible input curves and the set u + of all possible final time values to the output space Rm. For a final time tf and an input curve u, y(tf) ) ioD(p)(u,tf). A specific instance of this mapping is referred to as an experiment on D. Two such mappings, io(r) and io(p) with r and p in P, are (pointwise) equal if they are equal for every input in U and final time in u +. That is, the two mappings are equal if they generate identical outputs for every possible experiment. The input-output behavior is used to define the notion of global identifiability. A system is identifiable if the relationship between the set of possible parameter values (points in P) and the set of possible input-output behaviors is one-to-one (injective) almost everywhere in P. The meaning of “almost” is made clear in the following definition. Definition 2 (Global Identifiability6). System D is globally identifiable if and only if, for almost any two candidate parameter vector values, r and p, both in P
ioD(p) ) ioD(r) S p ) r By almost any p and r in P we mean all p and r in an open and dense subset of P.
Locally, about some point p in P we have the following definition. Definition 3 (Local Identifiability6). For a prespecified input u, positive final time tf in u +, and parameter value p in P, system D is strongly locally identifiable around x0 if and only if
ioD(p) ) ioD(r) S p ) r for all parameter values r in a subset P h of P, where P h is an neighborhood of p given by the set of parameter vectors that are arbitrarily close to p and are in P, so that P h ) {r ∈ P| |r - p| < } for arbitrarily small positive . This definition is equivalent to algebraic identifiability as defined by Xia and Moog.5 All systems that are globally identifiable are also locally identifiable. However, local identifiability does not imply global identifiability.10 A comprehensive discussion of the various definitions associated with identifiability for ODE systems has been presented by Walter and Pronzato.6 3. Differential-Algebraic Approach In this work, an attractive alternative to the GMP approach10 for assessing the identifiability of LTI DAE systems is presented. This approach makes use of the techniques and concepts from differential algebra to convert the DAE system to an ODE system involving only the inputs and outputs of the system. The proposed approach is beneficial for several reasons. First, it is readily implemented using computer algebra software such as Maple. Second, the technique produces an equalor lower-dimensional ODE realization of the LTI DAE system that uses only the minimal number of parameters, and states, necessary to produce the observed input-output behavior. Finally, the new realization is an ODE system, which allows classical techniques for control and identification to be applied. This last fact is significant because it naturally leads to approaches for identification and experimental design for DAE systems. Differential-algebra based techniques have been used for control and identification of DAE as well as ODE systems. Fliess and co-workers12,13 pioneered much of the work, in the control literature, on the use of differential algebra (developed by Ritt14 and others). Important contributions to the use of differentialalgebraic techniques for testing the identifiability of mathematical models, made by Glad and co-workers15,16 and Diop,17 have been applied to models in the biosciences.18 Most recently, differential algebra has been used to treat the identifiability problem of ODE systems whose vector fields and output mappings are meromorphic (all of the singular points of the mappings are poles19).5 The method proposed in this work has several advantages over existing methods. It is easily implemented using computer algebra and is suitable for medium-scale systems (e.g., containing tens of states and parameters). In addition, because our method provides an ODE realization for the DAE system, it provides the user with intuition as to the underlying dynamic process. Furthermore, this realization is identifiable even if the original DAE system is not. In contrast with the work of Xia and Moog,5 we consider systems whose vector field and output equations are not necessarily meromorphic functions of the parameters. Finally, in contrast
Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004 1253
with the work done by Saccomani et al.,20 we do not require that the initial conditions be independent of the input. Our approach uses a combination of elementary row operations and differentiation of equations to eliminate the semistate variables from a subset of the model equations. This resulting set of equations is then expressed only in terms of the input and output variables (and their time derivatives). The task of testing a system for identifiability is simplified because the set of input-output equations generated by an LTI DAE system corresponds to an LTI ODE system. Identifiability tests for LTI ODE systems are well-known,3,6,21 and minimal realizations can be constructed.22 Some new notation needed for the discussion of differential algebra is presented here in the context of LTI DAE systems. Our goal is to use a repeated combination and differentiation of the equations defining system D to remove the dependence on the state variables. The state and output equations that make up system D are linear functions of the variables x′, x, y, and u. Consequently, the model can be considered to be a set of functions A, whose elements Ai(x′,x,y,u,p), i ) 1, 2, ..., n + m, map to zero along the trajectories of system D.
{
Ai(x′,x,y,u,p) ) n
γ
∑[E(p) x′ (t) + M(p) x (t)] + ∑B(p) u (t) ij
j
ij j
j)1
ij j
i ) 1, ..., n
j)1
n
yi-n(t) -
∑[C(p)
(i-n)jxj(t)]
i ) n + 1, ..., n + m
j)1
For a fixed set of parameter values, system D has now been defined as a set, A, of linear functions (mappings) in the variables x, u, and y and their time derivatives. As a result, each element of A can be viewed as the action of a linear function on x, y, and u. Let δ ) d/dt denote the differentiation operator. We may represent any element Ai of A as
Ai(x′,x,y,u,p) ) ωi1(p,δ) x1 + ωi2(p,δ) x2 + ... + ωi,n(p,δ) xn + ωi,n+1(p,δ) u1 + ... + ωi,n+γ(p,δ) uγ + ωi,n+γ+1(p,δ) y1 + ... + ωi,n+γ+m(p,δ) ym where the coefficients ωij, j ) 1, 2, ..., n + γ + m, are rational functions in δ, with coefficients that are C1 functions defined on P. Letting K be the field of all rational functions in δ, we can represent {ωi1, ..., ωi,n+γ+m} as a row vector Ωi ) [ωi1, ωi2, ..., ωi,n+γ+m], in the vector space Kn+γ+m so that
([ })
x Ai(x′,x,y,u,p) )Ωi(p,δ) y u
We can therefore associate with the model A a set of row vectors Ω ) {Ω1, Ω2, ..., Ωn+m}, which spans a vector subspace of Kn+γ+m. It is important to distinguish between the vectors Ωi and the functions Ai because each element Ai is identically zero along solutions of system D, while Ωi is simply a vector in Kn+γ+m. By focusing on Ωi, we have transformed a problem dealing with a set of linear functions in the variables x′, x, y, and u to a set of vectors in Kn+m+γ, which we can treat
using the linear algebraic notions of linear independence, span, and range space. Using this approach, the dynamics for system D are associated with a set of vectors in Kn+γ+m whose action on solution trajectories [x, y, u] are identically zero. Because Ω is a vector space over the field K of rational functions in δ that are also functions of the parameters, linear combinations of these vectors consist specifically of differentiation and addition using coefficients that may be functions of the parameters. The first n entries of the coordinate representation of the vectors correspond to the semistates. Consequently, vectors in Ω whose first n entries are zero are associated with functions Ae of the form
Ae ) ωe,n+1(p,δ) u1 + ... + ωe,n+γ(p,δ) uγ + ωe,n+γ+1(p,δ) y1 + ... + ωe,n+γ+m(p,δ) ym (3.1) As Ljung and Glad22 noted for the nonlinear case, there exists an infinite number of functions that can be generated by combination and differentiation of the elements of A. The challenge is to find the subset of these equations that can be associated with the inputoutput behavior. For the linear case, this can be done systematically by determining a basis for the subspace Ωio of Ω consisting of all vectors whose first n elements are zero. This subspace, in turn, has associated with it a set of functions Aio that have no dependence on the semistates. Elements of this set will consist of functions of the form given in eq 3.1. The set Aio represents the input-output behavior of system D. For the linear problem, we can also consider this set as a basis. By a basis for Aio we mean that the set of mappings Ae in the set has associated vectors Ωe that provide a basis for Ωio. Definition 4 (Differential Basis). A subset of Aio is a differential basis for Aio if its associated vector set is a basis for Ωio. This differential basis is a set of equations that completely represent the input-output behavior of system D, using a finite number of independent coefficients. Example 1. Consider the following system:
p1x′1(t) + x2(t) + p1p2u1(t) ) 0
(3.2a)
1 x′ (t) ) 0 p5 2
(3.2b)
p4x2(t) + ep1x3(t) + p6u1(t) ) 0
(3.2c)
p5x3(t) - y(t) ) 0
(3.2d)
x(0) ) 0
(3.2e)
p3x1(t) + p2x2(t) +
First consider the set P of allowable parameter values. It has been shown10 that system 3.2 is solvable if there exists a value of λ such that the function
1 (p λ2 + p1p2p5λ - p5p3)ep1 p5 1 is nonzero. In addition, P must be defined so that 1/p5 is defined for every vector of parameter values p ∈ P. Assuming that the parameters can physically take on any value, we choose P ) {p ∈ R5|p5 * 0, p1 * 0, p3 * 0}. Note that for solvability we only require that one of
1254 Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004
p1 or p3 be nonzero, but for simplicity we restrict P so that neither can be zero. The elements of set A are given by the left-hand sides of eqs 3.2a-3.2d. The set Ω is given by
[
{
]
δ Ω ) [p1δ 1 0 0 p1p2 ], p3 p2 + p 0 0 0 , 5
}
[0 p4 ep1 0 p6 ], [0 0 p5 -1 0 ]
The input-output behavior is found by eliminating the semistate variables from a subset of eqs 3.2a-3.2d. This is accomplished by finding a basis for the subspace Ωio consisting of vectors in Ω whose first three entries are zero. A canonical basis for Ω can be computed by finding the row-reduced echelon form of a matrix whose rows are elements of Ω. Recall that the nonzero row vectors of a row-reduced echelon matrix form a canonical basis for the row space of the matrix.24 This computation can easily be done using commonly available computer algebra software such as Maple. In this example, this basis for Ω is given by
{[
1 0 0 0
[
[
0 0 1 0
[
0 0 0 1
]
p1p2(δ + p2p5) , -p3p5 + p1p2p5δ + p1δ2 -p1p2p3p5 0 1 0 0 , -p3p5 + p1p2p5δ + p1δ2
]
p1p2p3p4p5 - p6p3p5 + p1p2p5p6δ + p1p62δ2 ep1(-p3p5 + p1p2p5δ + p1δ2)
]
,
-p5(p1p2p3p4p5 - p6p3p5 + p1p2p5p6δ + p1p62δ2) ep1(-p3p5 + p1p2p5δ + p1δ2)
]}
Because the first three coordinates of Ωio correspond to the semistates, the basis for Ωio is given by the last row vector
[
0 0 0 1
-p5(p1p2p3p4p5 - p6p3p5 + p1p2p5p6δ + p1p62δ2) ep1(-p3p5 + p1p2p5δ + p1δ2)
]
(3.3)
This vector corresponds to the linear mapping given by
p3p5 dy d2y (t) + p2p5 (t) + y(t) + 2 dt p1 dt p1p2p4 + p6 du p3 p52 u(t) - p2p52p6e-p1 (t) p1 dt p1e
(
)
p5p6e-p1
d2u (t) ) 0 (3.4) dt2
The left-hand side of eq 3.4 forms a differential basis for Aio. The elements of this basis will be used to determine if system 3.2 is identifiable. In this section, a computationally simple algorithm has been presented that produces a set of equations called a differential basis for Aio that completely describe the input-output behavior of system D. This linear differential basis is a special case of the more general notion of a characteristic set used to identify the input-output behavior of nonlinear DAE systems.15
The approach presented here takes advantage of the linear nature of system D to associate a vector space structure with the model equations. The input-output behavior of the system then corresponds to a subspace of the vector space associated with the system (Ωio ⊆ Ω). A basis for the input-output subspace Ωio is then used to generate a set of mappings called a differential basis for system D. This differential basis completely characterizes the input-output behavior of system D up to the value of the initial conditions for the outputs and inputs, which are assumed known. Furthermore, the elements of Aio are all linear with respect to the input variables, output variables, and their time derivatives. Indeed, the elements of Aio are nothing but ordinary linear differential equations whose identifiability properties can be tested most directly by the transfer function approach.6,21 This result is especially useful in light of the fact that every solvable (on P) DAE system of the form of system D, with a sufficiently differentiable input, can be transformed into an ODE system by the method outlined in this work. The condition of solvability may be tested by examining whether the matrix pencil of the system in question is regular (this also implies that the system in question has a representation where the matrices E(p) and M(p) are square). The fact that the input must, in general, be suitably differentiable comes from the fact that the transfer function corresponding to the transformed ODE system may not be strictly proper. For example, in eq 3.4 the highest order of the derivative is the same in the output and input. As a result, if the input is not at least once continuously differentiable, the output may not be well defined as a function of time. Other methods for computing ODE realizations for nonlinear DAE systems have been proposed by Pantelides,25 Kumar and Daoutidis,4 and Reich.7 However, these methods may be far too involved if the system under study is linear. In addition, Luenberger’s shuffle algorithm26 allows the construction of an ODE realization for LTI DAE systems. However, the proposed method is more computationally efficient. This is because, unlike the proposed method, the Luenberger algorithm does not explicitly identify the input-output behavior of the system. 4. Test for Identifiability In this section the proposed transfer-function-based identifiability test for LTI ODE systems6,21 will be discussed and then applied to transfer functions generated by LTI DAE systems. In general, LTI ODE systems can be in the inputoutput form
y - Q(δ,p) u ) 0 y(0) ) 0 where Q(δ,p) ∈ Km×γ is called the input-output transfer function matrix whose entries are rational functions of δ, with coefficients that are C1 functions defined on a dense subset of P. Let qij(δ,p) be the ijth entry of Q(δ,p), for i ) 1, 2, ..., m and j ) 1, 2, ..., γ. We can assume without any loss of generality that every such nonzero function qij(δ,p) has a denominator that is monic (has a leading coefficient equal to 1) and there are no pole-
Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004 1255
zero cancellations in qij. As a result, we may write each such function uniquely as
qij(δ,p) )
∆lijaij(p)
(4.1)
∆kijbij(p)
where ∆lij ) [1, δ, δ2, ..., δlij] and ∆kij ) [1, δ, δ2, ..., δkij] with lij and kij the highest powers of δ that appear (with nonzero coefficients) in the numerator and denominator, respectively, of qij(δ,p). Also, a(p) ) [a0,ij(p), a1,ij(p), ..., alij,ij(p)]T and b(p) ) [b0,ij(p), b1,ij(p), ..., bkij-1,ij(p), 1]T are vectors with entries that are coefficients of δ in qij(δ,p). Theorem 1. Consider the mapping
S:p f [a11, a12, ..., am,γ, b11, b12, ..., bm,γ]T System D is identifiable if and only if the mapping S, of the parameters into the coefficients in each transfer function, is injective on a dense and open subset of P. Proof. First, recall that, by definition, the inputoutput behavior of the original system is fully represented by the equations making up the basis for Aio. Also, because the elements of Q(δ,p) can be uniquely written as in eq 4.1, the mapping
S h :(a11, a12, ..., am,γ, b11, b12, ..., bm,γ) f ioD is injective. If the mapping S is itself injective, then the composition S h oS is also injective and the system is identifiable. Note that this composite mapping can always be defined because the range of S is the domain of S h . If the mapping S is not injective, then neither is the composition mapping because, by definition, the range of S is the domain of S h. 9 Locally, one can check whether the mapping S is injective by checking if the rank of the Jacobian matrix J ) ∂S/∂p is equal to s, the number of parameters in the system. In our example, we use eq 3.4 to obtain the mapping S given by
S(p) ) [a1,1, b1,1]T ) [a0(p)
a1(p)
a2(p)
b0(p)
b1(p)]T
where
a0(p) ) p3 p52
p1p2p4 + p6 , a1(p) ) -p2p52p6e-p1 p1 p1e
a2(p) ) -p5p6e-p1, b0(p) ) p2p5, b1 )
p3p5 p1
The value of S at any parameter value p depends on only five independent coefficients, which are expressed in terms of six parameters. As a result, the mapping S cannot be injective. Therefore, we conclude that system 3.2 is neither locally nor globally identifiable. 5. Reparametrization of Unidentifiable Systems In general, each entry qij(p,δ) of Q(p,δ) can be expressed as the product of a strictly proper rational function in δ and some nonnegative integer power of δ. Let lj be the lowest nonnegative power of δ such that if j ij(p,δ) δlj for i ) 1, 2, ..., m, then q j ij(p,δ) is qij(p,δ) ) q strictly proper as a function of δ. This means that if we
Table 1. Summary of the Proposed Procedure step
action
1 2 3 4 5 6
express system D as A, a set of linear mappings compute Ω, the set of associated vectors compute a canonical basis for span(Ω) identify elements on this basis that make up Ωio use Ωio to define a differential basis for Aio check the LTI ODE given by this differential basis for identifiability compute a minimal realization for the differential basis that is identifiable
7
h (p,δ) whose entries are q j ij(p,δ) let vj ) ujlj, the matrix Q ) qij(p,δ) δ-lj is a matrix of strictly proper functions. That is, the degree of δ in the denominator in each element exceeds the degree of δ in the numerator in that element. It has been shown22 that, given any such matrix Q h (p,δ), there exist constant matrices A(p), B(p), and C(p) such that C(p) [Iδ - A(p)]-1B(p) ) Q h (p,δ). Furthermore, there exists a minimal time-invariant realization for Q h (p,δ) with v(t) ) [v1(t), v2(t), ..., vγ(t)]T as the input. This realization provides a state-space representation for the input-output behavior of the LTI DAE system using a minimal number of states. By examination of the coefficients of the input-output representation, an identifiable and minimal realization can be produced. This is done by assigning the new parameters to be the coefficients appearing in the transfer functions. This procedure is illustrated using the previous example. It has been shown in section 4 that system 3.2 is not globally identifiable. As a result, we reparametrize system 3.2 in such a way that the reparametrized system is identifiable. We reparametrize system 3.2 by setting the new parameter vector r as follows: r1 ) a0(p), r2 ) a1(p), r3 ) a2(p), r4 ) b0(p), and r5 ) b1(p). Note that the elements of the vector 3.3 are not strictly rational and so lj ) 1 and v(t) ) [du/dt](t). Using eq 3.4, we may now write the reparametrized system as
d3y d2y dy d2v dv (t) + r (t) + r (t) - r1 (t) (t) r 4 3 2 3 2 2 dt dt dt dt dt r0v(t) ) 0 which admits the state-space representation
[
] [] []
0 1 0 0 1 x′(t) ) 0 0 x(t) + 0 v(t) 1 0 -r3 -r4 r0 y(t) ) r1 x(t) r2 This LTI ODE system is locally identifiable, controllable, and observable and, in addition, has the exact same input-output behavior as system 3.2. Our proposed method is summarized in Table 1. 6. Gas-Phase Reactor Example10 Consider the gas-phase reactions k1
A 98 D k2
D 98 L occurring in a continuously isothermal reactor. The equations describing these reactions are
1256 Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004
First, we must place system 6.1 in the form of system D. This means expressing system 6.1 in deviation variables. The steady-state equation for this system is given by
dCA(t) ) yA,in(t) F - k1CA(t) V - FyA(t) dt dCD(t) V ) k1CA(t) V - k2CD(t) V - FyD(t) dt dCL(t) V ) k2CD(t) V - FyL(t) dt dCN(t) V ) FyN,in(t) - FyN(t) dt yA(t) + yD(t) + yL(t) + yN(t) - 1 ) 0 V
M(p) xs + R(p) + B(p) us ) 0
Subtracting this equation from eq 6.1, we have
E(p) xj′(t) + M(p) xj(t) + R(p) + B(p) u j (t) ) 0 yj(t) ) C(p) xj(t) (6.2)
V )0 yD(t) - CD(t) ntotal
j ) u - us, and yj ) y - ys. where xj ) x - xs, u We wish to know if the parameter set p ) {k1, k2} in system 6.2 is identifiable. As before, we start with step 1 in Table 1 and we let Ai, i ) 1, 2, ..., 9, be the rows of system 6.2. We can therefore consider Ω to be the row span of the matrix
V )0 yL(t) - CL(t) ntotal
Eδ + M 0 B -C 1 0
V )0 yA(t) - CA(t) ntotal
[
yN,in(t) + yA,in - 1 ) 0
(I)
where Ci and yi, with i taking on values in the set {A, D, L, N}, are the volumetric concentrations and mole fractions, respectively, of the ith species. The volume is V, and the unknown parameters are k1 and k2. The input is the inlet mole fraction of A, yA,in(t). The inlet contains only inert gas, N, and feed, A. The inlet flow rate F is known and constant. Also, ntotal may be calculated according to the ideal gas law. We assume that the pressure and the temperature are known. In addition, we assume that ntotal is a positive integer and k1, k2, and V are all positive real numbers. Finally, we assume that only a measurement of CL(t) is available. Letting x(t) ) [CA(t), CD(t), CL(t), CN(t), yA(t), yD(t), yL(t), yN(t), yN,in(t)], y(t) ) CL(t), and u(t) ) yA,in(t), we can rewrite system I as
Ex′(t) + Mx(t) + R + Bu(t) ) 0 where
[
ys ) C(p) xs
-V 0 0 0 E) 0 0 0 0 0
0 -V 0 0 0 0 0 0 0
0 0 -V 0 0 0 0 0 0
0 0 0 -V 0 0 0 0 0
y(t) ) Cx(t)
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
]
0 0 0 0 0 , 0 0 0 0
(6.1)
Proceeding with steps 2-5 in Table 1, we find a canonical basis for Ω, given by Chart 1, where F ) F/ntotal, and use it to obtain a basis for Aio given by the linear function
d3yj d2yj dyj (t) + b (t) + b1 (t) + b0y(t) + a0u j (t) ) 0 2 dt dt3 dt2 (6.3) where
F a0 ) k1k2 V b0 ) F(Fk1 + F2 + Fk2)
[
b1 ) 2(Fk1 + F2 + Fk2) + k1k2 b2 ) 3F + k1k2
]
Note that the order of eq 6.3 is 3, which is far less than
-k1V k1V 0 0 0 V M) -n total
[ ] []
0 0 0 0 R ) -1 , B ) 0 0 0 -1
]
0 -k2V k2V 0 0
0 0 0 0 0
0 0 0 0 0
0
0
0
-
V ntotal 0
0
0
0
0
-F 0 0 0 1
0 -F 0 0 1
0 0 -F 0 1
0 0 0 -F 1
0 0 0 F 0
0 1
0
0
0
0 ,
0 0
1
0
0
0
0
1
0
0
0
0
0
1
V ntotal 0 0 0 0 0 -
F 0 0 0 0 , C ) [0, 0, 1, 0, 0, 0, 0, 0, 0] 0 0 0 1
Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004 1257 Chart 1
9, the number of semistates in system I, and even less than 4, the number of ODEs in the original model. We now proceed with step 6 in Table 1. Recalling that p ) [k1, k2], we may check locally if the system is identifiable by letting J ) ∂[a0,b0,b1,b2]T/∂p and checking the rank of J. In this case, rank(J) ) 2, and so the system is locally identifiable. However, the system is not globally identifiable because we can switch k1 and k2 and obtain the same transfer function.10 Note that the number of coefficients used in this approach is lower than the number of GMPs that must be considered to determine the identifiability in our prior approach.10 The fact that this system is locally identifiable but not globally identifiable suggests that we cannot reduce the number of parameters in the system. However, we may still construct a minimal ODE realization of the original DAE system. In this case, there are two parameters k1 and k2 and so only two of a0, b0, b1, and b2 are independent. Specifically, we choose r1 ) k1k2 and r2 ) k1 + k2. With this parametrization, we obtain the minimal realization
[
]
ξ˙ (t) ) 0 1 0 0 0 1 ξ(t) + -F2r2 - F3 -3F2 - 2Fr2 - r1 -3F - r2
[
yj(t) ) r1
F V
0
[]
0 0 u j (t) 1
]
0 ξ(t)
which is locally identifiable with the two parameters r1 and r2. Note that this parametrization does not change the number of parameters in the system. Instead, it changes the manner in which the parameters appear, making the system that was originally locally identifiable into a system that is globally identifiable. This is because r1 ) k1k2 and r2 ) k1 + k2 are the quantities that can be directly estimated from the input-output behavior of system I. This simple example illustrates the economy in states achieved by our method. In addition, as illustrated by this example, the proposed method provides insight as to the dynamical structure
of the model. The original model, although containing nine equations, has only three independent dynamical states. 7. Promise for Analyzing the Identifiability of Large-Scale DAE Systems Several techniques can be used to test the identifiability of large-scale DAE systems. Because of the computational complexity of nonlinear identifiability test methods5,15,17,20 and restrictions on the types of nonlinear systems that can be treated, it is advisable to first attempt a linearization-based approach. Local identifiability of the linearized system is a sufficient condition for the local identifiability of the corresponding nonlinear DAE system.9 We have proposed two methods for assessing the identifiability of LTI DAE systems. The first method10 involves partitioning the system into an ordinary part and an algebraic part, followed by the generation of GMPs, which are examined to determine the identifiability. The second method, proposed here, relies on differential algebra to generate a high-order ODE or transfer function model. The coefficients in this model are examined to determine the identifiability. For large-scale systems, the second approach is more appropriate because it is easier to implement using computer algebra software. In particular, the partitioning of the LTI DAE system, which may, in practice, involve complicated functions of the parameters, is avoided. The method described in this work relies on the computation of a canonical basis for Ω. The proposed procedure is attractive because it can be readily implemented using computer algebra software using built-in functions to determine the span and to compute bases for sets of vectors. Finally, it should be noted that the differential-algebraic approach generates a minimal ODE realization of the original DAE system. The approach outlined in this work can be computationally intensive for systems with more than 12 states because the computer resources required to deal with the many symbolic expressions can grow dramatically. However, there are several ways of mitigating the computational burden. First, when it is recognized that we are dealing with a linear dynamic system at this point, the identifiability test can be partitioned by examining the relationship between the parameters and subsets of the input-output transfer functions and their
1258
Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004
coefficients. The procedure involves determining which parameters are unidentifiable for each input-output pair. If a parameter or group of parameters is unidentifiable for every input-output pair, then that parameter or group of parameters is unidentifiable for the entire system. Such a procedure lends itself to the use of parallel computing, where the tasks of assessing the identifiability of parameters in each input-output pair can be delegated to parallel processors. The computational burden could also be alleviated by using numerical approximations for the Jacobian ∂S/∂p evaluated at some nominal parameter value, thereby avoiding the need to compute large-scale symbolic expressions with their associated storage requirements. An alternative to examining the transfer function coefficients would be to convert the ODE model generated using the differential-algebraic row-reduction approach to a linear first-order system of ODEs, so that the Markov parameter approach proposed by Grewal and Glover27 could be used to test for identifiability. This approach, however, is not recommended for large-scale systems because, like the GMP approach, it may involve very complicated functions of the parameters. Care must be taken when performing identifiability tests for DAE systems, especially when the number of parameters is large, because each element of the set Ω must be defined on P. This restriction affects the parameter set under which the identifiability test is performed. The set P is often hard to compute. If a nominal parameter value is used for a local test, it is necessary to ensure that each element of Ω is defined at p. 8. Conclusion An LTI DAE system can be considered a set of linear mappings in the states, inputs, and outputs. We may associate with this system a vector space, a subspace of which effectively involves only the input and output variables. By generating a basis for this subspace, we can assign a set of linear equations in the inputs, outputs, and their derivatives, which completely characterizes the input-output behavior of the original system. This set of equations is always linear and therefore has a state-space representation, which may be tested for identifiability using well-known approaches.6 In addition, this realization is minimal realization and has a solution that is easily computed. The proposed technique has been used to test both the local and global identifiability of a simple gas-phase reactor model. The method shows promise for testing the identifiability of larger-scale DAE models and is an attractive alternative to the GMP approach that we previously developed because it is easier to implement using computer algebra software and produces a minimal realization of the DAE model that can be used for simplification and reparametrization of DAE models that are not identifiable. Nomenclature Curves and Functions alij, bkij ) coefficients in the transfer function matrix Q; a set of functions defined on P, l ) 1, 2, ..., lij, k ) 1, 2, ..., kij, i ) 1, ..., m, j ) 1, ..., γ Aioi ) linear function in u, y, and their derivatives, i ∈ {1, 2}
Ae ) linear function in u, y, and their derivatives; an element of A Ai ) state or output equation in system D bi ) coefficients in the input-output equations of a system; a function defined on P, i ∈ Z+ Ci ) volumetric concentration of the ith species in the gasphase reactor example i ∈ {A, D, L, N} D ) linear differential-algebraic model ioD ) input-output map for system D u, x, y ) inputs, semistates, and outputs of system D u j , xj, yj ) inputs, semistates, and outputs of system 6.2 v ) vector of new inputs defined by vj ) ujδlj, j ) 1, 2, ..., γ; a continuous curve defined on u, whose range lies in Rγ yi ) mole fraction of the ith species in the gas-phase reactor example i ∈ {A, D, L, N, Ain} δ ) differentiation operator ζ ) states of the gas-phase reactor example in deviation form ν ) rational function defined on δ whose coefficients are C1 functions defined on P ξ ) set of states that are used in the minimal and locally identifiable realizations computed in this paper Vectors and Matrices aij, bij ) vectors whose entries are coefficients of powers of δ in qij, i ) 1, ..., m and j ) 1, ..., γ B, E, M ) coefficient matrices of the input, derivative of the semistate, and the semistate, respectively, in the state equation of system D; matrices whose entries are functions of p defined on P C ) coefficient matrix of the semistates x in the output equation of system D; an m × n matrix whose entries are functions of p defined on P J ) Jacobian of the coefficients in the input-output transfer function in example 3.2; functions of p defined on P J ) Jacobian of the coefficients in the input-output transfer function in the gas-phase reactor example; functions of p defined on P p ) vector of parameters that lies in P r ) vector of parameters that lies in P x0 ) initial value of x, assumed to be identically zero ∆k ) vector whose entries are increasing powers of δ, i.e., ∆kij ) [1, δ, δ2, ..., δkij] ωi ) coefficients of the state, input, and output variables in an element of A; also, by definition, these are the entries of an element in Ω; mathematically, these are a rational function in δ with coefficients that are C1 functions defined on P, i ) 1, 2, ..., n + γ + m Ωi ) associated vector for Ai, i ) 1, 2, ..., n + m Ωe ) associated vector for Ae Sets A ) set whose elements, A1 to An+m, are linear functions in the variables x, y, and u Aio ) set of equations generated by combination and differentiation of elements of A that involve only the inputs, outputs, and their time derivatives K ) field of rational functions in δ whose coefficients are C1 functions defined on P P ) open and dense subset of Rs P h ) open subset of P containing all of the points that are arbitrarily close to p Q ) matrix of transfer functions for which system D is a realization Q h ) matrix of strictly rational transfer functions relating the input v to the output y R ) set of all real numbers u ) open interval in R that contains the origin u+ ) set whose elements are the positive elements of u
Ind. Eng. Chem. Res., Vol. 43, No. 5, 2004 1259 U ) set of Cϑ inputs that allow x0 ) 0 Z+ ) set of all strictly positive integers Ω ) set of vectors associated with the elements of A; a set whose elements are ωi, i ) 1, 2, ..., n + m Ωio ) set of vectors associated with the elements of Aio Scalars ki ) kinetic parameters in the gas-phase reactor example, i ) 1, 2 lj ) degree to which the power of δ in the numerator exceeds the degree of δ in the denominator of an improper rational function in δ, j ) 1, 2, ..., γ F ) ratio of the feed rate to the total moles in the gasphase reactor example; F ) F/ntotal F ) feed rate in the gas-phase reactor example m ) number of outputs in system D n ) number of semistates in system D ntotal ) total number of moles in the gas-phase reactor example s ) number of parameters in system D t ) value of time in the set u tf ) element of u+; a time at which the system output may be observed V ) reactor volume in the gas-phase reactor example ) sufficiently small positive constant λ ) complex variable used to define the matrix pencil, λE(p) + M(p), of system D γ ) number of inputs in system D ϑ ) number of times the inputs to system D can be differentiated with respect to time
Literature Cited (1) Pantelides, C. C.; Gritsis, D.; Morison, K. R.; Sargent, R. W. H. The Mathematical Modelling of Transient Systems Using Differential-Algebraic Equations. Comput. Chem. Eng. 1988, 12, 449. (2) McLellan, P. J. A Differential-Algebraic Perspective on Nonlinear Controller Design Methodologies. Chem. Eng. Sci. 1994, 49, 1663. (3) Vajda, S. Structural Identifiability and the Realization Problem. Sci. Pap. Inst. Tech. Cybernetics Tech. Univ. Wroclaw 1985, 29, 228. (4) Kumar, A.; Daoutidis, P. Control of Nonlinear Differential Algebraic Equation Systems. In Research Notes in Mathematics; Brezis, H., Douglas, R. G., Jeffrey, A., Eds.; Chapman and Hall/ CRC: London, 1999. (5) Xia, X.; Moog, C. H. Identifiability of Nonlinear Systems with Application to HIV/AIDS Models. IEEE Trans. Autom. Control 2003, 2, 330. (6) Walter, E.; Pronzato, L. Identifiabilities and Nonlinearities in Nonlinear Systems. Model. Estim. 1995, 1, 111. (7) Reich, S. On the Local Qualitative Behavior of DiffernetialAlgebraic Equations. IEEE Trans. Autom. Control 1995, 14, 427. (8) Campbell, S. L.; Linearization of DAEs Along Trajectories. Z. Angew. Math. Phys. 1995, 46, 70. (9) Ben-Zvi, A. Identifiability of Differential Algebraic Systems. Ph.D. Dissertation, Queen’s University, Kingston, Ontario, Canada, in progress.
(10) Ben-Zvi, A.; McLellan, P. J.; McAuley, K. B. Identifiability of Linear Time-Invariant Differential-Algebraic Systems. 1. The Generalized Markov Parameter Approach. Ind. Eng. Chem. Res. 2003, 42, 6607. (11) Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solutions of Initial-Value Problems in Differential-Algebraic Equations; North-Holland: New York, 1989. (12) Fliess, M. Nonlinear Control Theory and Differential Algebra. In Lecture Notes in Control and Information Sciences; Descusse, J., Fliess, M., Isidori, A., Leborgne, D., Eds.; SpringerVerlag: Berlin, 1988. (13) Fliess, M.; Glad, T. An Algebraic Approach to Linear and Nonlinear Control. In Essays on Control: Perspectives in the Theory and its Applications; Trentelman, H. L., Willems, J. C., Eds.; Birkha¨user: Boston, 1993. (14) Ritt, J. F. Differential Algebra; American Mathematical Society: New York, 1950. (15) Glad, T. Nonlinear State Space and Input Output Descriptions Using Differential Polynomials. In Lecture Notes in Control and Information Sciences; Descusse, J., Fliess, M., Isidori, A., Leborgne, D., Eds.; Springer-Verlag: Berlin, 1988. (16) Glad, T.; Sokolov, A. Structural Identifiability: Tools and Applications. Proceedings of the 38th Conference on Decision & Control, Phoenix, AZ, Dec 1999; WePO6 16:40. (17) Diop, S. A State Elimination Procedure for Nonlinear Systems. In Lecture Notes in Control and Information Sciences; Descusse, J., Fliess, M., Isidori, A., Leborgne, D., Eds.; SpringerVerlag: Berlin, 1988. (18) Margaria, G.; Riccomogano, E.; Chappell, M. J.; Wynn, H. P. Differential Algebra methods for the study of the Structural Identifiability of Rational Function State-Space Models in the Biosciences. Math. Biosci. 2001, 174, 1. (19) Marsden, J. E. Basic Complex Analysis; W. H. Freeeman and Company: New York, 1973. (20) Saccomani, M. P.; Audoly, S.; D’Angio´, L. Parameter Identifiability of Nonlinear Systems: the Role of Initial Conditions. Automatica 2003, 39, 613. (21) Bellman, R.; Astrom, K. J. On Structural Identifiability. Math. Biosci. 1970, 7, 329. (22) Brockett, R. W. Finite Dimensional Linear Systems. In Series in Decision and Control; Howard, R. A., Ed.; John Wiley and Sons: New York, 1969. (23) Ljung, L.; Glad, T. On Global Identifiability for Arbitrary Model Parameterizations. Automatica 1994, 30, 265. (24) Hoffman, K.; Kunze, R. Linear Algebra, 2nd ed.; Prentice Hall: Englewood Cliffs, NJ, 1971. (25) Pantelides, C. C. The Consistent Initialization of Differential-Algebraic Systems. SIAM J. Sci. Stat. Comput. 1988, 9, 213. (26) Lewis, F. L. A Survey of Linear Singular Systems. Circuits Syst. Signal Process. 1986, 5, 3. (27) Grewal, M. S.; Glover, K. Identifiability of Linear and Nonlinear Dynamical Systems. IEEE Trans. Autom. Control 1976, 21, 833.
Received for review June 30, 2003 Revised manuscript received November 11, 2003 Accepted December 23, 2003 IE030534J
