Ind. Eng. Chem. Res. 2003, 42, 6607-6618
6607
Identifiability of Linear Time-Invariant Differential-Algebraic Systems. I. The Generalized Markov Parameter 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 types of experiments that are performed. A new method for testing the identifiability of linear time-invariant (LTI) differential-algebraic systems is presented. This method is an extension of the Markov parameter method, proposed by Grewal and Glover (IEEE Trans. Autom. Control 1976, 21, 833) and Vajda (Sci. Pap. Inst. Tech. Cybernetics Tech. Univ. Wroclaw 1985, 29, 228) for LTI ordinary differential equation systems. In the proposed method, the differential-algebraic system is partitioned into an ordinary differential equation subsystem and an algebraic subsystem. This allows for the computation of structural invariants called the generalized Markov parameters (GMPs). A system is identifiable if and only if the dependence of the GMPs on the parameters is one-to-one almost everywhere on the allowable parameter set. Tests for global and local identifiability are developed. The application of this method is demonstrated using two examples including an idealized gas-phase reactor. 1. Introduction Differential-algebraic equation systems (DAEs), also called implicit or semistate systems, are common in the study of mechanical and electrical systems, as well as chemical processes.2-6 The accurate estimation of model parameters is important for the design, scale-up, and optimization of these systems.2,7 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 might be ill-posed. 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 well-designed experiments. The class of models (or systems) under discussion in this work 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)
(D)
where x(t) is the semistate vector, which lies in a subset of Rn. The output y(t) is an m-dimensional vector of outputs. The input, u, is a Cϑ curve (a ϑ-times-differentiable curve) for some sufficiently large positive integer ϑ and maps a time value, t, in an open interval * To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: (613)533-6343. Fax: (613)5336637.
u containing the origin to a point u(t) in Rγ. We will assume that system D has been expressed in deviation variables about a steady state and that the system semistates are initially at this steady-state point, so that x(0) ) 0, u(0) ) 0, and y(0) ) 0. Let U denote the set of inputs that are Cϑ curves having u(0) ) 0. Note that not all sufficiently differentiable inputs allow the semistate to be at the origin at time zero. The input u must be ϑ times differentiable with respect to time so that the output, y, is always defined on u . This can be seen explicitly when the DAE system is partitioned into differential and algebraic subsystems, as shown section 3, and is in contrast to the case for ordinary differential equation (ODE) systems in which the input can be piecewise continuous. Any experiment done on system D will start at t ) 0. However, for the time derivatives of the input to exist at t ) 0, the input must be defined on an open interval containing the origin. The output will be observed from t ) 0 to some strictly positive final time tf lying in u +, defined as the set of all strictly positive elements of u . The parameter vector p is assumed to be in an open and dense subset P of Rs. To ensure that system D is solvable, parameter vectors must be excluded if they cause the matrix pencil (to be defined in section 2) for system D to be singular. P must be open in Rs so that, at any point p in P , the parameter vector can be perturbed by a small amount in any direction. As we shall see later, there can be other restrictions on the possible values of the parameter vector. Matrices E(p), M(p), B(p), and C(p) contain entries that are continuously differentiable (C1) functions defined on the parameter space P . These functions must be continuous and differentiable so that we can take the partial derivative of each entry with respect to each parameter. Also, we assume that the n × n matrix E(p) is singular; if E(p) were not singular, then system D could be written as an ODE system simply by leftmultiplying the state equation by E-1(p). Finally, we
10.1021/ie030317i CCC: $25.00 © 2003 American Chemical Society Published on Web 11/08/2003
6608 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003
assume that C(p) has full row rank and B(p) has full column rank, which places additional restrictions on P . Although DAEs that are linear in their input-output behavior are relatively rare in chemical engineering applications, such systems have important application in the study of nonlinear DAEs, which are much more common. It has been shown that the linearization of a nonlinear DAE about rest points8 and along trajectories9 retains much of the qualitative behavior of the original nonlinear DAE system. Although, in general, the linearization of DAE systems along trajectories is not timeinvariant, the linearization of a time-invariant nonlinear DAE evaluated at a rest point can be expressed as a linear time-invariant (LTI) DAE system. In addition, local identifiability of the linearized system about a rest point is a sufficient condition for local identifiability of the original nonlinear system.10 In section 2, we introduce some of the theory and terminology used in the study of DAEs. In section 3, system D is transformed into two subsystems, one of which is simply an ODE system, while the other is completely algebraic. An algorithm for doing this is presented in section 4. In section 5, the concept of Markov parameters that have been used to test the identifiability of LTI ODE systems1 is generalized to LTI DAEs. These generalized Markov parameters are used to determine necessary and sufficient conditions for the identifiability of system D. Sections 6 and 7 apply the generalized Markov parameter methods to a simple illustrative example and a more complicated example involving an idealized gas-phase reactor. Finally, the application of the generalized Markov parameter method to large-scale systems is discussed in section 8. 2. Identifiability and Input-Output Behavior Identifiability concerns the effect of parameters on the input-output behavior of a model, so first, we must define what we mean by the input-output behavior of system D. In this work, we are interested only in the forced response of the system. Some dynamical systems, such as batch reactors, are generally not started from a steady state and knowledge of the initial conditions can be critical for estimating parameters. Incorporating unknown initial semistate values into the identifiability analysis is the topic of ongoing research10 and is beyond the scope of this article. To ensure that system D is solvable, it must have a regular matrix pencil.11,12 The matrix pencil for system D is given by λE(p) + M(p), where λ is a complex variable. This matrix pencil is regular if det[λE(p) + M(p)] is not identically zero as a function of λ. As noted in section 1, we consider only parameter values that allow system D to be regular and insist that there exists at least one such parameter vector in the set P . The notion of the regularity of the matrix pencil is often used to mean solvability.11 This notion of solvability is weaker than the notion of regularity used by Kumar and Daoutidis,5 whose definition of regularity implies both that the system is solvable and that the underlying algebraic constraints are independent of the input. In our work, we allow the input to appear in the algebraic equations. The input-output behavior is defined as a parametrized relationship between a specified input curve and the value of the output at a final time. Our identifiability test compares the outputs at an arbitrary final time rather than over a time interval because it is simpler
to compare two points in Rm rather than two curves. Two output curves are equal if and only if they are equal at every point in time, so our approach is to consider the output at every final time in u +. Definition 1 (Input-Output Behavior). The inputoutput behavior at x(0) ) 0 of system D is the parametrized vector-valued 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). According to definition 1, the mapping io(p) involves selecting an input curve and a strictly positive final time at which to observe the output. 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. Definition 2 (Global Identifiability7). System D at x(0) 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 a dense and open subset of P . This assumption is required because, for some unusual parameter values, the behavior described by ioD evaluated at that parameter value might be trivial. For example, if the parameter corresponding to the output gain of a dynamic system is 0, there is only one possible input-output behavior (y ) 0), regardless of what the other parameter values are. This is an example of a degenerate case, where the values of some subset of the parameters are such that the input-output behavior of the system does not depend on the other parameters. In practical situations where we might have a good initial guess for the values of the parameters or an idea of their relative scale, we can require only a local version of the above identifiability definition. If we are are interested only in local identifiability, we could restrict our search for parameter vector pairs that produce the same outputs to an arbitrarily small neighborhood around the nominal parameter vector value p. Definition 3 (Local Identifiability7). System D at x(0) is locally identifiable about p if and only if, for a parameter value p in P ,
ioD(p) ) ioD(r) S p ) r for all parameter values r in a subset P h of P given by an neighborhood of p, i.e., the set of parameter vectors that are arbitrarily close to p and are in P
P h ) {r ∈ P | ||r - p|| < } for sufficiently small positive . All systems that are globally identifiable are also locally identifiable. However (as will be illustrated by the example in section 7), local identifiability does not imply global identifiability. A comprehensive discussion of the various definitions associated with identifiability for ODE systems has been presented by Walter and Pronzato.7 3. Transformation of the DAE System To investigate the input-output behavior of system D, it is useful to separate the intrinsic differential
Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6609
equation part of system D from the algebraic part using a coordinate transformation.11 In the transformed semistate vector, z ) [z1, z2]T, with z1 ∈Rn1 and z2 ∈Rn2, the elements of z1 are states in an ODE system, and the elements of z2 are algebraic. It has been shown13,14 that there exist matrices P(p) and Q(p) that are nonsingular on P such that
P(p) E(p) Q(p) )
[ [
I 0 0 J(p)
] ]
computing the matrices Q(p) and P(p), the following algorithm is proposed on the basis of an algorithm previously proposed in the literature.13,14 This algorithm achieves the goal of partitioning system D into a differential equation subsystem and an algebraic subsystem. First, collect terms in x and rewrite system D as
[E(p) dtd + M(p)]x(t) + B(p) u(t) ) 0 (4.1)
H(p) 0 P(p) M(p) Q(p) ) 0 I
y(t) ) C(p) x(t)
where J(p) is a matrix of nilpotency j > 0 [meaning that Jj(p) ) 0 and Jj-1(p) * 0]. The nilpotency of J is exactly equal to the index of the DAE.11 The relationship between the index and the properties of a DAE system have been discussed by Pantelides.3 Determining these matrices can be computationally intensive, especially for large-scale systems or systems that contain complicated functions of the parameters. This problem can be posed as a Jordan block decomposition problem.13,14 Using P(p) and Q(p), system D can be rewritten as11
z′1(t) + H(p) z1(t) + B1(p) u(t) ) 0
(3.1a)
J(p) z′2(t) + z2(t) + B2(p) u(t) ) 0
(3.1b)
y(t) ) C1(p) z1(t) + C2(p) z2(t)
z(0) ) 0
Next, transform the system so that the coefficient of x(t) becomes I, rather than M(p). If M(p) is invertible, this can be accomplished by premultiplying eq 4.1 by M-1(p), as illustrated in the example in section 6. However, as illustrated in the example in section 7, the matrix M(p) is often not invertible. In this case, we can use the fact that, because system D is solvable, the matrix E(p)c + M(p) is invertible for some number c. Equation 4.1 can be rewritten as
[E(p) dtd - E(p)c + E(p)c + M(p)]x(t) + B(p) u(t) ) 0 which becomes
{(dtd - c)E(p) + [E(p)c + M(p)]}x(t) + B(p) u(t) ) 0
(4.2)
(3.1c)
Now, premultiply eq 4.2 by [E(p)c + M(p)]-1 to obtain the desired result
where
z(t) ) Q(p) x(t)
{(dtd - c)[E(p)c + M(p)]
[C1(p), C2(p)] ) C(p)Q(p)
[ ]
B1(p) ) P(p) B(p) B2(p)
Equation 3.1a is a set of ODEs and has a solution for any initial condition and continuous input.11 Equation 3.1b depends only on the input and can be expressed as j-1
∑ k)0
(-1)kJk(p) B2(p)
dku (t) dtk
}
E(p) + I x(t) + -1
[E(p)c + M(p)] B(p) u(t) ) 0 (4.3)
and
z2(t) )
-1
(3.2)
Next, compute a matrix T(p), invertible on P , such that T(p) [(d/dt) - c][E(p)c + M(p)]-1E(p) T-1(p) is in a blockdiagonal form such as the Jordan normal form. (Note that computer algebra software can be used to do the computation.) In addition, because block-diagonal matrices are unique up to a permutation in the block order,15 we can require a permutation matrix, K, such that
(
)
d - c KT(p)[E(p)c + M(p)]-1E(p) T-1(p)K-1 ) dt H h (p,c) 0 d -c J h (p) 0 dt
The semistate z2(t) is therefore a function of only the input curve and its time derivatives. In addition, it can be seen from eq 3.2 that, in general, the input has to be j - 1 times differentiable for a solution z2 to exist. If u is not at least j - 1 times differentiable, z2 might not be defined at every value of t ∈ u , and hence, y might not be defined at every t in u . In addition, for the curve z′2 to be continuous, the input must be j times continuously differentiable. Indeed, we will assume, from this point on, that ϑ, the number of times that we must be able to differentiate u, is at least equal to j.
{(
4. Partitioning the System
([
We now provide an algorithm for placing system D in the form of system 3.1. However, instead of directly
(
)[
]
Substituting the transformation z(t) ) KT(p) x(t) into eq 4.3 and premultiplying by KT(p), we have
)[
] }
H h (p,c) 0 d -c + I KT(p) x(t) + J h (p,c) 0 dt KT(p) [E(p)c + M(p)]-1B(p) u(t) ) 0
where I denotes an appropriately sized identity matrix. After rewriting, we obtain
] [
])
H h (p,c) 0 I - cH h (p,c) 0 d + z(t) + J(p,c) dt I - cJ h (p,c) 0 0 KT(p)[E(p)c + M(p)]-1B(p) u(t) ) 0
6610 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003
step
action
1 2
Collect terms in x in system D (see eq 4.1) Rearrange system 4.1 so that the coefficient matrix of x(t) is the identity matrix Solve for T(p) and K (if necessary) to place system 4.2 into block-diagonal form Premultiply the state equation by matrix 4.4 to place system in the form of eq 3.1
3 4
Premultiplying by
[
H h -1(p,c) 0 0 [I - cJ h (p,c)]-1
]
] [
(4.4)
[ ]
])
B (p) I 0 H(p,c) 0 d + z(t) + 1 u(t) ) 0 B2(p) 0 J(p,c) dt 0 I
where
h (p,c)] H(c,p) ) H h -1(p,c)[I - cH J(c,p) ) [I - cJ h (p,c)]-1J h (p,c) and
[ ] [
(-1)j-1C2(p) Jj-1(p) B2(p)
(-1)j-2C2(p) Jj-2(p) B2(p) · · · -C2(p) J(p) B2(p)
)
C2(p) B2(p)
gives
([
[ ] [ ]
Theorem 1 (Test for Global Identifiability). System D is globally identifiable if and only if, given any parameter value p in an open and dense subset of P
Table 1. Summary of Transformation Procedure for System D
]
B1(p) H h -1(p,c) 0 ) KT(p)[E(p)c + B2(p) 0 [I - cJ h (p,c)]-1 M(p)]-1B(p)
Finally, the outputs are related to the new states by substituting the transformation z(t) ) KT(p) x(t) into the output equation, yielding
y(t) ) C(p)T(p)-1K-1z(t) ) [C1(p) C2(p)]z(t) System D has now been placed in the form of system 3.1. The original system has been divided into a differential equation and an algebraic subsystem. The partitioned form of 3.1 is used to generate the GMPs for system D. Table 1 summarizes the procedure that was described in this section. 5. Identifiability The GMPs of system D are defined from its partitioned form (eq 3.1) using the matrices J(p), H(p), B1(p), B2(p), C1(p), and C2(p) computed in section 4. Definition 4 (Generalized Markov Parameters). The generalized Markov parameters (GMPs) for system D are elements in the following set
{(-1)i1C2(p) Ji1(p) B2(p), C1(p) Hi2(p) B1(p)}
C1(p) B1(p)
C1(p) H(p) B1(p) · · · 2n1-1 C1(p) H (p) B1(p)
(-1)j-1C2(r) Jj-1(p) B2(r)
(-1)j-2C2(r) Jj-2(p) B2(r) · · · -C2(r) J(p) B2(r) C2(r) B2(r)
(5.1)
C1(r) B1(r)
C1(r) H(r) B1(r) · · · 2n1-1 C1(r) H (r) B1(r)
implies that p ) r, where r is a vector in P . A detailed proof of this theorem is available in Appendix A and proceeds in the following way: First, noting that the LTI DAE system 3.1 has the same input-output behavior as system D, an LTI ODE system whose input-output behavior is the same as system 3.1 can be constructed by augmenting the states with additional states representing derivatives of the inputs. The Markov parameters of this augmented ODE system are exactly the GMPs of system D. It has been previously shown1 that the Markov parameters completely characterize the input-output behavior of ODE systems. This, in turn, implies that the GMPs completely characterize the input-output behavior of system D. Finally, it is shown that only the first j + 2n1 GMPs need to be checked to assess identifiability. This is because any higher-order GMPs can be reconstructed from the first j + 2n1 GMPs. To check for identifiability, we need to solve eq 5.1 for the parameter vector r in terms of the parameter vector p and determine whether the only solution in some dense and open subset of P is r ) p. This condition can be expressed as
i1 ) 0, 1, 2, ..., j - 1; i2 ) 0, 1, 2, ...
S(p) ) S(r) w p ) r
We now show that these parameters play an analogous role, in the linear time-invariant DAE system, to the Markov parameters in the LTI ODE case. Indeed, the GMPs of system D completely characterize the input-output behavior of system D and can be used to test the identifiability of system D using the following theorem:
where S denotes the matrix of GMPs in eq 5.1. Alternatively, if we are interested only in local identifiability around some nominal value of the parameter vector p ) p0 in P , we can use the following rank test: Theorem 2 (Test for Local Identifiability). System D is locally identifiable about some nominal parameter value p0 in P if and only if
(5.2)
( [ ])
Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6611
Because system 6.1 is already in deviation variables, step 1 in Table 2 is complete. Step 2 is to to split system 6.1 into differential and algebraic subsystems using the method outlined in Table 1. First, collect terms in x in system 6.1
(-1)j-1C2(p) Jj-1(p) B2(p)
rank
∂ ∂p
(-1)j-2C2(p) Jj-2(p) B2(p) · · · -C2(p) J(p) B2(p) C2(p) B2(p)
)s
[E(p) dtd + M(p)]x(t) + B(p) u(t) ) 0
(5.3)
C1(p) B1(p)
y(t) ) C(p) x(t)
C1(p) H(p) B1(p) · · · 2n1-1 (p) B1(p) C1(p) H
In this case, M is nonsingular on a dense and open subset of P . Therefore, step 2 in Table 1 involves premultiplying the state equation by M-1(p) and rewriting the system as
p)p0
[M
-1
where s is the number of parameters. A proof of this theorem is presented in Appendix A. Condition 5.3 implies that any other value of the parameter vector sufficiently close to p0 will produce a different input-output behavior than that produced by p0. This type of analysis has an advantage over the global test for identifiability given in theorem 1, because it involves a matrix rank test and, hence, is easier to implement using symbolic computation software. The primary challenge is expansion creep, in which expressions become more complicated for higher-order GMPs. The method proposed in this work is summarized in Table 2. The most difficult steps are step 2, wherein the system is partitioned, and step 4b, where global identifiability of the system is checked. This procedure is illustrated using two examples.
(p) E(p)
action
1 2 3 4a 4b
Place system in deviation form Partition system into a dynamical and an algebraic part Compute the system GMPs Check local identifiability using the rank of the Jacobian ∂S/∂p Check global identifiability using condition 5.2
6. Illustrative Example Consider the single-input-single-output system
[ ] [ ] [ ]
0 1 0 p1p2 1 0 0 0 p5-1 0 x(t)′ + p3 0 p2 x(t) + 0 u(t) ) 0 p6 ep1 p4 0 0 0 0 (6.1) y(t) ) [0, 0, p5]x(t)
(6.2)
which is nonlinear in the parameters. First, note that
Note that this is equivalent to letting c in eq 4.2 be 0. Next, transform system 6.3 into the block-diagonal form used in eq 3.1. Using computer algebra software,16 we obtain the matrix
[
0 T(p) ) 1 1
p2p4p5e-p1 0 1
]
[ ]
0 0 1 K) 1 0 0 0 1 0
Details concerning the computation of the elements of T(p) and the remaining matrices for all examples are available in a Maple worksheet.16 Letting z(t) ) KT(p) x(t), system 6.1 can be rewritten as
[KT(p) M
-1
d + I z(t) + dt KT(p) M-1(p) B(p) u(t) ) 0
(p) E(p)T(p)-1K-1
]
y ) C(p)K-1T(p)-1z(t)
([
] )
or
-p1p4e-p1 0 0 d 0 0 1 dt + I z(t) + 0 0 0
det[λE(p) + M(p)] ) e p2 - (p1p2p4)λ
P ) {p ∈ R6| p1, p2, p4, p5 * 0}
p1-1 - p3p42p5e-2p1 p4e-p1 0
The matrix K used to reorder the Jordan block is given by
p1
which is nonzero on a dense and open subset of R6, and thus, this matrix pencil is regular. Also, note that system 6.1 has an index of 2 because the state x3 does not appear explicitly in the algebraic constraint. The restriction p5 * 0 is necessary for this system to be properly defined. As we shall see later, we also have to assume that p1, p2, and p4 are also nonzero. Assuming that the remaining parameters can physically take on any value, we have
]
y(t) ) C(p) x(t)
Table 2. Summary of Identifiability Test Procedure step
d + I x(t) + M-1(p) B(p) u(t) ) 0 dt (6.3)
[
[
]
(p6 - p1p2p4)e-p1 p2 - p3p4p5p6e-2p1 u(t) ) 0 p6e-p1 (6.4)
]
e2p1 - p1p3p42p5 ep1 p1p3p42p5 - e2p1 y) z(t) p2p4 p1p2p42 p1p2p42 with
z)
[
x1
(e2p1 - p1p3p42p5) x2 + p2p4p5e-p1x3 p4e-p1x2 2p1 p1e
]
T
6612 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003
We therefore have
J h)
[ ] 0 1 0 0
and H h (p) ) [-p1p4e-p1]. In addition, it can now be seen that, for y(t) to be defined for every value in P , the set P must be restricted to nonzero values of p1, p2, and p4. Step 4 in Table 1 is to premultiply eq 6.4 by matrix 4.4 to obtain
([
[ ]
])
] [
B (p) I 0 H(p) 0 d + z(t) + 1 u(t) ) 0 B2(p) 0 J(p) dt 0 I y ) [C1(p) C2(p)]z(t)
where
[ ]
H(p) ) -
ep1 p1p4
GMPs, the Jacobian ∂S/∂p can have at most a rank of 4, and this system is not locally identifiable. This means that it is impossible to estimate all six parameters simultaneously. Because local identifiability is a necessary condition for global identifiability, this system is not globally identifiable. In addition, the greatest number of independent parameters that this system can have is four, the rank of ∂S/∂p. By carefully examining the GMPs, we can reduce the number of parameters used in this model. Specifically, we can consider the ratio p6/p2 as a single parameter because the parameters p6 and p2 only appear in the GMPs as this ratio. Similarly, we can treat the product p3p5 as a single parameter. These two lumped quantities can replace the four parameters p2, p3, p5, and p6 and thus provide a system with only four parameters that has the same input-output behavior as system 6.1 and is identifiable. 7. Gas-Phase Reactor Example Consider the following gas-phase reaction
B1(p) ) [p1p4e-2p1(p1p2p4 - p6)] C1(p) )
[
[
]
e-2p1 - p1p3p42p5 p1p2p42
J(p) )
B2(p) )
k1
A 98 D
[
[ ]
k2
D 98 L occurring in a continuous isothermal reactor. The equations describing this reaction are
0 1 0 0
p2 - p3p4p5p6e-2p1 p6e-p1
]
2p1 2 ep1 p1p3p4 p5 - e C2(p) ) p2p4 ppp2 1 2 4
]
Note that J(p) is a 2-by-2 superdiagonal matrix, and therefore, it has nullity j ) 2. System 6.1 is therefore an index-2 problem, with n1 ) 1 and n2 ) 2, implying that there are four GMPs for this system. This is less than the six Markov parameters for an ODE of the same dimensionality and much less than the 10 Markov parameters for the augmented system (using the input as a state). Step 3 in Table 2 is to compute the system GMPs. The GMPs are given by
[
]
C2(p) B2(p) -C2(p) J(p) B2(p) S(p) ) ) C1(p) B1(p) C1(p) H(p) B1(p)
[
(p1p2p4 - p6)e
p1
p1p2p42
p6 p2p4 p1p3p4p5p6 -2p1 p6 -p12p3p42p5e-2p1 + e + p1 p2 p2p4 p13p3p43p5e-3p1 -
p12p3p42p5p6 -3p1 p1p6 -p1 e - p12p4e-p1 + e p2 p2
]
The final step in Table 2 is to determine whether the mapping S, whose entries are the system GMPs, is injective using condition 5.2. Because there are four
V
dCA(t) ) yA,in(t)F - k1CA(t)V - FyA(t) dt
V
dCD(t) ) k1CA(t)V - k2CD(t)V - FyD(t) dt V
dCL(t) ) k2CD(t)V - FyL(t) dt
V
dCN(t) ) FyN,in(t) - FyN(t) dt
yA(t) + yD(t) + yL(t) + yN(t) - 1 ) 0
(7.1)
V )0 yA(t) - CA(t) ntotal V )0 yD(t) - CD(t) ntotal V )0 yL(t) - CL(t) ntotal yN,in(t) + yA,in - 1 ) 0 where Ci and yi are the volumetric concentrations and mole fractions, respectively, of species i for i ) A, D, L, and N. The volume is V, and the unknown parameters are the kinetic rate constants 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 can be calculated according to the ideal gas law. The reactor pressure and temperature are assumed to be known. For the model to make physical sense, ntotal, k1, k2, and V must all be
Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6613
positive real numbers. Finally, we assume that only a measurement of the product concentration, 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) ) yL(t), and u(t) ) yA,in(t), we can rewrite system 7.1 as
Ex′(t) + Mx(t) + R + Bu(t) ) 0
where
[
M)
[
y(t) ) Cx(t)
-V 0 0 0 E) 0 0 0 0 0
-k1V k1V 0 0 0 V ntotal 0
0 -V 0 0 0 0 0 0 0
0 0 0 -V 0 0 0 0 0
0 -k2V k2V 0 0
0 0 0 0 0
0 0 0 0 0
0
0
-
V ntotal 0
0
0
0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
-F 0 0 0 1
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
]
det(λE + M) ) -V4 λ(ntotal λ + F)(ntotal λ + ntotal k2 + F)(F + ntotal λ + ntotal k1) ntotal3
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
-
F 0 0 0 B) 0 0 0 0 1
0 0 0 0 0 0 0 0 0
0 -F 0 0 1
V ntotal 0 0 0 0 0
[ ] [] 0
0 0 0 0 R ) -1 0 0 0 -1
0 0 -V 0 0 0 0 0 0
(7.2)
We wish to know whether system 7 with the parameter set p ) {k1, k2} is identifiable. Because both k1 and k2 are positive constants, we let P ) {(k1, k2) ∈ R2| k1 > 0, k2 > 0}. We now check the identifiability of system 7 using the generalized Markov parameter approach. Step 2 in Table 2 is to partition this system into intrinsic differential and algebraic parts as shown in eq 3.1. First, we check that the system is solvable. We have
]
which is not identically zero as a function of λ, given the physical constraints on the parameters and other system variables. Step 1 in Table 1 is to put system 7 in the form of system 4.1
[E(p) dtd + M(p)]xj(t) + B(p) uj (t) ) 0
(7.4)
yj(t) ) C(p) xj(t)
(7.5)
Next, transform eq 7.4 to replace the matrix M(p) with the identity matrix. Note that, unlike in system 6.1, the matrix M(p) in this example is not invertible. Instead, we use some c * 0 such that E(p)c + M(p) is nonsingular. We express the matrix [E(p)(d/dt) + M(p)] as (d/dt - c)E(p) + [E(p)c + M(p)]. We can now left-multiply system 7.1 by [E(p)c + M(p)]-1 to obtain the following system
{(dtd - c)[E(p)c + M(p)]
-1
}
E(p) + I xj(t) + [E(p)c + j (t) ) 0 M(p)]-1B(p) u
yj(t) ) C(p) xj(t) C ) [0, 0, 1, 0, 0, 0, 0, 0, 0]
Step 1 in Table 2 is to express system 7.2 in term of deviation variables. The steady-state semistate xs, output ys, and input us are related by
Mxs + R + Bus ) 0
(7.3)
ys ) Cxs Subtracting eq 7.3 from eq 7.2, we have
E xj′(t) + Mxj(t) + Bu j (t) ) 0 yj(t) ) Cxj(t) where xj(t) ) x(t) - xs, u j (t) ) u(t) - us, and yj(t) ) y(t) ys.
Step 3 in Table 1 is to compute the matrix T(p) such that T(p)[E(p)c + M(p)]-1E(p) T-1(p) is a matrix in Jordan form. This key step can be accomplished with the aid of computer algebra software such as Maple. However, the matrix produced by Maple is unique only up to a permutation of the Jordan blocks. By premultiplying by a permutation matrix K and postmultiplying by K-1, we can place the Jordan blocks in a unique manner. Using the matrices
[ ]
0 0 0 0 K) 1 0 0 0 0
1 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 1
[
6614 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 1
0
0 1
ntotal 0 V 0 0 0 0 0 0
0
0 0
0
0
0
0
1
1 0 n 0 total V
0
0
0
0
1
0
0
0
-1 1
0
0
0
-1 1
0
0
0
0 -
0
1 1
0 1 k1 - k2 1 k1 1 1 T(p) ) V -1 - n total V 0 n total V 0 n total V 0 -n total
0
V ntotal 0 0 V 2 ntotal 0 1 V ntotal 0 -1 -
0
0
0
0 0
0 0
0 0
-1 1
]
[]
0 0 B2(p) ) 0 0 1
C2 ) [0 0 0 0 0 ]
and following the steps described in section 4, we obtain the matrices
[
H h (p) ) ntotal 0 cntotal + ntotalk1 + F
0
0
0
ntotal 0 cntotal + F
0
0
0
0
0
ntotal 0 cntotal + ntotalk2 + F 1 0 c
[ ]
0 0 J h (p) ) 0 0 0
[
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
For this example, with c ) 1, we have
F ntotal 0
k1 +
H h (p) )
0
0
0
F ntotal 0
0
0
0
k2 +
0
0
0
[]
F ntotal 0 0
F V F B1(p) ) V F V 0
C1(p) )
[
k2 k1 1 0 k2 - k 1 k1 - k2
[ ]
0 0 J(p) ) 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
]
]
]
Remark 1. Note that, J(p) ) (0); indeed, J is always equal to (0) if the system is of index 1 because the nilpotency of J is exactly equal to the system index.11 This is a useful way to check that the matrix J(p) obtained is the correct one. Step 3 in Table 2 is to compute the GMPs of the system. The number of states in the dynamical subsystem is n1 ) 4, and the number of algebraic states is n2 ) 5. Because the index is j ) 1, we need to consider only the first j + 2n1 ) 9 GMPs of system 7.1. The first nine GMPs for this system are {C2B2, C1B1, C1HB1, C1H2B1, ..., C1H7B1}. Letting p ) [k1, k2], we can consider the Jacobian matrix
[ ]
C2B2(p) C1B1(p) C 1HB1(p) ∂ J ) 2 ∂p C1H B1(p) l 7 C1H B1(p)
Finally, step 4 in Table 2 is to check whether the mapping S is injective. Locally, this is done by checking the rank of the Jacobian matrix ∂S/∂p. Because C1H2B1(p) ) -(k2k1F/V) and C1H3B1(p) ) -(ntotalk1 + 3F + ntotalk2)k2k1F/ntotalV, the rank of J is 2, which is the number of parameters. The parameters k1 and k2 are therefore locally identifiable on a dense and open subset of P . Physically, this means that there is at most one vector of parameter values in a region around the initial parameter guess for any given input-output behavior. Globally, we can check to see whether S is injective using condition 5.2. Note that, in GMPs C1H2B1(p) and C1H3B1(p), exchanging the values of k1 and k2 does not affect the value of the GMP. As can be seen by examining the rest of the GMPs,16 k1 and k2 can be interchanged without altering the GMPs. The parameters are, therefore, not globally identifiable. 8. Conclusion The method proposed in this work allows for the testing of LTI DAE models for identifiability. Such testing allows modelers to determine a priori whether it will be possible to obtain unique estimates of the model parameters. The main effort in using the GMP approach to identifiability is the computation of the matrices H(p) and J(p) that are used to partition the system into a differential and an algebraic part. This might require considerable computational effort for large-scale systems because, in general, it constitutes a problem of finding the Jordan normal form of the matrix pencil λE(p) + M(p)13,14 in terms of the entries of E(p) and M(p), which might be complicated functions of the parameters. In addition, when dealing with large-scale systems, it might be difficult to check whether the mapping S, containing the first j + 2n1 GMPs, is injective. However, locally, S is injective if it has a Jacobian that is full rank
Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6615
at some nominal value of the parameter vector. Local injectivity is therefore relatively easy to check, even for large systems, except for the problem of expansion creep as higher-order derivatives are computed. Finally, it is important that the matrix T(p) be nonsingular on P . Hence, when using a local test for identifiability, one must always check that the value of the nominal parameter vector is in P . Acknowledgment The authors gratefully acknowledge the financial support of the Province of Ontario (OGS and OGSST scholarships for Amos Ben-Zvi), Queen’s University (Queen’s Graduate Awards for A.B.-Z.), and the Natural Sciences and Engineering Research Council of Canada (Discovery Grants). Nomenclature Curves and Functions 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 ) input, state, and output curves, respectively, in system D us, xs, ys ) steady-state input, semistate, and output, respectively, of the gas-phase reactor system u j ) input, semistate, and output, respectively, of the gasphase reactor system in deviation variables υ ) input signal in system A.2; 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} z ) image of the curve x under the transformation z(t) ) Q-1(p) x(t) z1, z2 ) dynamic and algebraic, respectively, semistates in system 3.1 β ) mapping expressing the system parameters, or a subset thereof, as a function of the GMPs or of other system parameters ξi ) ith time-derivative of the input for i ) 1, 2, ..., j - 1
K ) matrix that permutes the Jordan blocks of T(p)[E(p) + cM(p)]-1M(p)T-1; nonsingular n × n matrix whose entries are either 0 or 1 Q, P ) invertible n × n matrices whose entries are functions of p defined on P p ) parameters in system D; vector of parameters that lies in P p j , p˜ ) partitioning of the parameter set into two disjoint sets p j and p˜ R ) vector in the idealized gas-phase reactor example used to express the system in deviation variables; vector on R9 whose entries are functions of the system parameters r ) alternative parameter set for system D; vector of parameters that lies in P S ) matrix whose entries are the generalized Markov parameters of system D S h ) matrix whose entries are the Markov parameters of system A.2 T ) unitary matrix used to express the matrix [E(p) + cM(p)]-1M(p) in a Jordan form; n × n unitary nonsingular matrix whose entries are functions of the parameters defined on P V1 ) n1 × jγ matrix whose entries are functions of p defined on P V2 ) n2 × jγ matrix whose entries are functions of p defined on P V3 ) jγ × jγ matrix whose entries are either 0 or 1 V4 ) jγ × γ matrix whose entries are either 0 or 1 Sets C i ) set of all i times differentiable functions, i ∈ Z+ D ) LTI DAE system that is tested for identifiability in this work; set of equations involving the semistate x, the input u, and the output y P ) set of all possible parameter values, assumed to be open and dense in some subset of Rs P h ) open subset of P containing all points that are arbitrarily close to p R ) set of all real numbers u ) set of all possible values of time; open interval in R that contains the origin u + ) set whose elements are the positive elements of u U ) set of Cϑ inputs that expressed in deviation variables Z+ ) set of all integers
Vectors and Matrices
Scalars
A h, B h, C h ) matrices making up system A.2 A ˜, B ˜, C ˜ ) matrices making up an LTI ODE system B ) coefficient of the input u in system D; n × γ matrix whose entries are functions of p defined on P B1, B2 ) coefficient of the input u in eqs 3.1a and 3.1b, respectively C, E, M ) matrices making up system D, whose entries are functions defined on P C1, C2 ) coefficients of the dynamic state z1 and the algebraic state z2, respectively, in the output equation of system 3.1 H ) coefficient of the dynamic state, z1, in eq 3.1a; n1 × n1 matrix whose entries are functions of p defined on P H h ) first (top-left) matrix in the quasidiagonal form of [E(p) + cM(p)]-1M(p); n1 × n1 matrix whose entries are functions of p defined on P Iγ ) identity on Rγ J ) coefficient of the time derivative of the algebraic state z2 in eq 3.1b; n2 × n2 matrix whose entries are functions of p defined on P J h ) second (bottom-right) matrix in the quasidiagonal form of [E(p) + cM(p)]-1M(p); n2 × n2 matrix whose entries are functions of p defined on P J ) Jacobian of the GMPs in the gas-phase reactor example
c ) complex number such that the matrix E + cM is nonsingular i1, i2 ) indices used in the context of Markov parameters and GMPs; positive integers j ) degree of nilpotency of the matrix J ki ) kinetic parameters in the gas-phase reactor example, i ) 1, 2 F ) feed rate in the gas-phase reactor example m ) number of outputs in system D n ) number of semistates in system D n1 ) number of dynamic states in system 3.1 n2 ) number of algebraic states in system 3.1 ntotal ) total number of moles in the gas-phase reactor example n j ) dimension of the state space of system A.2; integer equal to n + jγ s ) number of parameters in system D sj ) rank of Sp(p)p0; nonzero and full rank t ) point in time; value of time in the set u tf ) time at which the system output can be observed; element of u + V ) Reactor volume in the gas-phase reactor example w ) fixed value for the GMPs of system D ) sufficiently small positive constant γ ) number of inputs in system D
6616 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 ϑ ) number of times the inputs to system D can be differentiated with respect to time λ ) complex number benditemize
V1(p) ) [B1(p) 0 · · · 0]
[ ]
V2(p) ) [0 (-1)0J 0(p) B2(p) · · · (-1)j-2Jj-2(p) B2(p)]
Appendix A. Proofs of Theorems
0 Iγ 0 0
A.1. Proof of Theorem 1. The proof proceeds as follows: Lemma 1 shows that there exists an ODE realization for system D and that this realization has the same input-output behavior as system D. After reviewing the definition of Markov parameters for an LTI ODE system, it is shown in theorem 3 that the GMPs of system D are the Markov parameters of the ODE realization of system D. In lemma 2, we show that, in the case of DAE systems, the input-output behavior of the system is completely characterized by the first j + 2n1 GMPs. It is then shown in lemma 3 that, to test for identifiability, we need to consider only the effect of the parameters on these first j + 2n1 GMPs defining the mapping S from the parameter space to this subset of GMPs. It is then shown that system D is strongly globally identifiable if and only if the mapping S is oneto-one on a dense and open subset of P . Lemma1 (ODE Realization of System D). System D has an augmented LTI ODE realization that has the same input-output behavior as system D. Proof. By differentiating eq 3.2, we can write system 3.1 as an augmented ODE model
0 0 Iγ 0 · · · · ·· V3 ) · 0 0 0 · · 0 0 0 · · · 0 0 0 0 · · · 0
(-1)kJk(p) B2(p) ξk+1(t) +
with Iγ denoting the identity on Rγ. The fact that systems D and A.2 have the same input-output behavior is a consequence of the fact that the transformation Q(p) is invertible on P and the definition of vk(t) in terms of uk(j)(t). 9 Theorem 3. The Markov parameters of system A.2 are given by
{(-1)i1C2(p) Ji1(p) B2(p), C1(p) Hi2(p) B1(p)} i1 ) 0, 1, 2, ..., j - 1;
(-1)j-1Jj-1(p) B2(p) v(t) i ) 0, 1, ..., j - 2
(A.1)
ξ′j-1(t) ) v(t) y(t) ) C1(p)z1(t) + C2(p)z2(t) z(0) ) 0 where ξi(t) ) u(i)(t) ) [u1i, u2i, ..., uγi]T, i ) 0, 1, 2, ..., j 1, and v ) u(j). We can rewrite system A.1 in matrix form as
[ ] [ ] [ ]
V3i ) 0
2.
V2(p)V3j-1V4 ) 0
3.
V1(p)V3iV4 ) 0
4.
V1(p)V3j-1V4 ) B1(p)
[
[
where j is the nilpotency of the matrix J in eq 3.1b. Now, consider the Markov parameters of the ODE system A.1
C h (p) B h (p) ) (-1)j-1C2(p) Jj-1(p) B2(p) h (p) ) C h (p) A h i(p) B
C1(p)Hi-l-1(p)V1V3lV4 + ∑ l)0 C2(p) V2(p)V3i-1V4 i ) 0, 1, 2, ...
z1(t) y(t) ) C h (p) z2(t) ξ(t)
H(p) 0 V1(p) A h (p) ) 0 0 V2(p) 0 0 V3
∀i ∈ Z+, i * j - 1
i-1
z1(t) z1(t) d z2(t) ) A h (p) z2(t) + B h (p) v(t) dt ξ(t) ξ(t)
where
i2 ) 0, 1, 2, ...
∀i ∈ Z+, i g j
1.
k)0
ξ′i(t) ) ξi+1(t)
0
Proof. Note the following facts
j-2
∑
0 · · · Iγ
V4 ) [0 0 · · · 0 Iγ ]T
z′1(t) ) - H(p) z1(t) + B1(p) ξ0(t) z′2(t) )
· · · 0
(A.2)
]
0 B h (p) ) (-1)j-1Jj-1(p) B2(p) V4 C h (p) ) [C1(p) C2(p) 0]
For i < j, we have
h (p) ) C2(p) V2(p)V3i-1V4 ) C h (p) A h i(p) B (-1)j-1-iC2(p) Jj-1-i(p) B2(p) i ) 0, 1, 2, ..., j - 1
]
For i ) j, we have j-1
C h (p) A h i(p) B h (p) )
C1(p) Hi-l-1(p)V1V3lV4 ∑ l)0
) C1(p) H 0(p)V1V3j-1V4 ) C1(p) B1(p)
Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003 6617
For i > j, we have j-1
C h (p) A h i(p) B h (p) )
C1(p) Hi-l-1(p)V1V3lV4 ∑ l)0
) C1(p) Hi-j(p)V1V3j-1V4 ) C1(p) Hi-j(p) B1(p) i ) j + 1, j + 2, ...
(A.3)
Thus, the Markov parameters for system A.2 are given by the Markov parameters of system 3.1a and the coefficients of u(i)(t), i ) 0, 1, ..., j - 1, in eq 3.2, which, together, are the GMPs for system D. 9 Recall that the Markov parameters for an LTI ODE system are given by
x′(t) ) A ˜ (p) x(t) + B ˜ (p) u(t) y(t) ) C ˜ (p) x(t) x(0) ) 0 where x(t) is the n-dimensional state vector. p, u(t), and y(t) are defined as in section 1, and A ˜ (p), B ˜ (p) and C ˜ (p) are appropriately sized matrices whose entries are C1 functions defined on P . The Markov parameters of this LTI ODE system are given by the set
˜ (p)} {C ˜ (p)A ˜ i(p)B i ) 0, 1, 2, ... These parameters correspond to the time derivatives of y as t approaches 0+ for an impulse input and zero initial condition.17 These parameters have previously been used to test the identifiability of LTI ODE systems.1,2 The Markov parameters of an LTI ODE system are the entries in the Hankel matrix of the system and are invariant with respect to choice of realization.18 As a result, the Markov parameters are often referred to as invariants for the system. For the LTE ODE realization of system D in A.1, we have the Markov parameters C h (p) A h i(p) B h (p), i ) 0, 1,..., 2(n) - 1, where (n) is the dimension of the LTI ODE realization. Now consider the mapping S h : P f Rm(2nj )×γ given by
S h : p |f
[
C h (p) B h (p) C h (p) A h (p) B h (p) · · · 2n j -1 C h (p) A h (p) B h (p)
for i ) 0, 1, 2, ..., 2n1 - 1, with known coefficients. The first j terms arise as a consequence of the algebraic equations in the DAE model. Therefore, using the first j + 2n1 Markov parameters, we can reconstruct any higher-order Markov parameter. 9 As a result, we need to consider only a finite subset of the GMPs to test system D for identifiability. Consider the mapping S: P f Rm(j+2n1)×γ given by
]
which takes a parameter vector and produces a set of Markov parameters for the ODE realization of system D. Grewal and Glover1 and Vajda2 have previously shown that an ODE system is identifiable if and only if the mapping S h is one-to-one (injective) on a dense and open subset of P . Lemma 2. Any Markov parameter of system A.2 can be uniquely computed with knowledge of only the first j + 2n1 parameters. Proof. By eq A.3 and the Cayley-Hamilton theorem, any Markov parameter for the differential subsystem of the form C1(p) Hk(p) B1(p) with k g 2n1 can be written as a linear combination of the elements C1(p) Hi(p) B1(p)
[ ] (-1)j-1C2(p) Jj-1(p) B2(p)
(-1)j-2C2(p) Jj-2(p) B2(p) · · · 1 (-1) C2(p) J1(p) B2(p)
S: p |f
C2(p) B2(p)
C1(p) B1(p)
C1(p) H1(p) B2(p) · · · 2n1-1 C1(p) H (p) B2(p)
Lemma 3. The mapping S is injective if and only if the mapping S h is injective. Proof. By lemma 2, the elements of S constitute the independent set of Markov parameters, and consequently, rank(S) ) rank(S h ). 9 Proof of Theorem 1. If each parameter vector produces a unique set of Markov parameters, then the mapping S is injective, so that by lemma 3, S h must also be injective, and the system is identifiable. However, if condition 5.1 does not hold, then two different parameter values produce identical Markov parameters, and the mapping S is not injective, so that, by lemma 3, neither is S h , and the system is not identifiable. 9 A.2 Proof of Theorem 2. This proof uses the implicit function theorem and follows an approach similar to the proof presented in ref 2. Let w ∈ Rm(2n1+j)×γ be any fixed set of GMPs for system D and consider the continuously differentiable (with respect to p) mapping S(p) - w ) 0. By assumption, the Jacobian of this mapping given by ∂S(p)/∂p|p)p0 is full rank. It follows from the implicit function theorem that, for some neighborhood P h of p0, p0 ) β(w) holds for some function β. Therefore, the parameters can be obtained by observing the GMPs of system D. Because these GMPs uniquely characterize the input-output behavior of system D, the parameters themselves can be obtained from the input-output behavior of system D on P h . Consequently, the parameters are locally identifiable. Now consider the case where the mapping evaluated at p0, ∂S(p)/∂pp)p0, is not full rank. If the rank is 0, then none of the parameters can be estimated from inputoutput behavior, because no parameters affect the input-output map. If the rank of ∂S(p)/∂pp)p0 is some constant sj < s, then by the implicit function theorem, we can divide the parameters into two disjoint sets p j and p˜ and obtain, on P h , the relation p j ) β(p˜ ) for some mapping β. Thus, a subset of the parameters can be expressed in terms of the remaining parameters, and the system is not identifiable. Taking the two results together, we have that system D is strongly locally identifiable if and only if ∂S(p)/∂pp)p0 is full rank. 9
6618 Ind. Eng. Chem. Res., Vol. 42, No. 25, 2003
Literature Cited (1) Grewal, M. S.; Glover, K. Identifiability of linear and nonlinear dynamical systems. IEEE Trans. Autom. Control 1976, 21, 833. (2) Vajda, S. Structural identifiability and the realization problem. Sci. Pap. Inst. Tech. Cybernetics Tech. Univ. Wroclaw 1985, 29, 228. (3) 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. (4) McLellan, P. J. A differential-algebraic perspective on nonlinear controller design methodologies. Chem. Eng. Sci. 1994, 49, 1663. (5) 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. (6) Dudukovic, M. P.; Xu, Z. Modeling and simulation of semibatch photo reactive distillation. Chem. Eng. Sci. 1999, 54, 1397. (7) Walter, E.; Pronzato, L. Identifiabilities and nonlinearities in nonlinear systems. Model. Estim. 1995, 1, 111. (8) Reich, S. On the local qualitative behavior of differnetialalgebraic equations. IEEE Trans. Autom. Control 1995, 14, 427. (9) Campbell, S. L. Linearization of daes along trajectories. Z. Angew. Math. Phys. 1995, 46, 70.
(10) Ben-Zvi, A. Identifiability of Differential Algebraic Systems. Ph.D. Dissertation, Queens University, Kingston, ON, Canada, in progress. (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) Yip, E. L.; Sincovec, R. F. Solvability, controllability and observability of continuous descriptor systems. IEEE Trans. Autom. Control 1981, 26, 702. (13) Gantmacher, F. R. Matrix Theory; Chelsea Publishing Company: New York, 1960; Vol. I. (14) Gantmacher, F. R. Matrix Theory; Chelsea Publishing Company: New York, 1960; Vol. I. (15) Hoffman, K.; Kunze, R. Linear Algebra, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 1971. (16) Ben-Zvi, A. Maple sheets for selected examples. Available at http://www.chemeng.queensu.ca/People/PersonalWebs/pjm/, 2003. (17) Rugh, W. J. Linear System Theory; Prentice Hall: Upper Saddle River, NJ, 1993. (18) Brockett, R. W. Finite Dimensional Linear Systems. In Series in Decision and Control; Howard, R. A., Ed.; John Wiley and Sons: New York, 1969.
Received for review April 17, 2003 Revised manuscript received August 8, 2003 Accepted September 26, 2003 IE030317I
